遲緩愛德華氏菌刺激對牙鲆轉錄組差異基因表達的影響

周 密 ，鄭津輝 ，李慶亞 ，耿緒云 ，潘寶平 ，孫金生 ，，高 虹

(1.天津師范大學 生命科學學院，天津 300387；2.天津師范大學 天津市動植物抗性重點實驗室，天津 300387；2.天津市水產養殖病害防治中心，天津 300221)

牙鲆(Paralichthys olivaceus)別稱比目魚，隸屬于牙鲆科(Paralichthyidae)牙鲆屬(Paralichthys)，具有生長快速、體型龐大、肉質鮮美、營養價值高等優點，已成為亞洲海水養殖魚類的代表之一[1].牙鲆養殖過程中，養殖環境、生長速率以及抗病能力等是影響牙鲆產量和質量的重要因素.牙鲆的疾病種類很多，其中細菌性疾病會嚴重危害牙鲆養殖，如鰻弧菌病、愛德華氏綜合病等[2-3].遲緩愛德華氏菌(Edwardsiella tarda)屬于革蘭氏陰性菌，不僅能使較多的經濟魚類患病，給水產養殖業造成嚴重經濟損失，也是一種人、魚共患病原菌，直接威脅人類健康[4].在人工養殖環境條件下，由于養殖密度大、水體環境中代謝物及餌料沉積、抗生素藥物的使用等原因，水體水質差，極易引發疾病，因此選育抗病能力強的牙鲆魚種是推動牙鲆養殖業發展的重要途徑之一[5].利用遺傳學方法育種是目前牙鲆品種改良的重要手段，除了傳統的雜交育種之外，對牙鲆遺傳結構進行分析和改造也為培育抗病性能優良品種提供了理論依據，如Sekino等[6]和Coimbra等[7]分離了牙鲆大量的微衛星DNA，為種群遺傳分析、遺傳圖譜的構建奠定了基礎.高通量轉錄組測序(RNA-seq)是在轉錄水平上獲得大量基因序列信息的一種技術，目前廣泛應用于真核生物的基因測序[8-9].RNA-seq在研究基因功能、基因表達以及新基因的挖掘等方面應用越來越多，成為研究魚類遺傳學的一種重要方法[10-11].如陳阿琴等[12]對不同光周期條件下牙鲆的倒數8節尾椎骨進行轉錄組測序，功能富集分析表明差異基因主要富集在鈣離子信號轉導、光信號轉導等通路上.

本研究利用遲緩愛德華氏菌刺激牙鲆，對其腎臟組織進行轉錄組測序和生物信息學分析，包括基因功能注釋、SNP(單核苷酸多態性)分析、SSR(簡單重復序列)分析等，探討差異基因的表達情況，為研究牙鲆的基因功能、發掘新的免疫基因和克隆提供參考.

1 材料與方法

1.1 材料及處理

實驗用牙鲆購于天津市濱海新區水產養殖場，選擇生長良好的牙鲆6尾(免疫組、對照組各3尾，設置3個平行處理)，體長(10±2)cm.為了消除環境因素對牙鲆生長的影響，選擇可供氧、控溫和水過濾的水循環魚缸，并對魚缸進行消毒.鹵鹽鹽度為1.7%～1.8%(質量分數)，水溫20℃左右，充氧處理7 d.實驗開始前讓牙鲆在調配好的水環境下暫養1 d以適應環境.牙鲆營底棲生活，不適宜強光照射，因此對魚缸進行遮光處理，每天喂食(人工配置飼料)1次.

遲緩愛德華氏菌來自天津市水產養殖病害防治中心，將其接種于pH值為7.2的LB液體培養基中，放入28℃恒溫震蕩培養箱(200 r/min)中培養24 h.當菌液OD600值等于1.0時，在4℃、12 000 r/min條件下離心10 min，收集菌體.感染前用PBS離心洗滌菌體3次，利用血球計數板計數，用PBS將菌液稀釋到1×107CFU/mL，立即注射進免疫組牙鲆體內.

1.2 主要儀器和試劑

1.2.1 儀器

Hiseq 2000測序儀，美國Illumina公司；核酸蛋白檢測儀，美國Thermo公司；PCR儀、凝膠成像儀，美國Bio-Rad公司；電泳儀，北京六一儀器廠.

1.2.2 試劑

Total RNA提取試劑盒、PCR擴增試劑盒，上海生工生物有限公司；cDNA第一鏈合成試劑盒，美國Sangon BiotechR公司；DNA Maker，大連寶生物科技有限公司.

1.3 方法

1.3.1 樣品采集和總RNA抽提

對免疫組牙鲆注射遲緩愛德華氏菌，注射量為1.5mL；對照組牙鲆僅注射PBS緩沖液，注射量也為1.5mL.6 h后解剖魚的腎臟，將其剪切成小塊儲存于冷凍管中，立即放入液氮冷凍，之后存于-80℃冰箱中用于RNA提取.由北京諾禾致源有限公司提取RNA，每組3個樣品所提取的RNA等量混合，用于測序.

1.3.2 建庫測序

建庫測序由北京諾禾致源有限公司完成，基本流程如下：

①瓊脂糖凝膠電泳檢測分析后，選取未被降解、完整性好且沒有被污染的RNA，使用Nanodrop微量核酸定量儀檢測RNA的純度，使用Qubit熒光定量儀測定RNA的濃度.采用Agilent 2100對RNA的質量進行檢測[13-14]，之后分別把每組處理中的3個樣品所提取的RNA進行等量混合，用于后續實驗.

②牙鲆RNA質量檢測合格后，使用帶有Oligo(dT)的磁珠富集 mRNA.然后用 fragmentation buffer將mRNA打斷為200～700 nt(nt表示核苷酸數量)的片段，以其為模板合成第一條cDNA鏈，在反應體系中再加入緩沖液、DNA聚合酶I、RNase H、dNTPs以及隨機引物，合成雙鏈cDNA.用AMPure XP beads純化回收合成的雙鏈cDNA，進行粘性末端修復和加A尾處理，連接上測序接頭后對cDNA片段進行篩選[13，15].PCR擴增后，利用AMPure XP beads純化擴增產物，最后進行文庫構建.利用Qubit2.0、Agilent 2100和QPCR方法對建好的文庫進行質量檢測，質檢合格后，采用Illumina HiseqTM2000完成轉錄組測序[13，16].

1.3.3 序列拼接和組裝

對原始測序數據(raw reads)進行過濾，去除帶有接頭的序列、低質量序列(質量值Qphred≤20的堿基數占整個序列50%以上的序列)和不確定的堿基信息占整條序列10%的序列，得到clean reads.利用轉錄組拼接軟件Trinity對clean reads進行拼接形成轉錄本，利用Corest對拼接得到的轉錄本進行層次聚類，選取最長的Cluster序列進行功能注釋分析[17].

1.3.4 轉錄本功能的注釋

利用生物信息學方法對拼接好的Unigene的基因結構、表達水平和基因功能進行注釋.為了獲取更全面的基因功能信息，進行多個數據庫的比對，包括蛋白數據庫 NR(NCBI non-redundant protein sequences)、NT(NCBI nucleotide sequences)、PFAM(Protein family)、KOG/COG(clusters of orthologous groups of proteins/euKaryotic ortholog groups)、Swiss-Prot(a manually annotated and reviewed protein sequence database)、KEGG(Kyoto encyclopedia of genes and genomes)、GO(gene ontology)[14，18].通常用 blastx 的方式將 Unigene序列分別和NR、Swiss-Prot、KEGG和COG數據庫進行比對，從而得到給定Unigene的蛋白功能注釋信息.與NR數據庫進行比對后，可獲取牙鲆基因序列與其近緣物種的相似性，并且對牙鲆基因進行功能注釋[14]；依據NR等數據庫的注釋信息，對牙鲆基因進行GO功能注釋和分類，可以從宏觀水平上了解牙鲆基因的分布情況[13]；依據COG數據庫能夠對基因產物進行同源性分類[13]；目的基因和KEGG信號通路數據庫比對之后，可以對目的基因所參與的信號通路進行注釋，從而確定目的基因在其信號通路中的位置和功能，為確定基因的相關生物學途徑提供理論依據[13].

1.3.5 差異表達基因(DEGs)的篩選和分析

差異表達基因分析是為了找出免疫組和對照組相比表達發生變化的基因.比較牙鲆對照組與免疫組的基因表達水平，使用RSEM軟件對比對結果進行定量分析，采用FPKM方法對定量分析結果進行基因表達量的估算.采用TMM對read count標準化處理后，利用DESeq算法進行差異表達分析，以q＜0.005且|log2Fold Change|＞1為篩選差異基因功能的標準.篩選出遲緩愛德華氏菌刺激6 h后牙鲆腎臟的差異顯著表達基因，對差異基因進行GO功能富集分析，獲得GO term被差異基因富集的概率，著重分析顯著富集的GO term，對其中的基因總數進行統計.為了研究基因的功能，對每一個組合上調和下調的差異基因進行KEGG通路富集分析，從而對差異基因所參與的生化代謝途徑和信號轉導通路進行注釋.利用下列分析計算公式篩選出差異基因顯著富集的通路(pathway)，并統計出富集通路中差異基因的數量，以P＜0.05視為顯著富集.

式中：N表示所有的通路數目；M表示所有注釋為某一通路的基因數；n表示在N中的差異表達基因數；m表示某一通路中差異表達的基因數.

1.3.6 差異表達基因驗證

生物信息學分析因為組裝和數據處理的問題常會產生誤差，為了驗證本次轉錄組測序的組裝結果和差異基因篩選的準確性，選擇4個轉錄組中提示的4個差異表達基因(IL-21，白介素-21；IL-1，白介素-1；IL-15，白介素-15；TNFα，腫瘤壞死因子-α)，采用熒光定量PCR方法，分析4個基因在牙鲆注射遲緩愛德華氏菌6 h后的相對表達情況.實驗材料見第1.2.1部分，利用Trizol法提取樣品總RNA，電泳檢測RNA的質量，使用Sangon BiotechR公司的M-MuLV第一鏈cDNA合成試劑盒，參照說明書合成cDNA.以牙鲆β-actin基因為內參，設計4個基因的引物，如表1所示.反轉錄cDNA為模板，利用ABI 7500熒光PCR儀做基因表達定量分析.

表1 4個差異表達基因和β-actin基因的實時定量表達引物序列Tab.1 Real-time quantitative expression of four differentially expressed genes and β-actin gene

2 結果與分析

2.1 數據預處理

采用Illumina HiSeqTM測序技術獲得了牙鲆免疫組和對照組2個轉錄組的數據，對原始數據進行統計，對照組獲得了105 228 444條raw reads，過濾掉含有接頭、低質量的reads后獲得了98 735 284條clean reads，共14.81 G，GC含量為48.87%；免疫組共獲得了97 358 862條raw reads，其中 clean reads有90 362 310條，共13.55 G，GC含量為48.59%.

2.2 序列拼接

用Trinity軟件對clean reads進行拼接后，共得到222 980條轉錄本，99 808個基因.其中，轉錄本長度在 200～500 bp間的有 131 634條，500～1 000 bp間的42 702條，1 000～2 000 bp間的 23 275條，＞2 000 bp的有25369條；基因長度在200～500bp間的有21138個，500～1 000 bp間的 30 211個，1 000～2 000 bp間的23 091個，＞2 000 bp的有25 368個.轉錄本總長度為200 234 420 bp，長度范圍為 201～22 935 bp，平均長度為898 bp，中等長度為411 bp，N50(轉錄本長度相加超過總長度50%時的最后一條轉錄本的長度)為1 886 bp，N90(轉錄本長度相加超過總長度90%時的最后一條轉錄本的長度)為319 bp.基因總長度為158 960 029 bp，長度范圍為 201～22 935 bp，平均長度為1593bp，中等長度為966bp，N50(基因長度相加超過總長度50%時的最后一個基因的長度)為2633 bp，N90(基因長度相加超過總長度90%時的最后一個基因的長度)為694 bp.

2.3 基因功能注釋

對組裝得到的牙鲆轉錄組進行功能注釋，結果如表2所示.

表2 基因注釋成功率統計Tab.2 Statistics of gene annotation success rate

由表2可知，有51 348條Unigenes注釋到NT數據庫，注釋率最高(51.44%)；其次是NR數據庫，注釋了 43 839條(43.97%)；PFAM、Swiss-Prot、GO 數據庫的注釋率相近，約為39%；KO和KOG數據庫的注釋率最低.

2.3.1 NR注釋分析

在NR數據庫中有43 892條Unigenes得到注釋，分別對其進行相似性分布、E值分布(E值反映可信度的高低，E值越小代表可信度越高)和物種分布情況的分析，結果如圖1所示.

圖1 NR數據庫比對結果的統計Fig.1 Statistics of NR database comparison results

由圖1(a)可以看出，注釋到的 Unigenes約有66%與NR數據庫中相應蛋白的相似度超過80%，因此轉錄組所測得數據的可信度比較高.圖1(b)的E值分布圖中，NR所注釋到的E值小于1×10-15的Uingenes占94.1%；在圖1(c)的物種分布圖中，注釋到的Unigenes與大黃魚(Larimichthys crocea)的匹配度最高，33.5%的Unigenes比對到大黃魚中的相應蛋白，說明牙鲆與大黃魚的親緣性比較高；其次是深裂眶鋸雀鯛(Stegastes partitus)，牙鲆有28.7%的Unigenes比對到深裂眶鋸雀鯛中；有24.5%的Unigenes比對到其他魚中.

2.3.2 GO分析

對牙鲆轉錄組進行GO分析，結果如圖2所示.牙鲆轉錄組一共有39 390個Unigenes注釋到GO數據庫中，分為生物過程、細胞組分和分子功能3類.由圖2可以看出，在生物過程中注釋到的基因主要分布在代謝過程、單一生物過程、細胞過程及生物調節過程中；在細胞組分中主要分布在細胞和細胞器部分；分子功能中主要分布在粘合與催化活性方面.

2.3.3 KOG功能注釋分析

牙鲆轉錄組中KOG功能注釋的Unigenes共有19 444條，注釋到26個功能類別上，注釋比例為19.48%，如圖3所示.由圖3可以可以看出，大部分Unigenes都被注釋到KOG分類中的信號轉導機制、一般功能預測、翻譯后修飾、蛋白質周轉和分子伴侶、轉錄、細胞內運輸、分泌和水泡運輸、細胞骨架、RNA加工和修飾等功能類別.

圖2 GO分類Fig.2 GO classification

圖3 KOG分類Fig.3 KOG classification

2.3.4 KEGG分類

牙鲆轉錄組在數據庫中進行KO注釋后可進一步得到所在的通路注釋.本轉錄組所獲得的序列中有25 899條能在KO數據庫中得到注釋，參與到232條相應的通路上.將注釋到的基因所參與的KEGG代謝通路進行歸類，大致有5個分支，分別是細胞過程(CP，cellular processes)、環境信息處理(EIP，environmental information processing)、遺傳信息處理(GP，genetic information processing)、代謝(M，metabolism)和有機系統(OS，organismal systems)，每個分支的KEGG分類結果如圖4所示.免疫組牙鲆經過遲緩愛德華氏菌刺激后，同對照組相比，環境信息處理分支下的信號轉導(signal transduction)注釋了4 513個基因，在有機系統分支下的免疫系統(immune system)注釋了1 922個基因，這為愛德華氏綜合征的防治提供了新的依據.

2.4 差異基因的分析

2.4.1 差異基因的篩選

本研究共得到25 241個表達量不同的Unigenes，按照q＜0.005且|log2(foldchange)|＞1的原則篩選后，2 502個Unigenes視為存在顯著差異表達.以對照組為參照，存在差異表達的Unigenes在免疫組中有1 280個上調，1 222個下調.免疫組和對照組均有表達的基因有51 478個，免疫組和對照組特異性表達的基因個數分別是17 235個和16 694個.

2.4.2 差異基因GO富集分析

通過GOseq方法對差異表達基因進行GO富集分析，注釋到的差異表達基因有2 048個.校正后的P＜0.05時，視該功能為富集項，差異基因顯著富集在38個GO terms上，包括6個細胞組分、12個生物學過程和20個分子功能，如圖5所示.對圖5中的差異基因進行GO富集，差異基因主要富集在生物過程(BP)、細胞組分(CC)和分子功能(MF)3大類別上.差異基因的上調和下調情況如圖6所示.在圖6中，上調和下調的GO terms主要集中在生物學過程中的生物調節過程，包括細胞通訊、信號傳導、單細胞信號轉導、細胞組分中的膜的整體成分和膜的固有成分等.

圖4 KEGG分類Fig.4 KEGG classification

圖5 差異基因GO富集Fig.5 Differential gene GO enrichment

圖6 差異基因的表達變化Fig.6 Differential gene expression changes

2.4.3 差異基因的KEEG富集分析

通過Rich factor、q值和富集到此通路上的基因個數來衡量KEGG的富集程度.差異基因共獲得277個KEGG通路富集，其中627個差異基因顯著富集到20條信號通路上.20條信號通路分別是Osteoclast differentiation、NF-kappa B signaling pathway、Measles、Cytokine-cytokine receptor interaction、Influenza A、TNF signaling pathway、Malaria、Hematopoietic cell lineage、Jak-STAT signaling pathway、Leishmaniasis、Hepatitis C、NOD-like receptor signaling pathway、Herpes simplex infection、Asthma、RIG-I-like receptor signaling pathway、Cytosolic DNA-sensing pathway、Adipocytokine signaling pathway、Legionellosis、AGE-RAGE signaling pathway in diabetic complications、Rheumatoid arthritis，其中 8 個與信號通路有關，3個與發育有關，9個與相關疾病有關.下調基因顯著富集到20個KEGG通路，其中多達7條與代謝和發育有關，包括beta-Alanine metabolism、Osteoclast differentiation、Histidine metabolism、Transcriptional misregulation in cancer、Pyruvate metabolism、Phosphatidylinositol signaling system、Apoptosis-multiple species.上調基因顯著富集到20個KEGG通路，有7條與信號通路有關，包括NF-kappa B signaling pathway、TNF signaling pathway、NOD-like receptor signaling pathway、Cytosolic DNA-sensing pathway、Adipocytokine signaling pathway、RIG-I-like receptor signaling pathway、Jak-STAT signaling pathway.

2.4.4 實時定量PCR結果和轉錄組分析結果比較

牙鲆感染遲緩愛德華氏菌6 h后，利用實時熒光定量PCR方法，分析白介素21基因(IL-21)、白介素1基因(IL-1)、白介素 15基因(IL-15)和腫瘤壞死因子-a基因(TNFα)的相對表達情況，結果如表3所示.

表3 轉錄組數據的Real-time PCR驗證Tab.3 Real-time PCR validation of transcriptome data

由表3可知，基因的實時熒光定量PCR變化趨勢與轉錄組結果變化趨勢一致，在免疫刺激6 h后，4個基因都表現為上調趨勢，說明轉錄組測序結果可靠.

3 討論與結論

本研究用遲緩愛德華氏菌對牙鲆進行免疫刺激，采用Illumina HiSeq測序平臺對牙鲆轉錄組進行測序，對照組獲得了98 735 284條clean reads，共14.81 G，GC含量為48.87%；免疫組獲得了90 362 310條clean reads，共13.55 G，GC含量為48.59%.經過對數據進行處理和組裝后，在NR數據庫中進行注釋.所獲得的轉錄組共包含222 980條轉錄本，其中99 808條序列得到了注釋.分析注釋結果發現，有94.1%的Unigenes E值小于1×10-15，同時有66%的注釋Unigenes與NR數據庫中相應蛋白相似度大于80%，說明大部分數據的可信度高，而且與大黃魚的同源性很高，說明本次測序結果的質量可靠，數據覆蓋率也比較廣.

用遲緩愛德華氏菌對牙鲆進行免疫刺激后共篩選出2 502個差異表達基因，其中1 280個上調，1 222個下調.對差異表達基因分別進行GO和KECC富集分析，篩選出顯著富集的GO terms和KEGG通路富集.對其進行分析，可以確定差異表達基因的功能以及它們所參與的信號通路.對KEGG富集分析發現，牙鲆表達量上調基因的顯著富集主要與信號通路有關，遲緩愛德華氏菌刺激后，NF-κB信號通路中TNFα、IL-1β等免疫因子表達量上調，從而誘導多種免疫基因表達，參與炎癥反應；其余信號通路是TNF信號通路、NOD樣受體信號通路、RIG-I樣受體信號通路、Jak-STAT信號通路等，其中Jak-STAT信號通路中IL-6、SHP2(非跨模型蛋白酪氨酸磷酸酶)、原癌基因(Pim-1)表達上調，該通路參與細胞的免疫調節以及分裂分化等生物學過程.根據這些信號通路中所提供的差異基因表達量變化的信息，可以為研究牙鲆應對遲緩愛德華氏菌刺激時信號通路中各基因所起作用以及篩選免疫相關基因提供理論依據[19].

參考文獻：

[1]劉峰，陳松林，王磊，等.不同牙鲆群體遺傳力和育種值分析[J].中國水產科學，2013，20(4)：691-697.LIU F，CHEN S L，WANG L，et al.Analysis of heritability and breeding value of different stocks for Japanese flounder(Paralichthys olivaceus)[J].Journal of Fishery Sciences of China，2013，20(4)：691-697(in Chinese).

[2]吳戀，孫金生，耿緒云，等.牙鲆Toll樣受體1基因全長cDNA的克隆及特征分析[J].安徽農業科學，2013，40(26)：12754-12760.WU L，SUN J S，GENG X Y，et al.Molecular cloning and expression of toll-like receptors 1 cDNA in Japanese Flounder，Paralichthys oli-vaceus[J].Anhui Agricultural Sciences，2013，40(26)：12754-12760(in Chinese).

[3]GAO H，WU L，SUN J S，et al.Molecular characterization and expression analysis of Toll-like receptor 21 cDNA from Paralichthys olivaceus[J].Fish&Shellfish Immunology，2013，35(4)：1138-1145.

[4]敖弟書，吳中明，王玉，等.遲鈍愛德華氏菌感染大鯢的研究[J].四川動物，2010，29(3)：411-414.AO D S，WU Z M，WANG Y，et al.Study on infection in Andrias davidianus caused by Edwardsiella tarda[J].Sichuan Journal of Zoology，2010，29(3)：411-414.

[5]ANSORGE W J.Next-generation DNA sequencing techniques[J].New Biotechnol，2009，25(4)：195-203.

[6]SEKINO M，HRARA M，TANIGUCHI N.Loss of microsatellite and mitochondrial DNA variation in hatchery strains of Japanese flounder Paralichthys olivaceus[J].Aquaculture，2002，213(1)：101-122.

[7]COIMBRA M R M，KOBAYSSHI K，KORETSUGU S，et al.A genetic linkage map of the Japanese flounder，Paralichthys olivaceus[J].Aquachulture，2003，220(1/2/3/4)：203-218.

[8]李少波，傅國輝.轉錄組測序數據中cSNP和表達差異基因的分析方法[J].上海交通大學學報，2014，34(2)：129-133.LI S B，FU G H.RNA-Seq based analysis on cSNP and gene expression level[J].Journal of Shanghai Jiaotong University，2014，34(2)：129-133(in Chinese).

[9]宋文濤，張瀟峮，廖小林，等.牙鲆微衛星標記遺傳連鎖圖譜的構建[J].農業生物技術學報，2011，19(6)：981-987.SONG W T，ZHANG X Q，LIAO X L，et al.Construction of a microsatellite-based genetic linkage map in Japanese flounder(Paralichthys olivaceus)[J].Journal of Agricultural Biotechnology，2011，19(6)：981-987(in Chinese).

[10]袁靜.尼羅羅非魚性腺發育過程的轉錄組學研究及家族生物信息學分析[D].重慶：西南大學，2013.YUN J.Characterization of Gonadal Transcriptomes from Nile tilapia{Oreochromis niloticus)and Bioinformatic Analysis of Fox Gene Family[D].Chongqing：Southwest University，2013(in Chinese).

[11]XU Q Q，LI R G，MONTE M M，et al.Sequence and expression analysis of rainbow trout CXCR2，CXCR3a and CXCR3b aids interpretation of lineage-specific conversion，loss and expansion of these receptors during vertebrate evolution[J].Developmental&Comparative Immunology，2014，45(2)：201-213.

[12]陳阿琴，張影，陳松林，等.不同光周期條件下日本牙鲆尾部神經分泌系統轉錄組分析[J].水產學報，2016，40(6)：833-843.CHEN A Q，ZHANG Y，CHEN S L，et al.The transcriptome sequencing and functional analysis of CNSS in Paralichthys olivaceus treated with different photoperiods[J].Journal of Fisheries Science，2016，40(6)：833-843(in Chinese).

[13]許寶紅.感病草魚脾臟的比較轉錄組分析[D].長沙：湖南農業大學，2012.XU B H.Transcriptome Analysis of Grass Carp Infected GCHV[D].Changsha：Hunan Agricultural University，2012(in Chinese).

[14]北京諾禾致源科技股份有限公司.諾禾致源真核無參轉錄組生物信息分析結題報告[R].北京：北京諾禾致源科技股份有限公司，2013.Beijing Nuowui Source of Science and Technology Co Ltd.A Report on the Analysis of Bioinformatics in Non-transcriptome of Nucleus of[R].Beijing：Beijing Nuowui Source of Science and Technology Co.，Ltd，2013(in Chinese).

[15]COCK P J A，FIELDS C J，GOTO N，et al.The Sanger FASTQ file format for sequences with quality scores，and the Solexa/Illumina FASTQvariants[J].NucleicAcidsResearch，2010，38(6)：1767-1771.

[16]HANSENKD，BRENNERS E，DUDOIT S.Biases in illumina transcriptome sequencing caused by random hexamer priming[J].Nucleic Acids Research，2010，38(12)：e131.

[17]GRABHERR M G，HAAS B J，YASSOUR M，et al.Full-length transcriptome assembly from RNA-Seq data without a reference genome[J].Nature Biotechnology，2011，29(7)：644-652.

[18]FINN R D，TATE J，MISTRY J，et al.The Pfam protein families database[J].Nucleic Acids Research，2008，36(1)：281-288.