諾禾致源原核轉(zhuǎn)錄組生物信息分析結(jié)題報(bào)告2013年8月

一、建庫流TotalRNA上機(jī)差 差 文 列一、建庫流從RNA樣品到最終數(shù)據(jù)獲得，樣品檢測、建庫、每一個(gè)環(huán)節(jié)都會對數(shù)據(jù)質(zhì)量和數(shù)量產(chǎn)生影響，而數(shù)據(jù)質(zhì)量又會直接影響后續(xù)信息分析的結(jié)果。因此，獲得高質(zhì)量數(shù)據(jù)是保證生物信息分析正確、全面、可信的前提。為了從上保證數(shù)據(jù)的準(zhǔn)確性、可靠性，諾禾致源對樣品檢測、建庫、每一個(gè)生產(chǎn)步驟都嚴(yán)格把控，從根本上確保了高質(zhì)量數(shù)據(jù)的產(chǎn)出。流程圖如下：TotalRNA諾禾致源對RNA樣品的檢測主要包括4種方法瓊脂糖凝膠電泳分析RNA降解程度以及是否有污Nanodrop檢測RNA的純度（OD260/280比值Qubit對RNA濃度進(jìn)行精確定Agilent2100精確檢測RNA的完fragmentationbuffer將mRNA打斷成短片段，以mRNA為模板，用六堿基隨機(jī)引物（randomhexamers）合成一鏈cDNA，然后加入緩沖液、dNTPs（dTTP換為dUTP）和DNApolymeraseI合成二鏈cDNA，隨后利用AMPureXPbeads純化雙鏈cDNA。純化的雙鏈cDNA再進(jìn)行末端修復(fù)、加A尾并連接接頭，然后用AMPureXPbeads進(jìn)行片段大小選擇，最后消化二鏈后進(jìn)行PCR富集得到最終的cDNA文庫。構(gòu)建文庫構(gòu)建完成后，先使用Qubit2.0進(jìn)行初步定量，稀釋文庫至1ng/ul，隨后使用Agilent2100對文庫的insertsize進(jìn)行檢測，insertsize符合預(yù)期后，使用Q-PCR方法對文庫的有效濃度進(jìn)行準(zhǔn)確定量（文庫有效濃度＞2nM），以保證文庫質(zhì)量。上機(jī)庫檢合格后，把不同文庫按照有效濃度及目標(biāo)下機(jī)數(shù)據(jù)量的需求pooling后進(jìn)行 二、生物信息分析流程獲得原始序列(SequencedReads)后，在有相關(guān)物種參考序列或參 組的情況下，通過 程進(jìn)行生物信息分析三、項(xiàng)目結(jié)果說明1高通量(如illuminaHieqM200MiSq等平臺)得到的原始圖像數(shù)據(jù)文件經(jīng)堿基識別(Baseallng)分析轉(zhuǎn)化為原始序列(Sequenedeas)，我們稱之為RawDta或RawRes，結(jié)果以FASTQ(簡稱為fq)文件格式，其中包含序列(d)的序列信息以及其對應(yīng)的質(zhì)量信息。FASTQ格式文件中每個(gè)read由四行描述，如下@EAS139:136:FC706VJ:2:2104:15343:1973931:Y:18:ATCACG其中第一行以“@”開頭，隨后為illumina標(biāo)識符(SequeneIdetiies)和描述文字(選擇性部分)；第二行是堿基序列；第三行以“+”開頭，隨后為illumia標(biāo)識符(選擇性部分)；第四行是對應(yīng)序列的質(zhì)量(Cockta.)。illumina標(biāo)識符詳細(xì)信息如下UniqueinstrumentRunFlowcell2FlowcellTilenumberwithintheflowcell'x'-coordinateoftheclusterwithinthe'y'-coordinateoftheclusterwithinthe1Memberofapair,1or2(paired-endormate-pairreadsYYifthereadfailsfilter(readisbad),N0whennoneofthecontrolbitsareon,otherwiseitisanevenIndex公式一：Qphred=-10log10(e)illuminaCasava1.8版本 錯(cuò)誤率質(zhì)量值對應(yīng)字符.5?I2數(shù)據(jù)質(zhì)量評 錯(cuò)誤率分布檢每個(gè)堿基錯(cuò)誤率是通過Phred數(shù)值(Phredscore,Qphred)通過公式1轉(zhuǎn)化得到，而Phred數(shù)值是在堿基識別(BaseilluminaCasava1.8版本堿基識別與Phred分值之間的簡明對應(yīng)關(guān)Phred分不正確的堿基識別堿基正確識別率Q-9999錯(cuò)誤率與堿基質(zhì)量有關(guān)，受儀本身、試劑、樣品等多個(gè)因素共同影響。對于RNA-seq技術(shù)，錯(cuò)誤率分布具有兩個(gè)特點(diǎn)：(1)錯(cuò)誤率會隨著序列(SequencedReds)長度的增加而升高，這是由于過程中化學(xué)試劑的消耗而導(dǎo)致的，并且為illumina高通量平臺都具有的特征。(2)前6個(gè)堿基的位置也會發(fā)生較高的錯(cuò)誤率，而這個(gè)長度也正好等于在RNA-seq庫程中反轉(zhuǎn)錄所需要的隨機(jī)引物的長度。所以推測前6個(gè)堿基錯(cuò)誤率較高的原因?yàn)殡S機(jī)引物和RNA模版的不完全結(jié)合(Jiangetal.。圖2.1錯(cuò)誤率分布橫坐標(biāo)為reads的堿基位置，縱坐標(biāo)為單堿GCGC含量分布檢查用于檢測有無AT、GC分離現(xiàn)象，而這種現(xiàn)象可能是或者建庫所帶來的，并且會影響后續(xù)的定量分析在illumna平臺的轉(zhuǎn)錄組中，反轉(zhuǎn)錄成cDNA時(shí)所用的6bp的隨機(jī)引物會引起前幾個(gè)位置的核苷酸組成存在一定的偏好性。而這種偏與的物種和環(huán)境無關(guān)，但會影響轉(zhuǎn)錄組的均一化程度(Hansenetl.)。除此之外，理論上G和C堿基及A和T堿基含量每個(gè)循環(huán)應(yīng)分別相等，且整個(gè)過程穩(wěn)不變，呈水平線。對于DE來說，由于隨機(jī)引物擴(kuò)增偏差等原因，常常會導(dǎo)致在得到的每個(gè)ed67個(gè)堿基有較大的波，這種波動(dòng)屬于正常情況。對于鏈特異性文庫，由于只保留鏈，G和C堿基AT堿基含量不等。對于鏈特異性文庫，由建庫過程中只保留一鏈，G和C基及AT堿基含量不等。圖2.2GC含量分布橫坐標(biāo)為reads的堿基位置，縱坐標(biāo)為單堿基所占的比例；不同顏色代表不 數(shù)據(jù)過得到的原始序列，里面含有帶接頭的、低質(zhì)量的reads，為了保證信息分析質(zhì)量，必須對rawreads進(jìn)行過濾，得到cleanreads，后續(xù)分析都基于cleanreads。數(shù)據(jù)處理的步驟如下去除帶接頭(adapter)的去除N(N表示無法確定堿基信息)的比例大于10%的去除低質(zhì)量readsRNA-seq的接頭(Adapter,OligonucleotidesequencesforTruSeqTMRNAandDNASamplePrepKits)信息：RNA5’Adapter(RA5),part# 5’-AATGGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT-RNA3’Adapter(RA3),part 5’-GATCGGAAGAGCACACGTCTGAACTCCAGTCAC(6位index)ATCTCGTATGCCGTCTTCTGCTTG-圖2.3原始數(shù)據(jù)過濾結(jié)數(shù)據(jù)質(zhì)量情況匯

表2.4數(shù)據(jù)產(chǎn)出質(zhì)量情況一覽SampleRawCleanErrorGCGG數(shù)據(jù)質(zhì)量情Rawreads：統(tǒng)計(jì)原始序列數(shù)據(jù)，以四行為一個(gè)單位，統(tǒng)計(jì)每個(gè)文件的序列的個(gè)數(shù)Cleanreads：計(jì)算方法同RawReads，只是統(tǒng)計(jì)的文件為過濾后的數(shù)據(jù)。后續(xù)的生物信息分析都是基于CleanreadsCleanbases：序列的個(gè)數(shù)乘以序列的長度，并轉(zhuǎn)化為以G為單位Errorrate：通過公式1Q20、Q30：分別計(jì)算Phred數(shù)值大于20、30的堿基占總體堿基的百分GCcontent：計(jì)算堿基G和C的數(shù)量總和占總的堿基數(shù)針對細(xì)菌等密度較高的生物，我們用Bowtie2將過濾后的序列進(jìn)行組定位分析。如果參考組選擇合適，而且相關(guān)實(shí)驗(yàn)不存在污染，實(shí)驗(yàn)所產(chǎn)生的序列的定位的百分比正常情況下會高于70%(TotalMappedReadsorFragents)，其中具有多個(gè)定位的序列(MultipleMapedRadsrFraments)占總體的百分比通常不會超過10%。Reads與參考組比對情況統(tǒng)表3.1Reads與參考組比對情況一覽SampleTotalTotal(95MultiplniquelyRead-Read-(47ReadsmaptoReadsmapto'-Readsmappedinproper191477147063比對結(jié)果統(tǒng)Totalreads：序列經(jīng)過數(shù)據(jù)過濾后的數(shù)量統(tǒng)計(jì)(Cleandata)Totalmppd：能定位到組上的序列的數(shù)量的統(tǒng)計(jì)；一般情況下，如果不存在污染并且參考組選擇合適的情況下，這部分?jǐn)?shù)據(jù)的百分比大于70%。Multiplemapped：在參考序列上有多個(gè)比對位置的序列的數(shù)量統(tǒng)計(jì)；這部分?jǐn)?shù)據(jù)的百分比一般會小于10%Uniquelymapped：在參考序列上有唯一比對位置的序列的數(shù)量統(tǒng)計(jì)Readsmapto'+'，Readsmapto'- 序列比對 組上正鏈和負(fù)鏈的統(tǒng)計(jì)Reads在參考組不同區(qū)域的分布情根據(jù)組的注釋信息，對Totalmappedreads比對到組上的各個(gè)部分的情況進(jìn)行統(tǒng)計(jì)。正常情況下，區(qū)的reads定圖3.2Reads在參考組不同域的分布情Reads在上的密度分布情對Totalmapedreds的比對到組（分正負(fù)鏈）的密度進(jìn)行統(tǒng)計(jì)，下圖是隨機(jī)抽取部分reads，展示其在上的map圖3.3Reads在上的密度分布最外圈是組；中間的灰色背景區(qū)是實(shí)際的抽取的reads的分布情況，紅色map到正鏈，藍(lán)色到負(fù)鏈；最里面的圓圈，橘黃色為正鏈coverage分布，綠色為負(fù)鏈coverage分布，超過所有coverage集均值+2倍標(biāo)準(zhǔn)差的奇異點(diǎn)被舍棄Reads我們提供RNA-seqReads在組上比對結(jié)果的bam格式文件，部分物種還提供相應(yīng)的參考組和注釋文件，并推薦使用IGV(IntegraiveGnomicsViewr)瀏覽器對bam文件進(jìn)行可視化瀏覽。IGV瀏覽器具有以下特點(diǎn)：(1)能在不同尺度下顯示單個(gè)或多個(gè)讀段在組上的位置，包括讀段在各個(gè)上的分布情況和在注釋的外顯子、內(nèi)含子、剪接接合區(qū)、間區(qū)的分布情況等；(2)能在不同尺度下顯示不同區(qū)域的讀段豐度，以反映不同區(qū)域的轉(zhuǎn)錄水平；(3)能顯示及其剪接異構(gòu)體的注釋信息；(4)能顯示其他注釋信息；(5)既可以從服務(wù)器端各種注釋信息，又可以從本地加載注釋信息。IGV瀏覽器使用方法可參考我們提供的使用說明文檔(IGVQuickStart.pd)。圖3.4IGV瀏覽器4表達(dá)水平分一個(gè)表達(dá)水平的直接體現(xiàn)就是其轉(zhuǎn)錄本的豐度情況，轉(zhuǎn)錄本豐度程度越高，則表達(dá)水平越高。在RNA-seq分析中，我們可以通過定位到組區(qū)域或編碼區(qū)的序列(reads)的計(jì)數(shù)來估計(jì)的表達(dá)水平。Reads計(jì)數(shù)除了與的真實(shí)表達(dá)水平成正比外，還與的長度和深度成正相關(guān)。為了使不同、不同實(shí)驗(yàn)間估計(jì)的表達(dá)水平具有可比性，人們引入了RPKM的概念，RPKM(RadsPerKioassprillonreds)是每百萬reads中來自某一每千堿基長度的reads數(shù)目。RPKM同時(shí)考慮了測序深度和長度對reads計(jì)數(shù)的影響，是目前最為常用的表達(dá)水平估算方法(Mortaavietal,008)。結(jié)果文件分別統(tǒng)計(jì)了不同表達(dá)水平下的數(shù)量以及單個(gè)的表達(dá)水平。一般情況下，RPKM數(shù)值0.1或者1作為判斷是否表達(dá)的閾值，不同的文獻(xiàn)所采用的閾值不同。表4.1不同表達(dá)水平區(qū)間的數(shù)量統(tǒng)計(jì)RPKM0-1-3-15-2666(60表4.2表達(dá)水平統(tǒng)計(jì)76305RNA-seq表達(dá)水平的飽和曲線檢定量飽和曲線檢查反映了表達(dá)水平定量對數(shù)據(jù)量的要求。表達(dá)量越高的，就越容易被準(zhǔn)確定量；反之，表達(dá)量低的基因，需要較大的數(shù)據(jù)量才能被準(zhǔn)確定量。表達(dá)水平的飽和曲線的體算法描述如下：分別對10%、20%、30%……90%的總體數(shù)據(jù)單獨(dú)進(jìn)行定量分析，并把所有數(shù)據(jù)條件下得到的的表達(dá)水平作為最終的數(shù)值。用每個(gè)百分比條件下求出的單個(gè)的RPKM數(shù)值和最終對應(yīng)的表達(dá)水平數(shù)值進(jìn)行比較，如果差異小于15%，則認(rèn)為這個(gè)在這個(gè)條件下定量是準(zhǔn)確的。圖5.1定量飽和曲線檢查分布橫坐標(biāo)代表定位 組上的reads數(shù)占總reads數(shù)的百分比，縱坐標(biāo)代表定量誤差在15%以內(nèi)的的比RNA-Seq生物學(xué)重復(fù)是任何生物學(xué)實(shí)驗(yàn)所必須的，高通量技術(shù)也不例外(Hansenetl.)。生物學(xué)重復(fù)主要有兩個(gè)用途：一個(gè)是證明所涉及的生物學(xué)實(shí)驗(yàn)操作是可以重復(fù)的且變異不大，另一個(gè)為后續(xù)的差異分析所需要的。樣品間表達(dá)水平相關(guān)性是檢驗(yàn)實(shí)驗(yàn)可靠性和樣本選擇是否合理性的重要指標(biāo)。相關(guān)系數(shù)越接近1，表明樣品之間表達(dá)模式的相似度越高。Encode計(jì)劃建議生物學(xué)重復(fù)皮爾遜相關(guān)系數(shù)的平方(R2)大于0.92(理想的取樣和實(shí)驗(yàn)條件下)。具體的項(xiàng)目操作中，我們要求R2至少要大于0.8，否則需要對樣品做出合適的解釋，或者重新進(jìn)行實(shí)驗(yàn)。圖5.2RNA-Seq重復(fù)相關(guān)性檢R2：pearson6差異表達(dá)分 表達(dá)水平對通過所有的RPKM的分布圖以及盒形圖對不同實(shí)驗(yàn)條件下的表達(dá)水平進(jìn)行比較。對于同一實(shí)驗(yàn)條件下的重復(fù)樣品，最終RPKM為所有重復(fù)數(shù)據(jù)的平均值圖6.1不同實(shí)驗(yàn)條件下表達(dá)水平比對RPKM分布圖(圖一)的橫坐標(biāo)為log10(RPKM),縱坐標(biāo)為的密度。RPKM盒形圖(圖二)的橫坐標(biāo)為樣品名稱，縱坐標(biāo)為log10(RPKM)，每個(gè)區(qū)域的盒形圖對差異表達(dá)列etal,2010）進(jìn)行分析：該分析方法基于的模型是負(fù)二項(xiàng)分布，第i個(gè)在第j個(gè)樣本中的readcount值為Kij，則Kij～對于無生物學(xué)重復(fù)的樣品，先采用TMM對readcount數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，之后用DEGseq進(jìn)行差異分析。差異表達(dá)列表如下表6.2差異列3105.1372e-3.9323e-301.6497e-2.3834e-20344e-6.0811e-48119478e-1.2765e-差異列表主要包括的內(nèi)容Gene_id:編readcount_Sample1：校正后樣品1的readcountreadcount_Sample2：校正后樣品2的readcountlog2FoldChange:pvalue(pval):統(tǒng)計(jì)學(xué)差異顯著性qvalue(padj):校正后的pvalue。qvalue越小，表 表達(dá)差異越顯差異表達(dá)篩用火山圖可以推斷差異的整體分布情況，對于無生物學(xué)重復(fù)的實(shí)驗(yàn)，為消除生物學(xué)變異，我們從差異倍數(shù)和顯著水平兩個(gè)水平進(jìn)行評估，對差異進(jìn)行篩選，閾值設(shè)定一般為:|log2(FodChange)|>1且qvalue<0.005。對于有生物學(xué)重復(fù)的實(shí)驗(yàn)，由于DESeq已經(jīng)進(jìn)行了生物學(xué)變異的消除，我們對差異篩選的標(biāo)準(zhǔn)一般為:padj<005。圖6.3差異火山有顯著性差異表達(dá) 用紅色點(diǎn)表示；橫坐標(biāo)代 在不同樣本中表達(dá)倍數(shù)變化；縱坐標(biāo)代 表達(dá)量變化差異的統(tǒng)計(jì)學(xué)顯著差異聚類分聚類分析用于判斷差異在不同實(shí)驗(yàn)條件下的表達(dá)模式；通過將表達(dá)模式相同或相近的成類，從而識別未知的功能或已知的未知功能；因?yàn)檫@些同類的可能具有相似的功能，或是共同參與同一代謝過程或細(xì)胞通路。以不同實(shí)驗(yàn)條件下的差異的RPKM值為表達(dá)水平，做層次聚類(hierarchicallusering)分析，不同的顏色的區(qū)域代表不同的聚類分組信息，同組內(nèi)的表達(dá)模式相近，可能具有相似的功能或參與相同的生物學(xué)過程。除了差異表達(dá)量rpkm層次聚類分析，我們對還分別用H-cluster、K-means和SOM等三種方法對差異的相對表達(dá)水平值log2(ratios)進(jìn)行聚類。不同的聚類算法分別將差異分為若干cluster，同一cluster中的在不同的處理?xiàng)l件下具有相似的表達(dá)水平變化趨勢。圖6.4差異聚類圖一為整體rpkm層次聚類圖，以log10RPKM值進(jìn)行聚類，紅色表示高表達(dá)，藍(lán)色表示低表達(dá)。顏色從紅到藍(lán)，表示log10(RPKM)從大到??；圖二為log2(ratios)折線圖，每個(gè)子圖中的灰色線條表示一個(gè)cluster中的在不同實(shí)驗(yàn)條件下相對表達(dá)量的折線圖，藍(lán)色線條表示這個(gè)cluster中的所有在不同實(shí)驗(yàn)條件下相對表達(dá)量的平均值的折線圖，x軸表示實(shí)驗(yàn)條件，y軸表示相對表達(dá)量差異GO富集分GeneOnology（簡稱GO,）是功能國際標(biāo)準(zhǔn)分類體系。根據(jù)實(shí)驗(yàn)?zāi)康暮Y選差異后，研究差異在GeneOntology中的分布狀況將闡明實(shí)驗(yàn)中樣本差異在功能上的體現(xiàn)。GO富集分析方法為GOseq（Youngetal,2010）,此方法基于Walleniusnon-entralhype-geomtricistributon。相對于普通的超幾何分布(Hyper-geomericdistribuion)，此分布的特點(diǎn)是從某個(gè)類別中抽取的概率與從某個(gè)類別之外抽取一個(gè)的概率是不同的，這種概率的不同是通過對長度的偏進(jìn)行估計(jì)得到的，從而能更為準(zhǔn)確地計(jì)算出GOterm被差異富集的概率。差異GO富集列表7.1差異GO富集列cell1.2529e-1.7778e-localizationof1.2529e-1.7778e-ciliaryorflagellar4.2065e-39794e-cellularcomponent1.717e-12182e-結(jié)果表格詳細(xì)內(nèi)容如下GO_accession：GeneOntology數(shù)據(jù)庫中唯一的Description：GeneOntology功能的描述Term_type：該GO的類別( onent：細(xì)胞組分；biological_prcess：生物學(xué)過程；molecular_funtion：功能)Over_represented_pValue：富集分析統(tǒng)計(jì)學(xué)顯Corrected_pValue：矯正后的P-Value，一般情況下，P-value<0.05該功能為富集DEG_item：與該GO相關(guān)的差異的數(shù)DEG_list：GO注釋的差 數(shù)差異GO富集DAG有向無環(huán)圖(DirectdAycicGrah，DAG)為差異GO富集分析結(jié)果的圖形化展示方式，分支代表包含關(guān)系，從上至下所定義的功能范圍越來越小，一般選取GO富集分析的結(jié)果前10位作為有向無環(huán)圖的主節(jié)點(diǎn)，并通過包含關(guān)系，將相關(guān)聯(lián)的GOTerm一起展示，顏色的深淺代表富集程度。我們的項(xiàng)目中分別繪制生物過程(bioogiclroess)、分子功能(moleclarfucton)和細(xì)胞組分(celulrompnet)的DAG圖。圖7.2GO富集有向無每個(gè)節(jié)點(diǎn)代表一個(gè)GO術(shù)語，方框代表的是富集程度為TOP10的GO，顏色的深淺代表富集程度，顏色越深就表示富集程度越高，每個(gè)節(jié)點(diǎn)上展示了該TERM及富集分析的p-value差異GO富集柱狀差異GO富集柱狀圖，直觀的反映出在生物過程(biologicalprocess)、細(xì)胞組分(cellularcomponent)和分子功能(molecularfunction)富集的GOterm上差異的個(gè)數(shù)分布情況。我們挑選了富集最顯著的30個(gè)GOterm在圖中展示，如果不足30圖7.3GO富集柱狀圖一：縱坐標(biāo)為富集的GOterm，橫坐標(biāo)為該term中差異個(gè)數(shù)。不同顏色用來區(qū)分生物過程細(xì)胞組分和分子功能，帶“*”為富集的GOterm圖二：對圖一中的GO，按生物過程細(xì)胞組分和分子功能三大類別及差異 差異KEGG富集分在生物體內(nèi)，不同相協(xié)調(diào)行使其生物學(xué)功能，通過Pathway顯著性富集能確定差異表達(dá)參與的最主要生化代謝途徑和信號轉(zhuǎn)導(dǎo)途徑。KEGG（KyotoEncyclopediafGenesndGenos）是有關(guān)Pathway的主要公共數(shù)據(jù)庫(Kanehisa,2008）。Pathwa性富集分析以KEGGPathway為單位，應(yīng)用超幾何檢驗(yàn)，找出與整個(gè)組背景相比，在差異表達(dá)中顯著性富集的Pathway。差異KEGG富集列表8.1差異KEGG富集列SampleP-CorrectedP-FlagellarKEGGponentKEGG553887613552e-252018864166e-HistidineKEGG221572423431e-6.72103017741e-BacterialKEGG0結(jié)果表格詳細(xì)內(nèi)容如下#Term：KEGG通路的描述信息Id：KEGG數(shù)據(jù)庫中通路唯一Samplenumber：該通路下差異的個(gè)數(shù)Backgroundnumber：該通路下的個(gè)數(shù)P-value：富集分析統(tǒng)計(jì)學(xué)顯著CorrectedP-value：矯正后的統(tǒng)計(jì)學(xué)顯著水平，一般情況下，P-value0.05該功能為富集差異KEGG富集散點(diǎn)散點(diǎn)圖是KEGG富集分析結(jié)果的圖形化展示方式。在此圖中，KEGG富集程度通過Richfactr、Qvlue和富集到此通的個(gè)數(shù)來衡量。其中Richfactr指差異表達(dá)的中位于該pathway條目的數(shù)目與所有有注釋中位于該pathway條目的總數(shù)的比值。Richfacor越大，表示富集的程度越大。Qvalue是做過多重假設(shè)檢驗(yàn)校正之后的Pvalue，Qalue的取值范圍為[0,1]，越接近于零，表示富集越顯著。我們挑選了富集最顯著的0條typwy0圖8.2差異KEGG富集散點(diǎn)縱軸表示pathway名稱，橫軸表示 factor，點(diǎn)的大小表示此pathway中差異表達(dá)個(gè)數(shù)多少，而點(diǎn)的顏色對應(yīng)于不同的Qvalue范富集KEGG通路將差異富集出的通路圖展示出來，通路圖中，包含上調(diào)的KO節(jié)點(diǎn)標(biāo)紅色，包含下調(diào)的KO節(jié)點(diǎn)標(biāo)綠色，包含上下調(diào)的標(biāo)黃色。鼠標(biāo)懸停于標(biāo)記的KO節(jié)點(diǎn)，彈出差異細(xì)節(jié)框，標(biāo)色同上，括號中數(shù)字為log2(Foldchang)。以上步驟可脫機(jī)實(shí)現(xiàn)，如連接互聯(lián)網(wǎng)，點(diǎn)擊各個(gè)節(jié)點(diǎn)，可以連接到KEGG數(shù)據(jù)庫中各個(gè)KO的具體信息頁。圖8.3著富集的KEGGpathway代謝通路SNP和InDelSNP全稱SingleNucleoidePlymorhisms，是指在組上由單個(gè)核苷酸變異形成的遺傳標(biāo)記，其數(shù)量很多，多態(tài)性豐富。從理論上來看每一個(gè)SNP位點(diǎn)都可以有4種不同的變異形式，但實(shí)際上發(fā)生的只有兩種，即轉(zhuǎn)換和顛換，二者之比為1:2。SNP在CG序列上出現(xiàn)最為頻繁，而且多是C轉(zhuǎn)換為T，原因是CG中的C常為甲基化的，自發(fā)地脫氨后即成為胸腺嘧啶。一般而言，SNP是指變異頻率大于1%的單核苷酸變異。InDel(insertio-deleion)是指相對于參考組，樣本中發(fā)生的小片段的插入缺失，該插入缺失可能含一個(gè)或多個(gè)堿基。我們通過samtools和picard-tools等工具對比對結(jié)果進(jìn)行坐標(biāo)排序、去掉重復(fù)的reads等處理，最后通過變異檢測軟件GATK(McKennaetal2010)分別進(jìn)行SNPCalling和InDelCalling，并對原始結(jié)果進(jìn)行過濾，得到如下表形式的分析結(jié)果。其中表9SNP分析結(jié)TCTGAGAG#CHROM：SNP位點(diǎn)所在othercoloums：每 該位點(diǎn) 型(0與REF一致；1與ALT一致；.缺少數(shù)據(jù)支持用Rockhopper軟件(R.McClur,etal,213)將結(jié)果根據(jù)參照參考組進(jìn)行組裝，并與已注釋的模型進(jìn)行比較，發(fā)現(xiàn)新的未知轉(zhuǎn)錄本區(qū)域。通過Blastx與nr庫作比對，對新預(yù)測的轉(zhuǎn)錄本區(qū)域進(jìn)行注釋，將能注釋上的轉(zhuǎn)錄本區(qū)域作為具有編碼潛能的新轉(zhuǎn)錄本區(qū)域。表10新轉(zhuǎn)錄本注釋結(jié)- -- +轉(zhuǎn)錄本編轉(zhuǎn)錄起始位轉(zhuǎn)錄終止位鏈方NR_GI:比對到NR庫的的GenbankNR_ID:比對到NR庫 的序列 結(jié)構(gòu)分原核生物功能上關(guān)的幾個(gè)往往串聯(lián)排列在一起，構(gòu)子(Operon)結(jié)構(gòu)作為的表達(dá)單元，受上游共同的調(diào)控區(qū)和下游轉(zhuǎn)錄終止信號的調(diào)控。錄時(shí)，幾個(gè)轉(zhuǎn)錄在一條mRNA鏈上，再分別翻譯成各自不同的蛋白質(zhì)。原核生物結(jié)構(gòu)和調(diào)控模式如下圖所示：我們通過Rockhopper軟件，根據(jù)reads在組上的分布情況，對子、轉(zhuǎn)錄起始位點(diǎn)(TranscriptionStartSite,TSS)和轉(zhuǎn)錄終止位點(diǎn)(TranscriptionTerminationSite,TTS)進(jìn)行預(yù)測。然后提取轉(zhuǎn)錄起始位點(diǎn)上游700bp序列，通過軟件Virtual

表11.1子預(yù)測結(jié)Numberof+4-3pdxJ,XAC0013,+2XAC0015,+2XAC0021,Start：第一個(gè)的起始坐Stop：最后一個(gè)的終止坐Strand：鏈方NumberofGenes：個(gè)Genes：名TSS和TTS預(yù)表11.2TSS和TTS預(yù)測結(jié)++++TSS：轉(zhuǎn)錄起始位點(diǎn)TTS：轉(zhuǎn)錄終止位點(diǎn)Strand：鏈方Genes：名表11.3啟動(dòng)子預(yù)Sequence(species)StartAbrB|Bacillussubtilis(strain+7AlgU(-10)|Pseudomonasaeruginosa(strainATCC15692/+5AlgU(-10)|Pseudomonasaeruginosa(strainATCC15692/+5AlgU(-10)|Pseudomonasaeruginosa(strainATCC15692/+5SequenceID入的序列(species)：位置矩陣StartPosition：啟動(dòng)子EndPosition動(dòng)子終止位點(diǎn)坐Strand方Score:精準(zhǔn)Sequence序UTRUTR我們根據(jù)轉(zhuǎn)錄起始位點(diǎn)(轉(zhuǎn)錄終止位點(diǎn))和翻譯起始位點(diǎn)(翻譯終止位點(diǎn))信息，提取5'UTR(3'UTR)序列，并對其長度分布情況進(jìn)行統(tǒng)計(jì)。針對5’UTR，用RBSfinder軟件對SD序列進(jìn)行預(yù)測；針對3’UTR，用TransTermHP軟件對不依賴σ因子的終止子進(jìn)行預(yù)測。圖12.1UTR長度分布橫軸表示UTR長度區(qū)間，縱軸是不同區(qū)間UTR密度的統(tǒng)計(jì)，紅色虛線代表UTR5’UTRSD序列預(yù)表12.2SD序列預(yù)測結(jié)++-+gene_id：編Start：起始坐Stop：終止坐Strand：鏈方Pattern：SD序列信Position：SD序列起始坐3’UTR不依賴σ因子的終止子預(yù)表12.3終止子預(yù)----GCGCGC-gene_id：編Term_start：終止子起始坐Term_end：終止子終止坐strand：鏈方5'_tail：莖環(huán)結(jié)構(gòu)5'5'_stem：莖環(huán)結(jié)構(gòu)5'loop：莖環(huán)結(jié)構(gòu)環(huán)序3'_stem：莖環(huán)結(jié)構(gòu)3'3'_tail：莖環(huán)結(jié)構(gòu)3'Cis-natualanisensetrancript(cis-ATs)反義轉(zhuǎn)錄本是由源DNA鏈相同區(qū)域轉(zhuǎn)錄的內(nèi)生RNA分子，與正義轉(zhuǎn)錄本存在部分收斂或分散方向的重復(fù)。據(jù)目前研究發(fā)現(xiàn)，反義轉(zhuǎn)錄本是重要的生物機(jī)理，主要通過表觀遺傳學(xué)上的改變，對進(jìn)行調(diào)控。反義轉(zhuǎn)錄本分為三種類型：enclosed（全部包含）、convergent(3'3')和divergent(5'-5')。對于鏈特異性建庫的RNA-seq數(shù)據(jù)，可以鑒定其反義轉(zhuǎn)錄本在組上的位置、種類以及數(shù)量等。表13反義轉(zhuǎn)錄本預(yù)測結(jié)(ISxac3(ISxac3(ISxac3(ISxac3plus_transcript_id：正鏈轉(zhuǎn)錄本編plus_start：正鏈轉(zhuǎn)錄本起始位plus_end：正鏈轉(zhuǎn)錄本終止plus_length：正鏈轉(zhuǎn)錄本長(11)types：反義轉(zhuǎn)錄本類型(包括：enclosed、convergent和divergent) sRNA細(xì)菌中，長度在50～500nt的非編碼RNA通常定義為小RNA(smallNA,sRN)。用Rockhopper軟件發(fā)現(xiàn)新的間區(qū)轉(zhuǎn)錄本，通過Blastx與nr庫作比對，對新預(yù)測的轉(zhuǎn)本區(qū)域進(jìn)行注釋，將注釋不上的轉(zhuǎn)錄本作為候選的非編碼sRNA。通過RNAfold軟件和IntaRNA對候選的sRNA分別進(jìn)行二級結(jié)構(gòu)預(yù)測和靶預(yù)測。sRNA圖14.1sRNA長度分布橫軸表示sRNA長度區(qū)間，縱軸是不同區(qū)間sRNA密度的統(tǒng)計(jì)，紅色虛線代表sRNA長度的均sRNA圖14.2sRNA二級結(jié)sRNA 預(yù)

表14.3sRNA靶預(yù)測結(jié)11--382--1--416--1--234--93--196--sRNA_id：sRNA編mRNA_id：靶編energy(kcal/mol)：自由sRNA_position：sRNA互補(bǔ)位mRNA_position： 互補(bǔ)位四、參考文獻(xiàn)Anders,S.(2010).HTSeq:ysinghigh-throughputsequencingdatawithAnders,S.,andHuber,W.(2010).Differentialexpressionysisforsequencecountdata.GenomeBiol.(DESeq)Anders,S.andHuber,W.(2012).DifferentialexpressionofRNA-Seqdataatthegenelevel-theDESeqpackage.(DEseq)Busch,A.,A.S.Richter,etal.(2008).IntaRNA:efficientpredictionofbacterialsRNAtargetsincorporatingtargetsiteaccessbilityandseedregions.Hofacker,I.L.andP.F.Stadler(2006).MemoryefficientfoldingalgorithmsforcircularRNAsecondarystructures.Bioinformatics.(RNAfold)Kanehisa,M.,M.Araki,etal.(2008).KEGGforlinkinggenomestolifeandtheenvironment.Nucleicacidsresearch.(KEGG)Kingsford,C.L.,K.Ayanbule,etal.(2007).Rapid,accurate,computationaldiscoveryofRho-independenttranscriptionterminatorsilluminatestheirrelationshiptoDNAuptake.Genomebiology.(TransTermHP)Langmead,B.,Trapnell,C.,Pop,M.&Salzberg,S.L.(2009).Ultrafastandmemory-efficientalignmentofshortDNAsequencestothehumangenome.GenomeBiol.(Bowtie)Langmead,B.andS.L.Salzberg(2012).Fastgapped-readalignmentwithBowtie2.Naturemethods.(BowtieMao,X.,Cai,T.,Olyarchuk,J.G.,Wei,L.(1995).AutomatedgenomeannotationandpathwayidentificationusingtheKEGGOrthology(KO)asaco