全長轉(zhuǎn)錄組有參報賬用結(jié)題報告

合同關(guān)鍵指標

2+3全長有參轉(zhuǎn)錄組服務(wù)合完成1個樣品的全長轉(zhuǎn)錄組，樣品產(chǎn)出不少于8GbCleanData。完成可變剪接分析。完成長鏈非編碼RNA預(yù)測和長鏈非編碼RNA靶預(yù)測。完成轉(zhuǎn)錄本功能注釋分析。實驗流全長轉(zhuǎn)錄組實驗流程包括樣品檢測、文庫構(gòu)建和上機。實驗流程見下圖全長轉(zhuǎn)錄組實驗流程樣品檢高質(zhì)量的RNA是整個項目成功的基礎(chǔ)，為保證數(shù)據(jù)準確性，我們使用以下方法對樣品進行檢測，檢測結(jié)果達到要求后方可進行庫Nanodrop檢測RNA的純度（OD260/280）、濃度、核酸吸收峰是否正常Agilent2100精確檢測RNA的完整性，檢測指標包括:RIN值、28S/18S、圖譜基線有無上抬、5S峰電泳檢測RNA樣品是否有組DNA的污文庫構(gòu)樣品檢測合格后，進行文庫構(gòu)建，主要流程如下使用SMARTer?PCRcDNASynthesisKit合成mRNA的cDNABluePippin篩選全cDNA片段，構(gòu)建不同大小cDNA文庫:。(3)再次PCR擴增放大篩選的全長cDNA。對全長cDNA進行末端修復連接SMRT啞鈴型接頭進行核酸外切酶消化使用BluePippin進行二次篩選，獲得文庫文庫質(zhì)文庫構(gòu)建完成后，對文庫質(zhì)量進行檢測，檢測結(jié)果達到要求后方可進行上機，檢測方法如下使用Qubit2.0進行準確定量使用Agilent2100對文庫大小進行檢測，文庫大小符合預(yù)期后才可進行上機上庫檢合格后，按照目標下機數(shù)據(jù)量，使用PacBio儀器進行全長轉(zhuǎn)錄組生物信息學轉(zhuǎn)錄組研究是理解生命過程必不可少的工具之一，然而基于第二代高通量平臺的R-e2.0技術(shù)往往不能準確得到或組裝出完整轉(zhuǎn)錄本，無法識別Isoform、同源、超、等位表達的轉(zhuǎn)錄本，使人們難以理解這一生命活動更次的含義?；赼cioRT單分子實時技術(shù)的全長轉(zhuǎn)錄組無須打斷A片段，使用AE技術(shù)反轉(zhuǎn)錄得到的全長cDA。該平臺的超長讀?。╡dian10kb）包含了單條完整轉(zhuǎn)錄本序列信息，后期分析無需組裝，所測即所得[1][2][3]。獲取全長轉(zhuǎn)錄組的分析過程主要包括3個階段[4]，全長序列識別、isoform水平聚類得到一致性序列和一致性序列polishin，詳細步驟如下：從原始下機序列中提取ROI（ReadsOfInsert）序列，過濾序列中的cDNA引物和polyA，并根據(jù)序列中是否存在3'引物、5'引物和PolyA（可選）將序列分成全長序列與非全長序列、嵌合序列與非嵌合序列。使用ICE（Iterativeisoform-clustering）算法將來自同一isoform的全長序列進行聚類，將具有相似序列的全長序列聚成一簇（cluster），每個cluster得到一條一致性序列使用Quiver算法聚類非全長序列，對得到的一致性序列進行校正（polishing），篩選高質(zhì)量的序列進行后續(xù)分析?？紤]到cNA建庫的局限性[4]，我們篩選得到的高質(zhì)量序列由于建庫時5'端序列的缺失可能是非全長序列，因此，僅5'末端外顯子存在差異，其余外顯子均一致的序列進行合并，取其中最長的序列作為最終的轉(zhuǎn)錄本序列最終得到的轉(zhuǎn)錄本序列可直接用于后續(xù)的Isoform、同源、、等位、SS、可變剪接、lncNA等分析。引導人們更次的理解位于中心法則中心地位的這一生命活動，另外還可用于對所在組的注釋升級，完善組數(shù)據(jù)庫。轉(zhuǎn)錄組生物信息分析流程見下圖全長轉(zhuǎn)錄組生物信息分析流程2.2.1數(shù)基于單分子（SequencingBySynthesis，SMRT）技術(shù)，Pacbio高通量平臺對cDNA文庫進序，產(chǎn)出大量的數(shù)據(jù)，稱為原始數(shù)據(jù)（RawData）。結(jié)果文件列(1)pacbio文件命名規(guī)則pacbio產(chǎn)出文件或序列命名方式示例如下pacbio命名示"m"=TimeofRunStart(yymmdd_hhmmss)上機時InstrumentSerialNumber儀器序列編SMRTCellBarcodeSMRTcellbarcode編SetNumber(a.k.a."LookNumber".Deprecatedfield,usedinearlierversionofRS)6PartNumber(usually"p0","X0"whenusingexpiredreagents)ZMWholenumberSubreadRegion(start_stopusingpolymerasereadcoordinates)?subread在原始序列中的起始和終止(2)主要文件介紹主要生成*.bas.h5,*.1.bax.h5,*.2.bax.h5,和*.3.bax.h5文件。bax.h5文件包括每個循環(huán)中堿基識別信息。bas.h5文件bax.h5文件的描述（pointer）信息。bas.h5文件說明請考 /software/instrument/2.0.0/bas.h5+Reference+Guide.pdfetadat.xl文件主要包括數(shù)據(jù)的基本信息，包括所使用的酶和化學試劑信息、樣品名稱和其他實驗信息。詳細說明請參考:https: /software/instrument/2.0.0/Metadata+Output+Guide.pdf。sts.xml文件主要包括每輪（movie）中的統(tǒng)計信息說明，詳細說明請考 /software/instrument/1.3.1/Statistics+Output+Guide.pdfsubread序列文件：原始序列（Polymerasereads）從接頭處打斷并過濾掉接頭序列后得到subread，該序列FASTQFASTA格式保存，其中FASTQ主要包括序列和質(zhì)量信息，F(xiàn)ASTQ格式文件示意圖如下:FASTQ格式文件示意注:FASTQ文件中通常每4行對應(yīng)一個序列單元:第一行以@開頭，后面接著序列標識（I）以及其它可選的描述信息;第二行為堿基序列，即es;第三行以開頭，后面接著可選的描述信息;第四行為as每個堿基對應(yīng)的質(zhì)量打分編碼，長度必須和es的序列長度相同，該行中每個字符對應(yīng)的ASII值減去，即為對應(yīng)第二行堿基的質(zhì)量值。中每個ZMW(zero- waveguide)中得到的原始序列如下圖所示序列介Polymerase聚合酶序列（polymeraseread）:DNA聚合酶以SMRTbell?環(huán)狀模板鏈合成的核酸序列，可用于過程中每輪（run）的質(zhì)控。olyerasereds經(jīng)過濾后僅剩余高質(zhì)量片段，包含接頭序列和通過環(huán)狀模板鏈合成的含多個序列的拷貝，如圖6中'ra'所示。每個聚合酶序列（polymeraseread）可以分割成一個或多個子序列（Subread），subread是聚合酶以SMRTbell?一條模板鏈經(jīng)過一輪（passes）合成的，不包括接頭序列。每個subreads包含質(zhì)量值和相關(guān)酶活參數(shù)。圖6共有2個半的subread。其中fullpasses（完整的subread）數(shù)為2。numberoffull指原始序列中存在兩端均含有SMRTbell?接頭（adapter，圖6中黑域）的子序列（接頭間的序列）個數(shù)，圖6中fullpasses2CircularConsensus(CCS)CCS序列是通過每個ZMW孔中的（subreads）子序列得到的一致性序列，無需進行參考序列的比對。不同于ROI序列，序列要求每個插入序列中至少含2個完整（full-pass）的subreadsReadsofInsert根據(jù)需求篩選原始序列中滿足最小numberoffullpasses（<低序列準確性（0.70~1.00）的原始序列，通過子序列獲得一致序列即為每個單分子反應(yīng)器ZMW的ROI序列。ROI是每個單分子反應(yīng)器ZMW里插入序列的最高質(zhì)量序列。ROI序列可以從一定程度上評估建庫質(zhì)量和SMRT?Cell上樣時序列的長度。為充分利用數(shù)據(jù)量，本項目不考慮完整的subread個數(shù)（即設(shè)置最小fullpasses=0），最低序列準確度為0.75。（詳細內(nèi)容請參見：/acificiocince/c_prier/iki/ndersanding-acio-tansripoe-daa#radxplined）。例如，如果根據(jù)模板鏈最終合成了一個半的subread，最終將合并成一條插入序列（ReadofInsert，ROI），而CCS序列要求最少有兩個完整的subreads，所以它是ROI序列的特殊情況。ROI序列可以很好的評估SMRTcell上樣時的插入序列的長度。對于較長的模板鏈，最終產(chǎn)生的ROI序列可能和polymeraseread序列等長。Full-Length(FL)ReadversusNon-Full-Length(nFL)SMRTLink軟件中IsoSeq2定義兩端同時含有3'引物和5'3'引物前polyA尾（可選）的序列稱為全長序列（Full-Length(FL)Read）。5'或3'primer可以是Clontech或其他全長cDNA建庫引物，或特異性的RT-PCR引物。反之，則為非全（non-full-lengthread）Full-Lengthnon-chimericRead(FLNC)：建庫過程中因接頭濃度或SMRTbell濃度過低造成兩個cDNA模板鏈直接相連而生的嵌合序列稱為人工嵌合序列，如下圖所示全長序列中的人工嵌合序SMRT所使用的PacBio儀器每個cell含有ZMWs，reads進入ZMW孔中被，一個ZMW中含一條的reads（P1）為有據(jù)本項目F01下機數(shù)據(jù)過濾PolymeraseRead片段長度小于50bp、序列準確性小于0.75的序列，將剩余序列從接頭處打斷并過濾掉接頭序列后得到subreads，過濾長度小于50bp的sbred，剩余subred即為cleandata。 數(shù)據(jù)統(tǒng)計SampleSamplecDNASMRTData1-22-13-1根據(jù)條件fullpasses>=0且序列準確性大于0.75從原始序列提取ROI序列。統(tǒng)計各建庫中ROI（ReadsofInsert）序列數(shù)、ROI的堿基數(shù)和插入序列的平均長度對下機數(shù)據(jù)進行評估。該項目各建庫數(shù)據(jù)產(chǎn)出統(tǒng)計見下表:ROI數(shù)據(jù)統(tǒng)計部分結(jié)果展cDNAReadsofReadBasesofMeanReadLengthofMeanReadQualityof1-2-3-(1)ROI序列長度分布:ROI序列長度與建庫時的cDNA長度選擇有關(guān)，ROI序列隨cDNA長度的增加而增加cDNA長度建庫ROI序列長度分布如下圖所示cDNA建庫ROI長度分布示意注:橫坐標:ROI長度分布;左側(cè)縱坐標：ROI序列長度頻數(shù)分布直方圖；右側(cè)縱坐標：ROI序列長度累積頻率曲線(3)ROIfullpasses數(shù)分布ROI的各孔fullpassescDNA的長度有關(guān)，一般cDNA長度的增加而減少。ROI序列的準確性fullpasses的影響，fullpasses越高，序列準確性越高，各cDNA建庫數(shù)據(jù)的fullpasses數(shù)分布如下圖所示:生成ROI序列的fullpasses數(shù)分布示注:橫坐標:生成ROI序列的fullpasses數(shù);縱坐標:相應(yīng)fullpasses數(shù)的ROI序列數(shù)(3)ROI質(zhì)量分布 質(zhì)量值反映序列的準確性，各建庫數(shù)據(jù)的ROI序列質(zhì)量分布如下圖所示ROI質(zhì)量值分布示意注:橫坐標:ROI質(zhì)量值分布;左側(cè)縱坐標:ROI序列質(zhì)量頻數(shù)分布直方圖;右側(cè)縱坐標:大于相應(yīng)質(zhì)量值的堿基通過檢測ROI序列中是否包含5'引物，3'引物及polyA尾，可將序列分成全長序列（包含5'引物，3'引物及polyA尾）和非全長序列。classify過程去除ROI序列中cDNA引物序列及polyA序列獲得建庫時的插入序列，同時根據(jù)建庫時兩端引物的差別確定鏈合成方向，并將序列分為全長序列和非全長序列、嵌合序列和非嵌合序列。全長序列數(shù)據(jù)統(tǒng)計cDNAReadsofNumberoffiveprimeNumberofthreeprimeNumberofpoly-A1-2-3-全長序列長度反映了建庫時cDA序列長度，可通過統(tǒng)計全長序列的長度評估建庫質(zhì)量，本次數(shù)據(jù)中的全長序列長度與建庫大小一致。不同cDNA建庫時FLNC序列的長度分布如下圖所示：全長非嵌合序列長度分布示意注:橫坐標:FNC序列長度分布;左側(cè)縱坐標：FNC序列長度頻數(shù)分布直方圖；右側(cè)縱坐標：FLNC序列長度累積頻率曲線。LC長度除受建庫時cDA的長度影響外，還受試劑和時間等因素的影響，故全長序列的百分比（L%）在不同cell中會有很大變化。因建庫過程中SMRTbell或adaptor濃度過低，會產(chǎn)生嵌合（ArtificialConcatemers）序列，本次全長序列中ArtificialConcatemers比例0.38%，可認為該數(shù)SMRTbell濃度適中。各庫中所有ROI序列分類結(jié)果如下圖所示:序列分類（classify）示意isoformSMRTysis(v2.3.0)軟件使用ICE（IterativeClusteringforErrorCorrection）算法對序列進行迭代聚類，將相似的序列（即同一轉(zhuǎn)錄本的多個拷貝）聚類到一簇cluster，每個cluster得到一個一致性序列（consensusisoform）。結(jié)合非全長序列，使用quiver程序?qū)Ω鱟luster中的一致序列進行校正,最終得到準確度大于99%的高質(zhì)量轉(zhuǎn)錄（HQ，high- isoforms）。不同建庫長度得到的HQ和LQ(低質(zhì)量轉(zhuǎn)錄本，low-quality)轉(zhuǎn)錄本統(tǒng)計結(jié)果見下表ICE聚類結(jié)果統(tǒng)計NumberofconsensusAverageconsensusisoformsreadNumberofpolishedhigh-quality本項目F01得到consensusisoform65,257個，其中HQ轉(zhuǎn)錄本數(shù)為24,889，LQ轉(zhuǎn)錄本40,345個，得到的consensus序列長度分布如下所示consensusisoform長度分布示注:橫坐標:consensusisoform序列長度分布;左側(cè)縱坐標：consensusisoform序列長度頻數(shù)分布直方圖；右側(cè)縱坐consensusisoform序列長度累積頻率曲線利用proovread軟件[5]通過RNA-Seq數(shù)據(jù)對各個樣品得到的低質(zhì)量一致性序列進行校正，提高序列的準確性。轉(zhuǎn)錄本去冗在全長轉(zhuǎn)錄本cluster過程中參數(shù)設(shè)置較嚴格，為得到質(zhì)量較高的一致性序列，同一轉(zhuǎn)錄本的多拷貝序列分到不同cluster的可能性比隨機將兩個不屬于同一轉(zhuǎn)錄本拷貝的序列分到同一cluster中的可能性大，不可避免的產(chǎn)生了冗余序列。同時，全長轉(zhuǎn)錄本過程中，3'端因存在poly-A結(jié)構(gòu)，可以確定3'端比較完整，而5'端序列可能存在降解，導致同一轉(zhuǎn)錄本的不同拷貝分到不同的cluster中，如下圖綠色圈中所示，5'端差異造成不同轉(zhuǎn)錄本，導致冗余序列的產(chǎn)生。5'序列差異造成轉(zhuǎn)錄本序列通過GMAP（GenomicMapandAlignmentProgram）[6]將得到的校正后的一致序列與參考組進行序列比對（設(shè)置參數(shù)--cross-species--allow-close-indels0），使用cDNA_Cupcake（ /Magdoll/cDNA_Cupcake/wiki）軟件對比對結(jié)果去冗余，過濾identity小于0.9，coverage小于0.85的序列，合并僅5’端外顯子有差異的比對。過濾的轉(zhuǎn)錄本列1.轉(zhuǎn)錄本序列1.轉(zhuǎn)錄本分類結(jié)1.轉(zhuǎn)錄本注釋結(jié)果1.轉(zhuǎn)錄組eads與參考組序列比對結(jié)果文件（通常為BAM格式）、物種參考組序列和注釋文件，推薦使用整合組瀏覽器（IGV，IntegrativeGenomicsViewer）進行可視化瀏覽。IGV具有以下特點:能在不同尺度下顯示單個或多個eads在參考組上的位置，包括eads在各個上的分布情況和在注釋的外顯子、內(nèi)含子、剪接接合區(qū)、間區(qū)的分布情況等;能在不同尺度下顯示不同區(qū)域的Reads豐度，以反映不同區(qū)域的轉(zhuǎn)錄水平能顯示及其剪接異構(gòu)體的注釋信息能顯示其他注釋信息既可以從服務(wù)器端各種注釋信息，又可以從本地加載注釋信息IGV瀏覽器界利用BUSCO[7]對去冗余后的轉(zhuǎn)錄組進行完整性評估，評估結(jié)果如圖所示轉(zhuǎn)錄組完整性評估結(jié)果示意轉(zhuǎn)錄本融合發(fā)因為去冗余分析過程中，通過設(shè)置最小coverage為0.85過濾了嵌合序列，故融合轉(zhuǎn)錄本分析中需要使用去冗余前的一致序列，按下條件進行融合轉(zhuǎn)錄本篩選比對到2個或多個位點每個位點必須比對上至少5%的轉(zhuǎn)錄本長度,最小比對長度為所有位點比對到的總長度必須占轉(zhuǎn)錄本總長度的95%以上兩位點間距離必須達10kb以上。gff注釋文件:融合轉(zhuǎn)錄本GFF文件1.提供組注釋補充信息的同時，也提供以FATA格式的新轉(zhuǎn)錄本序列。融合轉(zhuǎn)錄本序列的FSTA文件:融合轉(zhuǎn)錄本FASTA文件1.轉(zhuǎn)錄本結(jié)構(gòu)分轉(zhuǎn)錄生成的前體mN（pre-RA），有多種剪接方式，選擇不同的外顯子，產(chǎn)生不同的成熟mN，從而翻譯為不同的蛋白質(zhì)，構(gòu)成生物性狀的多樣性。這種轉(zhuǎn)錄后的mNA加工過程稱為可變剪接或選擇性剪接（lternativespliing）。通過Astalavista軟件[8]獲取每個樣品存在的可變剪接類型，主要的可變剪接類型如下圖所示可變剪接類注：(A外顯子跳躍；(B)可變轉(zhuǎn)錄終止位點；(C)可變外顯子；(D)可變轉(zhuǎn)錄起始位點；(E內(nèi)含子保留我們從Astalavista軟件分析結(jié)果中,對轉(zhuǎn)錄本發(fā)生上述5種可變剪接事件情況進行統(tǒng)計,各樣品中預(yù)測的可變剪接事件數(shù)量統(tǒng)計見下可變剪接事件數(shù)量統(tǒng)注：Alternative3'splicesite:可變轉(zhuǎn)錄終止位點；Alternative5'splicesite:可變轉(zhuǎn)錄起始位點；Exonskip外顯子跳躍；Intronretention內(nèi)含子保留；Mutuallyexclusiveexon:可變外顯子可變剪接事件統(tǒng)計1.由于使用的軟件或數(shù)據(jù)本身的局限性，導致所選參考組的注釋往往不夠精確，這樣就有必要對原有注釋的結(jié)構(gòu)進行優(yōu)化。如果在原有邊界之外的區(qū)域有MappedReads支持，將的非翻譯區(qū)（UntranslatedRegion，UTR）向上下游延伸，修正的邊界。結(jié)構(gòu)優(yōu)化結(jié)果見下面文件：結(jié)構(gòu)優(yōu)化結(jié)注：GeneID：ID；Locus：座，格式為“編號:起點坐標-終點坐標”；Strand：正負鏈；Site：優(yōu)化的位置，3'或5'UTR；OriginalSite：原來注多聚腺苷酸化是指多聚腺苷酸與信使RNA（mRNA）分子的共價鏈結(jié)。在蛋白質(zhì)生物合成的過程中，這是產(chǎn)生準備作翻譯的成mRNA的方式的一部份。在真核生物中，多聚腺苷酸化是一種機制，令mRNA分子于它們的3'端中斷。多聚腺苷酸尾（A尾）保護mRNA，免受核酸外切酶，并且對轉(zhuǎn)錄終結(jié)、將mRNA從細胞核輸出及進行翻譯都十分重要。在原核生物中，前體mRNA的可變多聚腺苷酸化（alternativepolyadenylation,APA）可能貢獻于轉(zhuǎn)錄組多樣性，組的編碼能力以及的調(diào)控機制。我們采用TAPISpipeline[9]來識APA。各個樣品所識別出的APA如下：可變多聚腺苷酸化結(jié)1.多聚腺苷酸化位點個數(shù)分布示意注:橫坐標:多聚腺苷酸化位點個數(shù);縱坐標：個數(shù)利用MEME對所有轉(zhuǎn)錄本polyA位點上游50bp的序列進行分析，鑒定得到的motif如下圖所示polyA位點上游motifSSR分MISA(MIcroSAliteidentificationtool)是一款鑒定簡單重復序列的軟件，其參考見附表。它可以通過對轉(zhuǎn)錄本序列的分析，鑒定出7種類型的SSR:Mono-nucleotide（單堿基）、Di-nucleotide（雙堿基）、Tri-nucleotide（三堿基）、Tetra-nucleotide（四堿基）、enta-nuceotde（五堿基）、exa-nuclotie（六堿基）、copoundR（混合微，兩個SR距離小于100bp）從新轉(zhuǎn)錄本中篩選500bp以上的轉(zhuǎn)錄本，利用MISA軟件做SSR分析，結(jié)果見下表SSR分析結(jié)果統(tǒng)SearchingTotalnumberofsequencesTotalsizeofexaminedsequencesTotalnumberofidentifiedNumberofSSRcontainingNumberofsequencescontainingmorethan1NumberofSSRspresentincompoundMonoDiTriTetraPentaHexaSSR分析結(jié)注：每一行為一個SSR標記及其引物設(shè)計結(jié)果，每個SSR標記最多設(shè)計三對引物;eI：轉(zhuǎn)錄本編號；SSn：同一轉(zhuǎn)錄本上的SSR序號；SSty：SSR類型，包括完美單堿基重復（）、完美雙堿基重復（p）、完美三堿基重復（p）、完美四堿基重復（p）、完美五堿基重復（）、完美六堿基重復（6）和混合SS（c，即包含至少兩個完美SS，且之間距離小于1b）；SSR：SSR序列，括號內(nèi)為重復單元，括號外數(shù)字表示重復次數(shù)；Siz：SSR的長度；SSStat：SSR在轉(zhuǎn)錄本上的開始位置；SSREd：SSR在轉(zhuǎn)錄本上的結(jié)束位置;FP1'-：第一條正向引物序列；Tm：第一條正向引物序列的退火溫度，單位為C；Siz：第一條正向引物序列的長度；P1'-')：第一條反向引物序列；Tm：第一條反向引物序列的退火溫度，單位為C；Siz：第一條反向引物序列的長度；Psiz：產(chǎn)物的長度；PStt：產(chǎn)物在上的開始位置；PEd：產(chǎn)物在上的結(jié)束位置。對不同SSR類型的密度分布進行統(tǒng)計，結(jié)果見下圖SSR類型分布示意新編碼區(qū)序列預(yù)TransDecoder[10](v3.0.0)軟件基于開放閱讀框（OpenReadingFrame，ORF）長度、對數(shù)似然函數(shù)值（Log-likelihoodcore）、氨基酸序列與fam數(shù)據(jù)庫蛋白質(zhì)結(jié)構(gòu)域序列的比對等信息，能夠從轉(zhuǎn)錄本序列中識別可靠的潛在編碼區(qū)序列（odingequence，CDS），是常用的CDS預(yù)測軟件對可變剪接分析中得到的新轉(zhuǎn)錄本使用TransDecoder軟件對其編碼區(qū)序列及其對應(yīng)氨基酸序列的預(yù)測。本次共獲得orf個，其中完整orf20,456條，預(yù)測的CDS文件結(jié)果如下圖CDS分析結(jié)果文件示意注:文件為標準的FASTA格式，每個序列單元以開始到下一個“<”之前結(jié)束?！?lt;”后面緊接編碼區(qū)序列編號（由轉(zhuǎn)錄本編號和開放閱讀框序號構(gòu)成），轉(zhuǎn)錄本編號，之后是序列描述信息，包括序列類型（OFty）、序列長度（l）、序列在轉(zhuǎn)錄本上的位置（lc）和鏈類型。預(yù)測的編碼區(qū)序列類型有完整（cmlt，即同時預(yù)測到起始子和終止子）、’端部分（im_til，即僅預(yù)測到起始子）、’端部分（pim_atil，即僅預(yù)測到終止子）和內(nèi)部區(qū)段（itrl，即起始子和終止子都沒有預(yù)測到）。從第二行開始直到出現(xiàn)之前為止為編碼區(qū)序列。預(yù)測得到的完整ORF區(qū)編碼蛋白序列長度分布如下預(yù)測的CDS編碼蛋白長度分布示意LncRNA預(yù)因lncA不編碼蛋白，因此，通過對轉(zhuǎn)錄本進行編碼潛能篩選，判斷其是否具有編碼潛能，從而可以判定該轉(zhuǎn)錄本是否為lncRNA。百邁客綜合目前應(yīng)用最廣泛的編碼潛能分析方法對新發(fā)現(xiàn)的轉(zhuǎn)錄本進行l(wèi)ncRNA的預(yù)測，主要包括:CPC[11]分析、CNCI[12]分析、pfam蛋白結(jié)構(gòu)域分析、CPAT[13]分析四種方法。CPC(CodingPotentialCalculator)是一種基于序列比對的蛋白質(zhì)編碼潛能計算工具。通過將轉(zhuǎn)錄本與已知蛋白數(shù)據(jù)庫比對根據(jù)轉(zhuǎn)錄本各個編碼框的生物學序列特征評估其編碼潛能。Score<0時，為noncodingRNA。CPC分析結(jié)果CPC分析結(jié)果統(tǒng)注:第1列:transcript_id:轉(zhuǎn)錄本ID;第2列 length:ORF長度;第3列:type:轉(zhuǎn)錄本類型;第4列:score:轉(zhuǎn)錄本得分，當score＜0時，為NoncodingCNCICNCI(Coding-Non-CodingIndex)分析是一種通過相鄰核苷酸三聯(lián)體特征區(qū)分編碼-非編碼轉(zhuǎn)錄本的方法。該工具不依賴于已知的注釋文件，可以有效對不完整的轉(zhuǎn)錄本和反義轉(zhuǎn)錄本進行預(yù)測。NI工具提供兩種比對模式：ve（脊椎物種）；pl（植物物種），本項目選擇pl參數(shù)。當score<0時，為 CI分析結(jié)果如下:CNCI分析結(jié)果統(tǒng)注:第1列:transcript_id:轉(zhuǎn)錄本ID;第2列:type:轉(zhuǎn)錄本類型;第3列:score:轉(zhuǎn)錄本得分，當score<0時，為Noncoding;第4列:start:轉(zhuǎn)錄本起始位置;第5列:end:轉(zhuǎn)錄本終止位置。CPATCPAT(CodingPotentialAssessmentTool)分析是一種通過構(gòu)建邏輯回歸模型，基于ORF長度、ORF覆蓋度，計算Fickett得分和Hexamer得分來判斷轉(zhuǎn)錄本編碼和非編碼能力的分析方法。CPAT分析結(jié)果如下:CPAT分析結(jié)果統(tǒng)fam數(shù)據(jù)庫是最全面的蛋白結(jié)構(gòu)域注釋的分類系統(tǒng)。蛋白質(zhì)是由一個或多個結(jié)構(gòu)域組成的，而每個特定結(jié)構(gòu)域的蛋白序列具有一定保守性。fam將蛋白質(zhì)的結(jié)構(gòu)域分為不同的蛋白，通過蛋白序列的比對建立了每個的氨基酸序列的HM統(tǒng)計模型。百邁客將轉(zhuǎn)錄本各個編碼框上的蛋白序列與pfam數(shù)據(jù)庫做hmscan同源搜索，能比對上的轉(zhuǎn)錄本即為具有某個蛋白結(jié)構(gòu)域的轉(zhuǎn)錄本，即認為具有編碼能力，而無比對結(jié)果的轉(zhuǎn)錄本被認為是潛在的lncRNA。結(jié)果如下:pfam分析結(jié)果統(tǒng)注:第1列:transcript_id:轉(zhuǎn)錄本ID;第2列hmmacc:比對到pfam結(jié)構(gòu)域ID;:hmmstart:比對到結(jié)構(gòu)域的起始位置;第3列:hmmname:pfam結(jié)構(gòu)稱;第4列:hmmstart:比對到結(jié)構(gòu)域的起始位置;5列:hmmend:比對到結(jié)構(gòu)域的終止位置;6列:hmmlength:pfam結(jié)構(gòu)域的長度;7列bitscore比對打分值;第八列:E-value:比對的E值，pfam結(jié)構(gòu)域篩選的條件E-value<0.001。為直觀展示分析結(jié)果，將以上4種分析軟件鑒定得到的noncdingtranscripts進上4種分析結(jié)果取交集，用于后續(xù)lncRNA4種篩選方法維恩根據(jù)lncRNA在參考組注釋信息（gff）上的位置，對lncRNA進行分類繪圖，結(jié)果如下lncRNA位置分類示lncRNA靶預(yù)對預(yù)測得到的lncRNA序列進行靶預(yù)測?；趌ncRNA與其靶的作用方式，我們采用2種預(yù)測方法第一種，lncRNA調(diào)控其鄰近的表達，主要根據(jù)lncRNA與mRNA的位置關(guān)系預(yù)測，定義中每100kbp范圍內(nèi)存在差異表達lncRNA與差異表達的mRNA;第二種，lncA與mNA由于堿基互補配對而產(chǎn)生作用，主要利用Lncar[14]靶預(yù)測工具對我們的LncRA進行靶預(yù)測。兩種方法共預(yù)測得到的靶結(jié)果請見文件：靶預(yù)測結(jié)果注:第一列:lncRNA的id號;第二列:對應(yīng)的編號lncRNA；靶id號;基于位置關(guān)系的LncRNA靶預(yù)基于物理位置的靶預(yù)測結(jié)注:第一列:lncRNA的id號;第二列:lncRNA對應(yīng)的靶id號;編號基于互補序列的靶預(yù)注:#LncRNA_ID:lncRNA的id號;TargetGene_id:lncRNA對應(yīng)的靶ID轉(zhuǎn)錄因子分轉(zhuǎn)錄因子（Transcriptionfactor）是指能夠結(jié)合在某上游特異核苷酸序列上的蛋白質(zhì)，這些蛋白質(zhì)可以調(diào)控RNA聚合酶與DNA模板的結(jié)合，從而調(diào)控的轉(zhuǎn)錄。植物轉(zhuǎn)錄因子預(yù)測使用iTAK[15]軟件，動物轉(zhuǎn)錄因子鑒定使用動物轉(zhuǎn)錄因子數(shù)據(jù)庫—— 2.0[16]，共預(yù)測得到轉(zhuǎn)錄因子4,005個，預(yù)測結(jié)果如下:轉(zhuǎn)錄因子序列文使用BLAST[17]軟件(version2.2.26)將得到的新轉(zhuǎn)錄本序列與NR[18]、Swissprot[19]，GO[20]、COG[21]、KOG[22]、Pfam[23]、NR數(shù)據(jù)庫是NCBI中的非冗余蛋白質(zhì)數(shù)據(jù)庫，包含了Swissprot、PIR(ProteinInformationResource)、PRF(ProteinResearchFoundation)、PDB(ProteinDataBank)蛋白質(zhì)數(shù)據(jù)庫及從GenBank和RefSeq的CDS數(shù)據(jù)翻譯過來的蛋白質(zhì)數(shù)據(jù)信息。該數(shù)據(jù)庫網(wǎng)址見附表。issprot數(shù)據(jù)庫是由I（歐洲生物信息學）負責的數(shù)據(jù)庫，包含了有相關(guān)參考文獻且經(jīng)過校對的蛋白質(zhì)注釋信息數(shù)據(jù)庫，可信度很高。COG(ClustersofOrthologousGroups)數(shù)據(jù)庫是對產(chǎn)物進行同源分類的數(shù)據(jù)庫，是一個較早的識別直系同源的數(shù)據(jù)庫，通過對多種生物的蛋白質(zhì)序列大量比較而來的。KOG（euKaryoticOrthologGroups）數(shù)據(jù)庫是針對真核生物，基于直系同源關(guān)系，結(jié)合進化關(guān)系將來自不同物種的同源分為不同的Ortholog簇。來自同一Ortholog的具有相同的功能，這樣就可以將功能注釋直接繼承給同一KOG簇的其他成員。Pfam(Proteinfamily)是最全面的蛋白結(jié)構(gòu)域注釋的分類系統(tǒng)。蛋白質(zhì)是由一個個結(jié)構(gòu)域組成的，而每個特定結(jié)構(gòu)域的蛋白序列具有一定保守性。Pfam將蛋白質(zhì)的結(jié)構(gòu)域分為不同的蛋白，通過蛋白序列的比對建立了每個的氨基酸序列的HMM統(tǒng)計模型。GO(GeneOntology)數(shù)據(jù)庫是一個國際標準化的功能分類體系，提供了一套動態(tài)更新的標準詞匯表來全面描述生物體中和產(chǎn)物的功能屬性。該數(shù)據(jù)庫總共有三大類，分別是分子功能(molecularfunction)，細胞組分(cellularcomponent)和生物學過程(biologicalprocess)，各自描述了產(chǎn)物可能行使的分子功能，以及所處的細胞環(huán)境和參與的生物學過程。GO數(shù)據(jù)庫中最基本的概念是Term，每個條目都有一個Term名，比如“cell”、“fibroblastgrowthfactorreceptorbinding”或者“signaltransduction”，同時有一個唯一的編號，形如GO:nnnnnnn。KEGG(KyotoEncyclopediaofGenesandGenomes)數(shù)據(jù)庫是系統(tǒng)分析產(chǎn)物在細胞中的代謝途徑以及這些產(chǎn)物功能的數(shù)據(jù)庫。它整合了組、化學分子和生化系統(tǒng)等方面的數(shù)據(jù)，包括代謝通路(AHWA)、藥物(DU)、疾病(IES)、序列(GENES)及組(GENOME)等。利用該數(shù)據(jù)庫有助于把及表達信息作為一個整體的網(wǎng)絡(luò)進行研究。將可變剪接分析中得到的新轉(zhuǎn)錄本進行功能注釋，各數(shù)據(jù)庫注釋到的轉(zhuǎn)錄本數(shù)量統(tǒng)計見下表注釋的轉(zhuǎn)錄本數(shù)量統(tǒng)計AnnotatedNewIsoformSwiss-(genefamily)是來源于同一個祖先，由一個通過重復而產(chǎn)生兩個或的拷貝而構(gòu)成的一組，它們在結(jié)構(gòu)和功能上具有明顯的相似性，編碼相似的蛋白質(zhì)產(chǎn)物。通過Pfam注釋，轉(zhuǎn)錄本的 分類結(jié)果如下:分類結(jié)注:第一列:名稱;第二列:對應(yīng)轉(zhuǎn)錄本的idSNP/InDel分SNP（SingleNucleotidePolymorphisms）是指在組上由單個核苷酸變異形成的遺傳標記，其數(shù)量很多，多態(tài)性豐富。百邁客基于各樣品reads與參考組序列的opat2比對結(jié)果，使用AK軟件[25]識別樣品與參考組間的單堿基錯配，識別潛在的SNP位點。進而可以分析這些SNP位點是否影響了的表達水平或者蛋白產(chǎn)物的種類Inel(inertin-eleton)是指相對于參考組，樣本中發(fā)生的小片段的插入缺失，該插入缺失可能含一個或多個堿基。GAK也能夠檢測樣品的插入缺失（Inel）。Inel變異一般比SP變異少，同樣反映了樣品與參考組之間的差異，并且編碼區(qū)的Inel會引起移碼突變，導致功能上的變化。GATK識別標準如下：35bp范圍內(nèi)連續(xù)出現(xiàn)的單堿基錯配不超過3個經(jīng)過序列深度標準化的SNP質(zhì)量值大于2.0。各樣品分別按照以上條件篩選，最終獲得可靠的SNP位點SnpEff[26]是一款用于注釋變異（SNP、InDel）和預(yù)測變異影響的軟件。根據(jù)變異位點在參考組上的位置以及參考組上的位置信息，可以得到變異位點在組發(fā)生的區(qū)域（間區(qū)、區(qū)或DS區(qū)等），以及變異產(chǎn)生的影響（同義非同義突變等）。由于轉(zhuǎn)錄完成之后，RA除了需要加帽、加loy()和可變剪接之外，較少mNA會經(jīng)歷RA編輯（RAeditin），從而會產(chǎn)生單堿基的替換、插入、缺失。NA編輯能使同一產(chǎn)生序列多樣的mRA，但是這種多態(tài)性不是組固有的多態(tài)性。從比對結(jié)果來看，SNP和單堿基替換的RNA編輯結(jié)果是一樣的。因此，通過轉(zhuǎn)錄組數(shù)據(jù)識別出SNP不免會含有RNA編輯的產(chǎn)物。SNP位點信InDel位點信注：Chr：SNP/InDel位點所在 編號；Pos：SNP/InDel位點在 上的位置；Gene_id：SNP/InDel位點所在的或原來未注釋的區(qū)（表中用Intergenic表示）；Ref：所選參考 組中的SNP/InDel等位；Alt： 樣品中識別到的其他的SNP/InDel等位；T*：樣品T*該SNP/InDel位點的分型；Depth：樣品T*該SNP/InDel位點的 深度；AlleDp：樣品T*該SNP/InDel位點的各等位 深度；Effect：SNP/InDel所在區(qū)域或類型；Codon_change：編碼改變方式，未改變用點表示。Effect具體說明詳見： 根據(jù)SP位點堿基替換的不同方式，可以將SP位點分為轉(zhuǎn)換（ransition）和顛換（ransvesion）兩種類型。根據(jù)SP位點的等位（Allele）數(shù)目，可以將SNP位點分為純合型SNP位點（只有一個等位）和雜合型SNP位點（兩個或多個等位）。不同物種雜合型SNP所占的比例存在差異。對各樣品篩選出的SNP位點數(shù)目、轉(zhuǎn)換類型比例、顛換類型比例以及雜合型SNP位點比例進行統(tǒng)計，如下表SNP位點統(tǒng)計部分結(jié)果展BMK-SNPGenicIntergenicSNP突變類型統(tǒng)計分布如下圖所示SNP突變類型分布示意注：橫軸為SNP突變類型，縱軸為相應(yīng)的SNP將每個的SP位點數(shù)目除以的長度，得到每個的SP位點密度值，統(tǒng)計所有的NP位點密度值并做密度分布圖。部分的SNP位點密度分布圖如下：SNP密度分布注：橫軸為上平均每1000bp序列中分布的SNP數(shù)目，縱軸為數(shù)采用SNPEff分別對SNP，InDel注釋，SNP，InDel的注釋結(jié)果統(tǒng)計如下所示注：縱軸為SNP所在區(qū)域或類型，橫軸為分類數(shù)目。InDel注釋分類示注：縱軸為InDel所在區(qū)域或類型，橫軸為分類數(shù)目本項目使用指定的組作為參考進行序列比對及后續(xù)分析，利用HISAT2將CleanReads與參考組進行序列比對，獲取在參考組或上的位置信息。將比對到不同上的Reads進行位置分布統(tǒng)計，繪制MappedReads在所選參考組上的覆蓋深度分布圖二代數(shù)據(jù)與參考組比對部分結(jié)果展TotalMappedUniquelymappedMultiplemappedMappedReads在參考組上的位置及覆蓋深度分布示意注：橫坐標為位置；縱坐標為覆蓋深度以2為底的對數(shù)值，以10kb作為區(qū)間單位長度，劃分成多個小窗口（Window），統(tǒng)計落在各個窗口內(nèi)的Mapped Reads作為其覆蓋深度。藍色為正鏈，綠色為負鏈。統(tǒng)計MappedReads在指定的參考組不同區(qū)域（外顯子、內(nèi)含子和間區(qū)）的數(shù)目，繪制組不同區(qū)域上各樣品Mapped的分布圖，如下組不同區(qū)域Reads分布示意注：圖中將組分為外顯子區(qū)、間區(qū)、內(nèi)含子區(qū)，區(qū)域大小按Map到相應(yīng)區(qū)域的Reads在所有MappedReads中所占理論上，來自成熟RA的eads應(yīng)比對到外顯子區(qū)。eads比對到內(nèi)含子是由于RA前體和發(fā)生可變剪切的內(nèi)含子保留；Reads比對到間區(qū)是由于組注釋不完善。合格的轉(zhuǎn)錄組文庫是轉(zhuǎn)錄組的必要條件，為確保文庫的質(zhì)量，從以下3個不同角度對轉(zhuǎn)錄組文庫進行質(zhì)量評估。(1)片段化隨機性檢驗mNA片段化后的插入片段大小選擇，是從RA序列中獨立隨機地抽取子序列，mNA數(shù)目越大、打斷方式和時間控制得越合適，目的RA每個部分被抽取到的可能性就越接近，RA片段化隨機性越高，mRA上覆蓋的eads越均勻。通過apped eds在各mRNA轉(zhuǎn)錄本上的位置分布，模擬mRNA片段化結(jié)果，檢驗mRNA片段化的隨機程度。如果mRNA存在嚴重降解，被降解的堿基序列不能被，即無eads比對上。因此，通過查看appedeads在NA轉(zhuǎn)錄本上的位置分布可了解RA的降解情況。樣品MappedReads在mRNA轉(zhuǎn)錄本上的位置分布如下圖：MappedReadsmRNA上的位置分布示注：橫坐標為標準化后的mRNA位置，縱坐標為對應(yīng)位置區(qū)間內(nèi)Reads在總MappedReads中所占百分比。由于參考的mRNA長度不同，作圖時對把每個mRNA按照長度劃分成100個區(qū)間，進而統(tǒng)計每一區(qū)間內(nèi)的MappedReads數(shù)目及所占的比例，圖中反映的是所有mRNA各個區(qū)間內(nèi)的MappedReads比例的匯總。插入片段長度檢驗插入片段長度檢驗插入片段長度的離散程度能直接反映出文庫過程中磁珠純化的效果。通過插入片段兩端的eads在參考組上的比對起止點之間的距離計算插入片段長度。大部分的真核生物為斷裂，外顯子被內(nèi)含子隔斷，而轉(zhuǎn)錄組得到的是無內(nèi)含子的成熟RN。當RA中跨內(nèi)含子的片段兩端的eads比對到組上時，比對起止點之間的距離要大于插入片段長度。因此，在插入片段長度模擬分布圖中，主峰右側(cè)形成1個或多個雜峰。各樣品的插入片段長度模擬分布圖如下：插入片段長度模擬分布示意注：橫坐標為雙端Reads在參考組上的比對起止點之間的距離，范圍為0到800bp；縱坐標為比對起止點之間不同距離的雙端或插入片段數(shù)量轉(zhuǎn)錄組數(shù)據(jù)飽和度檢驗為了評估數(shù)據(jù)是否充足并滿足后續(xù)分析，對得到的轉(zhuǎn)錄本數(shù)進行飽和度檢測。由于一個物種的數(shù)目是有限的，且轉(zhuǎn)錄具有時間和空間特異性，因此隨著量的增加，檢測到的轉(zhuǎn)錄本數(shù)目會趨于飽和。對于表達量越高的轉(zhuǎn)錄本，越容易被檢測定量。因此，對于表達量越低的轉(zhuǎn)錄本，需要更大的數(shù)據(jù)量才能被準確定量。使用各樣品的appedata對檢測到的不同表達情況的轉(zhuǎn)錄本數(shù)目飽和情況進行模擬，繪制曲線圖如下：轉(zhuǎn)錄組數(shù)據(jù)飽和度模擬示意注：本圖為隨機抽取10%、20%、30%……90%的總體數(shù)據(jù)單獨進行定量分析的結(jié)果；橫坐標代表抽取數(shù)據(jù)定位到組上Reads數(shù)占總定位的reads數(shù)的百分比，縱坐標代表所有抽樣結(jié)果中表達量差距小于15%的Gene在各個FPKM范圍的百分比轉(zhuǎn)錄組可以模擬成一個隨機抽樣的過程，即從一個樣品轉(zhuǎn)錄組的任意一段核酸序列上獨立地隨機抽取序列片段。抽取自某一（或轉(zhuǎn)錄本）的片段數(shù)目服從負二項分布（BetaNegativeBinomialDistribution）。基于該數(shù)學模型，使用Cufflinks軟件[27]的Cuffquant和Cuffnorm組件，通過MappedReads在上的位置信息，對轉(zhuǎn)錄本的表達水平進行定量。抽取自一個轉(zhuǎn)錄本的片段數(shù)目與數(shù)據(jù)（或apped ata）量、轉(zhuǎn)錄本長度、轉(zhuǎn)錄本表達水平都有關(guān)，為了讓片段數(shù)目能反映轉(zhuǎn)錄本表達水平，需要對樣品中的appedeads的數(shù)目和轉(zhuǎn)錄本長度進行歸一化。uffuant和uffnorm采用FPKM（FragmentsPerKilobaseoftranscriptperMillionfragmentsmapped）作為衡量轉(zhuǎn)錄本或表達水平的指標，F(xiàn)PKM計算公式如下：FPKM計算公注：公式中，cDNAFragments表示比對到某一轉(zhuǎn)錄本上的片段數(shù)目，即雙端Reads數(shù)目；MappedFragmentsMillions)表示比對到轉(zhuǎn)錄本上的片段總數(shù)，以10^6為單位；TranscriptLength(kb)：轉(zhuǎn)錄本長度10^3個堿基為單位。所有轉(zhuǎn)錄本的注：第一列表示轉(zhuǎn)錄本ID,其余各列表示轉(zhuǎn)錄本在不同樣品中的表達量樣品轉(zhuǎn)錄本表達量總體分各樣品FPKM密度分布對比注：圖中不同顏色的曲線代表不同的樣品，曲線上點的橫坐標表示對應(yīng)樣品FPKM的對數(shù)值，點的縱坐標表示概率密為了進一步查看單個樣品轉(zhuǎn)錄本表達水平分布的離散程度，以及直觀的比較不同樣品的整體轉(zhuǎn)錄本表達水平。采用箱線圖展示FPKM分各樣品FPKM箱線注：圖中橫坐標代表不同的樣品；縱坐標表示樣品表達量FPKM的對數(shù)值。該圖從表達量的總體離散角度來衡量各樣品表達水平。對于設(shè)立生物學重復的項目，評估生物學重復的相關(guān)性對于分析轉(zhuǎn)錄組數(shù)據(jù)非常重要。生物學重復的相關(guān)性不僅可以檢驗生物學實驗操作的可重復性；還可以評估差異表達的可靠性和輔助異常樣品的篩查。將皮爾遜相關(guān)系數(shù)r（PearsonCorrelationCoefficient）作為生物學重復相關(guān)性的評估指標。r2越接近1，說明兩個重復樣品相關(guān)性越強。百邁客保證對同一條件的所有生物學重復樣品進行同人同批樣品提取、建庫，同Run同Lane。對異常樣品進行詳細分析，并根據(jù)分析結(jié)果與溝通共識決定重新進行實驗，還是剔除異常樣品進行后續(xù)分析。為研究在不同實驗處理下的表達模式，首先對選取樣品的FPKM進行標準化歸一處理，然后對標準化的FPKM值進行K-聚類分析，同一類的在不同的實驗處理下具有相似的變化趨勢，具有相似變化趨勢的往往具有相似的功能轉(zhuǎn)錄本表達聚類折線示意注：x軸表示實驗分組，y軸表示標準化的FPKM。黑色線條表示這個cluster中的所有在不同實驗條件下相表達量的平均值的折線圖差異表達分表達具有時間和空間特異性，在兩個不同條件下，表達水平存在顯著差異的，差異表達（EG）。差異表達分析得到的基因集合叫做差異表達集，使用“_vs_”的方式命名。根據(jù)兩（組）樣品之間表達水平的相對高低，差異表達可以劃分為上調(diào)（Up-regulatedTranscript）和下調(diào)（Down-regulatedTranscript）。上調(diào)在樣品（組）B中的表達水平高于樣（組）A中的表達水平；反之為下調(diào)。上調(diào)和下調(diào)是相對的，由所給A和B的順序決定對于有生物學重復的樣本，Eeq適用于進行樣品組間的差異表達分析，獲得兩個生物學條件之間的差異表達集；對于沒有生物學重復的樣本，使用EBSeq進行差異分析。在差異表達過程中，將FoldChange≥2且FDR<0.01作為篩選標準。差異倍數(shù)（FoldChange）表示兩樣品（組）間表達量的比值。錯誤發(fā)現(xiàn)率（alseiscoeryate，R）是通過對差異顯著性p值（p-value）進行校正得到的。由于轉(zhuǎn)錄組的差異表達分析是對大量的表達值進行獨立的統(tǒng)計假設(shè)檢驗，會存在假陽性問題，因此在進行差異表達分析過程中，采用了公認的Benjamini-Hochberg校正方法對原有假設(shè)檢驗得到的顯著性p值（p-value）進行校正，并最終采用FDR作為差異表達篩選的關(guān)鍵差異表達集數(shù)目統(tǒng)計如下表差異表達數(shù)目統(tǒng)計部分結(jié)果展DEGDEGup-down-通過火山圖（VolcanoPlot）可以快速地查看在兩個（組）樣品中表達水平的差異，以及差異的統(tǒng)計學顯著性。差異表達火山圖下圖差異表達火山示意注：差異表達火山圖中的每一個點表示一個，橫坐標表示某一個在兩樣品中表達量差異倍數(shù)的對數(shù)值；縱坐標表示表達量變化的統(tǒng)計學顯著性的負對數(shù)值。橫坐標絕對值越大，說明表達量在兩樣品間的表達量倍數(shù)差異越大；縱坐標值越大，表明差異表達越顯著，篩選得到的差異表達越可靠。圖中綠色的點代表下調(diào)差異表達，紅色的點代表上調(diào)差通過MA圖可以直觀地查看的兩個（組）樣品的表達水平和差異倍數(shù)的整體分布。差異表達MA圖見下圖差異表達MA示意注：差異表達MA圖中每一個點代表一個。橫坐標為A值：log2(FPKM)，即兩樣品中表達量均值的對數(shù)值；縱坐標為M值：log2(FC)，即兩樣品間表達量差異倍數(shù)的對數(shù)值，用于衡量表達量差異的大小。圖中綠色的點代表下調(diào)差異表達，紅色的點代表上調(diào)差異表達，黑色的點代表非差異表達。對篩選出的差異表達做層次聚類分析，將具有相同或相似表達模式的進行聚類，部分差異表達聚類結(jié)果如下圖差異差異表 聚類部分結(jié)果展注：橫坐標代表樣品名稱及樣品的聚類結(jié)果，縱坐標代表的差異及的聚類結(jié)果。圖中不同的列代表不同的樣品，不同的行代表不同的。顏色代表了在樣品中的表達量水平log2（FPKM+1）。對差異表達進行數(shù)據(jù)庫的功能注釋，各差異表達集注釋到的數(shù)量統(tǒng)計見下表注釋的差異表達數(shù)量部分結(jié)果展DEGSwiss-差異表達GO注釋及富GO數(shù)據(jù)庫是GO組織（GeneOntologyConsortium）于2000年構(gòu)建的一個結(jié)構(gòu)化的標準生物學注釋系統(tǒng)，旨在建立及其產(chǎn)物知識的標準詞匯體系，適用于各個物種。GO注釋系統(tǒng)是一個有向無環(huán)圖，包含三個主要分支，即：生物學過程（BiologicalProcess），分子功能（MolecularFunction）和細胞組分（CellularComponent）。差異表達GO分類統(tǒng)計結(jié)果見下圖差異表達GO注釋分類統(tǒng)計示意注：橫坐標為GO分類，縱坐標左邊為數(shù)目所占百分比，右邊為數(shù)目。此圖展示的是在差異表達背景和全部背景下GO各二級功能的富集情況，體現(xiàn)兩個背景下各二級功能的地位，具有明顯比例差異的二級功能說明差異差異表達topGO有向無環(huán)示意注：對每個GO節(jié)點進行富集，最顯著的10個節(jié)點在圖中用方框表示，圖中還包含其各層對應(yīng)關(guān)系。每個方（或橢圓）內(nèi)給出了該GO節(jié)點的內(nèi)容描述和富集顯著性值。不同顏色代表不同的富集顯著性，顏色越深，顯著性越差異表達KEGG注釋及富在生物體內(nèi)，不同的產(chǎn)物相互協(xié)調(diào)來行使生物學功能，對差異表達的通路（Pathway）注釋分析有助于進一步解讀的功能。KEGG（KyotoEncyclopediaofGenesandGenomes）是系統(tǒng)分析功能、組信息數(shù)據(jù)庫，它有助于研究者把及表達信息作為一個整體網(wǎng)絡(luò)進行研究。作為是有關(guān)Pathway的主要公共數(shù)據(jù)庫(Kanehisa,2008），KEGG提供的整合代謝途徑(pathay)查詢，包括碳水化合物、核苷、氨基酸等的代謝及有機物的生物降解，不僅提供了所有可能的代謝途徑，而且對催化各步反應(yīng)的酶進行了全面的注解，包含有氨基酸序列、PDB庫的等等，是進行生物體內(nèi)代謝分析、代謝網(wǎng)絡(luò)研究的強有力工具。差異表 的通路注釋結(jié)果見下圖差異表達KEGG通路注釋示注：相對于對照組來說，紅色框標記的酶與上調(diào)因有關(guān)，綠色框標記的酶與下調(diào)有關(guān)。藍色框標記的酶與上調(diào)和下調(diào)均有關(guān)，框內(nèi)的數(shù)字代表酶的編號（ECnumber），而整個通路由多種酶催化的復雜生化反應(yīng)構(gòu)成，此通路圖中與差異表達相關(guān)的酶均用不同的顏色標出，根據(jù)研究對象間的差異，重點研究某些代謝通路相關(guān)的對差異表達KEGG的注釋結(jié)果按照KEGG中通路類型進行分類，分類圖如下圖所示差異表達KEGG分類示意注：縱坐標為KEGG代謝通路的名稱，橫坐標為注釋到該通路下的個數(shù)及其個數(shù)占被注釋上的總數(shù)的比例分析差異表達在某一通是否發(fā)生顯著差異（over-presentation）即為差異表達的通路富集分析。Pathway顯著性富集分析以KEGG數(shù)據(jù)庫中Pathway為單位，應(yīng)用超幾何檢驗，找出與整個背景相比，在差異表達中顯著性富集的Pathway。差異表達KEGG通路富集分析結(jié)果見下圖，圖中呈現(xiàn)了顯著性Q值最小的前20個通路。差異表達KEGG通路富集散點示意注：圖中每一個圓表示一個KEGG通路，縱坐標表示通路名稱，橫坐標為富集因子（EnrientFactor），表示差異中注釋到某通路的比例與所有中注釋到該通路的比例的比值。富集因子越大，表示差異表達在該通路中的富集水平越顯著。圓圈的顏色代表qvalue，qvalue為多重假設(shè)檢驗校正之后的Pvalue,qvalue越小，表示差異表達在該通路中的富集顯著性越可靠；圓圈的大小表示通路中富集的數(shù)目，圓圈越大，表示越多。差異表達COG分COG（ClusterofOrthologousGroupsofproteins）數(shù)據(jù)庫是基于細菌、藻類、真核生物的系統(tǒng)進化關(guān)系構(gòu)建得到的，利用COG數(shù)據(jù)庫可以對產(chǎn)物進行直系同源分類。差異表達COG分類統(tǒng)計結(jié)果見下圖差異表達COG注釋分類統(tǒng)計示意注：橫坐標為COG各分類內(nèi)容，縱坐標為數(shù)目。在不同的功能類中，所占多少反映對應(yīng)時期和環(huán)境下代謝或者生理偏向等內(nèi)容，差異表達eggNOG分eggNOG（evolutionarygenealogyofgenes:non-supervisedorthologousgroups）數(shù)據(jù)庫對直系同源類群進行了功能描述和功能分類的注釋,包含1133個物種的直系同源類群。eggNOG分類統(tǒng)計結(jié)果見下圖:差異表達eggNOG注釋分類統(tǒng)計示意圖注：橫坐標為eggNOG各分類內(nèi)容，縱坐標為數(shù)目。差異表達蛋白互作網(wǎng)TIG是收錄多個物種預(yù)測的和實驗驗證的蛋白質(zhì)-蛋白質(zhì)互作的數(shù)據(jù)庫，包括直接的物理互作和間接的功能相關(guān)。結(jié)合差異表達分析結(jié)果和數(shù)據(jù)庫收錄的互作關(guān)系對，構(gòu)建差異表達互作網(wǎng)絡(luò)。對于數(shù)據(jù)庫中包含的物種，可直接從數(shù)據(jù)庫中提取出目標集的互作關(guān)系對構(gòu)建互作網(wǎng)絡(luò)；對于數(shù)據(jù)庫中未收錄信息的物種，使用LAT軟件，將目的與數(shù)據(jù)庫中的蛋白質(zhì)進行序列比對，尋找同源蛋白，根據(jù)同源蛋白的互作關(guān)系對構(gòu)建互作網(wǎng)絡(luò)。構(gòu)建完成的蛋白質(zhì)互作網(wǎng)絡(luò)可導入Cytoscape軟件進行可視化。Cytoscape可視化的差異表達蛋白質(zhì)互作網(wǎng)絡(luò)如下圖差異表達蛋白質(zhì)互作網(wǎng)絡(luò)示意注：圖中的節(jié)點為蛋白質(zhì)，邊為互作關(guān)系?；プ骶W(wǎng)絡(luò)點(node)的大小與此節(jié)點的度(degree)成正比，即與此節(jié)點相連的邊越多，它的度越大，節(jié)點也就越大。節(jié)點的顏色與此節(jié)點的系數(shù)(clusteringcoefficient)相關(guān)，顏色梯度由綠到紅對應(yīng)系數(shù)的值由低到高；系數(shù)表示此節(jié)點的鄰接點之間的連通性好壞，系數(shù)值越高表示此節(jié)點的鄰接點之間的連通性越好。邊(edge)的寬度表示此邊連接的兩個節(jié)點間的互相作用的關(guān)系強弱，互相作用的關(guān)系越參考文AuKF,SebastianoV,AfsharPT,DurruthyJD,LeeL,WilliamsBA,etal.CharacterizationofthehumanESCtranscriptomebyhybridsequencing.ProcNatlAcadSciUSA.2013;110(50):E4821–30. :10.1073/pnas. SharonD,TilgnerH,GrubertF,SnyderM.Asingle-moleculelong-readsurveyofthehumantranscriptome.NatBiotechnol.2013;31(11):1009–14. :10.1038/nbt.2705.TreutleinB,GokceO,QuakeSR,SudhofTC.Cartographyofneurexinalternativesplicingmappedbysingle-moleculelong-readmRNAsequencing.ProcNatlAcadSciUSA.2014;111(13):E1291–9. :10.1073/pnas. GordonS, gE,SalamovA,ZhangJ,MengX,etal.(2014)Widespreadpolycistronictranscriptsinmushroom-formingfungirevealedbysingle-moleculelong-readmRNAsequencing.BioRxiv, Hackl,T.,Hedrich,R.,Schultz,J.,Forster,F.(2014).proovread:large-scalehigh-accuracyPacBiocorrectionthroughiterativeshortreadconsensus.Bioinformatics,30(21),3004-3011.WuTD,WatanabeCK(2005)GMAP:agenomicmapandalignmentprogramformRNAandESTsequences.Bioinformatics21:1859–1875.PMID: FelipeA.Sim?o,RobertM.Waterhouse,PanagiotisIoannidis,EvgeniaV.Kriventseva,andEvgenyM.Zdobnov.BUSCO:assessinggenomeassemblyandannotationcompletenesswithsingle-copyorthologs.Bioinformatics31.19(2015):3210-FoissacS,SammethM.ASTALAVISTA:dynamicandflexibleysisofalternativesplicingeventsincustomgenedatasets.NucleicAcidsResearch2007,35(WebServerissue):W297-9.AbdelghanySE,HamiltonM,JacobiJL,etal.Asurveyofthesorghumtranscriptomeusingsingle-moleculelongreads[J].NatureCommunications,2016,7:11706.HaasBJ,PapanicolaouA,YassourM,GrabherrM,BloodPD,BowdenJ,etal.DenovotranscriptsequencereconstructionfromRNA-sequsingtheTrinityplatformforreferencegenerationandysis.NatProtoc.2013;8:1494–1512.:L.Kong,Y.Zhang,Z.Q.Ye,X.Q.Liu,S.Q.Zhao,L.Wei,andG.Gao.2007.CPC:assesstheprotein-codingpotentialoftranscriptsusingsequencefeaturesandsupportvectormachine.NucleicAcidsRes36:W345-349.LiangSun,HaitaoLuo,DechaoBu,GuoguangZhao,KuntaoYu,ChanghaiZhang,YuanningLiu,RunShengChenandYiZhao*Utilizingsequenceintrinsiccompositiontoclassifyprotein-codingandlongnon-codingtranscripts.NucleicAcidsResearch(2013),:10.1093/nar/gkt646.WangL,ParkHJ,DasariS,WangS,KocherJP,LiW.CPAT:Coding-PotentialAssessmentToolusinganalignment-freelogisticregressionmodel.NucleicAcidsRes.2013Apr1;41(6):e74.:10.1093/nar/gkt006.LiJ,MaW,ZengP,etal.LncTar:atoolforpredictingtheRNAtargetsoflongnoncodingRNAs[J].BriefingsinBioinformatics,2015,16(5):806.ZhengY,JiaoC,SunH,RosliHG,PomboMA,ZhangP,BanfM,DaiX,MartinGB,GiovannoniJJ,ZhaoPX,RheeSY,FeiZ(2016)iTAK:aprogramforgenome-widepredictionandclassificationofplanttranscriptionfactors,transcriptionalregulators,andproteinkinases.MolecularPlant9:1667-1670.Hong-MeiZhang,TengLiu,Chun-JieLiu,ShuangyangSong,XiantongZhang,WeiLiu,HaiboJia,YuXue,andAn-YuanGuo.AnimalTFDB2.0:aresourceforexpression,predictionandfunctionalstudyofanimaltranscriptionfactors.Nucl.AcidsRes.(28January2015)43(D1):D76-D81.AltschulSF,MaddenTL,Sch?fferAA,ZhangJ,etal.GappedBLASTandPSIBLAST:ANewGenerationofProteinDatabaseSearchPrograms.NucleicAcidsResearch1997,25(17):3389-3402.DengYY,LiJQ,WuSF,ZhuYP,etal.IntegratedNRDatabaseinProteinAnnotationSystemandItsLocalization.ComputerEngineering2006.,32(5):71-74.ApR,BairochA,WuCH,BarkerWC,etal.UniProt:theUniversalProteinknowledgebase.NucleicAcidsResearch2004Jan1;32(Databaseissue):D115-9.AshburnerM,BallCA,BlakeJA,BotsteinD,etal.Geneontology:toolfortheunificationofbiology.NatureGenetics2000,25(1):25-29.TatusovRL,GalperinMY,NataleDA.TheCOGdatabase:atoolforgenomescaleysisofproteinfunctionsandevolution.NucleicAcidsResearch2000,28(1):33-36.KooninEV,FedorovaND,JacksonJD,etal.Acomprehensiveevolutionaryclassificationofproteinsencodedcompleteeukaryoticgenomes.Genomebiology,2004,5(2):FinnRD,BatemanA,ClementsJ,etal.Pfam:theproteinfamiliesdatabase.NucleicAcidsResearch,2013:KanehisaM,GotoS,KawashimaS,OkunoY,etal.TheKEGGresourcefordecipheringthegenome.NucleicAcidsResearch2004,32(Databaseissue):D277-D280.MckennaA,HannaM,BanksE,etal.TheGenomeysisToolkit:aMapReduceframeworkforyzingnext-generationDNAsequencingdata[J].GenomeResearch,2010,20(9):1297-303.CingolaniP,PlattsA,LeLW,etal.Aprogramforannotatingandpredictingtheeffectsofsinglenucleotidepolymorphisms,SnpEff[J].Fly,2012,6(2):80.TrapnellC,WilliamsBA,PerteaG,etal.TranscriptassemblyandficationbyRNA-Seqrevealsunannotatedtranscriptsandisoformswitchingduringcelldifferentiation.[J].NatureBiotechnology,2010,28(5):511-5.AndersS,HuberW.Differentialexpressionysisforsequencecountdata[J].GenomeBiology,2010,LengN,DawsonJA,ThomsonJA,etal.EBSeq:anempiricalBayeshierarchicalmodelforinferenceinRNA-seqexperiments[J].Bioinformatics,2013,29(8):1035.FranceschiniA,SzklarczykD,FrankildS,et