諾禾致源lncrna測(cè)序生物信息分析結(jié)題報(bào)告

ProvidingProvidingadvancedgenomicLongnoncodingRNA生物信息分析20139RNA諾禾致源生物信息科技地址：市海淀區(qū)學(xué)清路38號(hào)金碼B座21層 ：010-8283lncRNA一、建庫(kù)流TotalRNA樣品檢文庫(kù)構(gòu)庫(kù)上原始序列數(shù)數(shù)據(jù)質(zhì)量評(píng)估2.1錯(cuò)誤率分布檢2.2GC含量分布檢2.3數(shù)據(jù)過(guò)2.4數(shù)據(jù)產(chǎn)出情況匯參考序列比對(duì)分Reads與參 組比對(duì)情況統(tǒng)Reads在參 組不同區(qū)域的分布情Reads 上的密度分布情Reads比對(duì)結(jié)果IGV可視化瀏表達(dá)分已知注釋類型含量分已知表達(dá)水平分RNA-seq整體質(zhì)量樣品間相關(guān)性檢樣品間聚類及PCA分均一性分布檢轉(zhuǎn)錄本拼cufflinks拼scripture拼候選lncRNA篩基本篩編碼潛能篩重現(xiàn)性篩候選lncRNA描述性統(tǒng)長(zhǎng)度分布統(tǒng)外顯子數(shù)目統(tǒng)已知和預(yù)測(cè)lncRNA統(tǒng)lncRNA保守性分序列保守性分位點(diǎn)保守性分lncRNA差異表達(dá)lncRNA表達(dá)水平lncRNA差異表達(dá)差異表達(dá)lncRNA篩lncRNA組織或表型特異性分lncRNA與mRNA表達(dá)聚類組織或表型特異性分lncRNA靶預(yù)cis作用 預(yù)trans作用 預(yù)特異lncRNA靶功能富集分GO富集分KEGG富集分一、建庫(kù)流從N樣品到最終數(shù)據(jù)獲得，樣品檢測(cè)、建庫(kù)、每一個(gè)環(huán)節(jié)都會(huì)對(duì)數(shù)據(jù)質(zhì)量和數(shù)量產(chǎn)生影響，而數(shù)據(jù)質(zhì)量又會(huì)直接影響后續(xù)信息分析的結(jié)果。因此，獲得高質(zhì)量數(shù)據(jù)是保證生物信息分析正確、全面、可信的前提。為了從上保證數(shù)據(jù)的準(zhǔn)確性、可靠性，諾禾致源對(duì)樣品檢測(cè)、建庫(kù)、每一個(gè)生產(chǎn)步驟都嚴(yán)格把控，從根本上確保了高質(zhì)量數(shù)據(jù)的產(chǎn)出。實(shí)驗(yàn)流程圖如下：TotalRNA諾禾致源對(duì)RNA樣品的檢測(cè)主要包括4種方法瓊脂糖凝膠電泳分析RNA降解程度以及是否有污Nanodrop檢測(cè)RNA的純度（OD260/280比值Qubit對(duì)RNA濃度進(jìn)行精確定Agilent2100精確檢測(cè)RNA的完整文庫(kù)構(gòu)RNA檢測(cè)合格后，通過(guò)epicetreRb-ZroM試劑盒去除rRNA。隨后加入fragmnttinufer將RA打斷成短片段，以短片段RNA為模板，用六堿基隨機(jī)引物（randomheamrs）合成一鏈cDNA，然后加入緩沖液、dNTPs（dTP、dTP、dTP和dTP）和DNAplyerseI合成二鏈cDNA，隨后利用AMPuePeas純化雙鏈cDNA。純化的雙鏈cDNA再進(jìn)行末端修復(fù)、加A尾并連接接頭，然后用AMPureXPbeds進(jìn)行片段大小選擇。之后用RUcARNA庫(kù)文庫(kù)構(gòu)建完成后，先使用Qubit2.0進(jìn)行初步定量，稀釋文庫(kù)至1ng/ul，隨后使用Agilent2100對(duì)文庫(kù)的insertsize進(jìn)行檢測(cè)，insert符合預(yù)期后，使用Q-PCR方法對(duì)文庫(kù)的有效濃度進(jìn)行準(zhǔn)確定量（文庫(kù)有效濃度＞2nM），以保證文庫(kù)質(zhì)量上機(jī)庫(kù)檢合格后，把不同文庫(kù)按照有效濃度及目標(biāo)下機(jī)數(shù)據(jù)量的需求pooling后進(jìn)行HiSeq/MiSeq獲得原始序列(SequencedReads)后，在有相關(guān)物種參考序列或參考組的情況下，通過(guò)如程進(jìn)行生物信息分析1原始序列數(shù)高通量(如illuminaHiSeqTM2000/MiSeq等平臺(tái))得到的原始圖像數(shù)據(jù)文件經(jīng)堿基識(shí)別(BaseCalling)分析轉(zhuǎn)化為原始序列(SequencedReads)，我們稱之為RawData或RawReads，結(jié)果以FASTQ(簡(jiǎn)稱為fq)文件格式，其中包含序列(reads)的序列信息以及其對(duì)應(yīng)的質(zhì)量信息。FASTQ格式文件中每個(gè)read由四行描述，如@EAS139:136:FC706VJ:2:2104:15343:1973931:Y:18:ATCACG其中第一行以“@”開(kāi)頭，隨后為illumina識(shí)(SequeceIdntiies)和描述文字(選擇性部分)；第二行是堿基序列；第三行以“+”開(kāi)頭，隨后為illumna標(biāo)識(shí)符(選擇性部分)；第四行是對(duì)應(yīng)序列的質(zhì)量(Cocktl.)。illumina標(biāo)識(shí)符詳細(xì)信息如下UniqueinstrumentRunFlowcell2FlowcellTilenumberwithintheflowcell'x'-coordinateoftheclusterwithinthe'y'-coordinateoftheclusterwithinthe1Memberofapair,1or2(paired-endormate-pairreadsYYifthereadfailsfilter(readisbad),N0whennoneofthecontrolbitsareon,otherwiseitisanevenIndex第四行中每個(gè)字符對(duì)應(yīng)的ASCII值減去33，即為對(duì)應(yīng)第二行堿基的質(zhì)量值。如果錯(cuò)誤率用e表示，illumiaiSqT200/iSq的堿基質(zhì)量值用d表示，則有下列關(guān)系：公式一：Qphred10log10(e)illuminaCasava1.8版本錯(cuò)誤率與質(zhì)量值簡(jiǎn)明對(duì)應(yīng)關(guān)系如下：對(duì)應(yīng)字.5?0I2數(shù)據(jù)質(zhì)量評(píng) 錯(cuò)誤率分布檢每個(gè)堿基錯(cuò)誤率是通過(guò)Phred數(shù)值(Phredscore,Qphred)通過(guò)公式1轉(zhuǎn)化得到，而Phred數(shù)值是在堿基識(shí)別(BaseCalling)過(guò)程中通illuminaCasava1.8版本堿基識(shí)別與Phred分值之間的簡(jiǎn)明對(duì)應(yīng)關(guān)Phred分Q-錯(cuò)誤率與堿基質(zhì)量有關(guān)，受儀本身、試劑、樣品等多個(gè)因素共同影響。對(duì)于RNA-seq技術(shù)，錯(cuò)誤率分布具有兩個(gè)特點(diǎn)錯(cuò)誤率會(huì)隨著序列(SequecedRad)長(zhǎng)度的增加而升高，這是由于過(guò)程中化學(xué)試劑的消耗而導(dǎo)致的，并且為illumna高通量平臺(tái)都具有的特征。(2)前6個(gè)堿基的位置也會(huì)發(fā)生較高的錯(cuò)誤率，而這個(gè)長(zhǎng)度也正好等于在RNA-seq建庫(kù)過(guò)程中反轉(zhuǎn)錄所需要的隨機(jī)引物的長(zhǎng)度。所以推測(cè)前6個(gè)堿基錯(cuò)誤率較高的原因?yàn)殡S機(jī)引物和RNA模版的不完全結(jié)合(Jiangtl.)。圖2.1錯(cuò)誤率分布橫坐標(biāo)為reads的堿基位置，縱坐標(biāo)為單堿基錯(cuò)誤GCGC含量分布檢查用于檢測(cè)有無(wú)AT、GC分離現(xiàn)象，而這種現(xiàn)象可能是或者建庫(kù)所帶來(lái)的，并且會(huì)影響后續(xù)的定量分析在illuina平臺(tái)的轉(zhuǎn)錄組中，反轉(zhuǎn)錄成cDNA時(shí)所用的6bp會(huì)前置苷存定偏這偏好性與的物種和環(huán)境無(wú)關(guān)，但會(huì)影響轉(zhuǎn)錄組的均一化程度(Hanenetal)。除此之外，理論上G和C堿基及A和T堿基含量每個(gè)循環(huán)上應(yīng)分別相等，且整個(gè)過(guò)穩(wěn)定不變，呈水平線。對(duì)E來(lái)說(shuō)，由于隨機(jī)引物擴(kuò)增偏差等因，常常會(huì)導(dǎo)致在得的每個(gè)ed-7圖2.2GC含量分布橫坐標(biāo)為reads的堿基位置，縱坐標(biāo)為單堿基所占的比例；不同顏色代表不同的堿基 數(shù)據(jù)過(guò)得到的原始序列，里面含有帶接頭的、低質(zhì)量的reads，為了保證信息分析質(zhì)量，必須對(duì)rawreds進(jìn)行過(guò)濾，得到cleaneas，后續(xù)分析都基于cleanrads。數(shù)據(jù)處理的步驟如下去除帶接頭(adapter)的去除N(N表示無(wú)法確定堿基信息)的比例大于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ù)過(guò)濾結(jié)表2.4數(shù)據(jù)產(chǎn)出質(zhì)量情況一覽SampleRawCleanErrorGC00000030300000數(shù)據(jù)質(zhì)量情況詳細(xì)內(nèi)容如Rawreads：統(tǒng)計(jì)原始序列數(shù)據(jù)，以四行為一個(gè)單位，統(tǒng)計(jì)每個(gè)文件的序列的個(gè)數(shù)Cleanreads：計(jì)算方法同RawReads，只是統(tǒng)計(jì)的文件為過(guò)濾后的數(shù)據(jù)。后續(xù)的生物信息分析都是基于CleanreadsCleanbases：序列的個(gè)數(shù)乘以序列的長(zhǎng)度，并轉(zhuǎn)化為以G為單位Errorrate：通過(guò)公式1計(jì)算得到Q20、Q30：分別計(jì)算Phred數(shù)值大于20、30的堿基占總體堿基的百分比GCcontent：計(jì)算堿基G和C的數(shù)量總和占總的堿基數(shù)量的百分比3參考序列比對(duì)分我們采用Tophat2(Kimetal,2013)對(duì)過(guò)濾后的序列進(jìn)行參考組的比對(duì)分析.TopHat2的算法主要分為三個(gè)部分將序列和轉(zhuǎn)錄組進(jìn)行比對(duì)(可選將序列整段比對(duì)到組外顯子將序列分段比對(duì)到組的兩個(gè)外顯子上下圖為T(mén)opHat2的算法示意圖(Kimetal2013)：TopHat2的算法主要分為三個(gè)部分將序列和轉(zhuǎn)錄組進(jìn)行比對(duì)(可選將序列整段比對(duì)到組外顯子將序列分段比對(duì)到組的兩個(gè)外顯子如果參考組選擇合適，而且相關(guān)實(shí)驗(yàn)不存在污染，實(shí)驗(yàn)所產(chǎn)生的序列的定位的百分比正常情況下會(huì)高于70%(TotalMappedReadsorFragments)，其中具有多個(gè)定位的序列(MultipleMappedReadsorFragments)占總體的百分比通常不會(huì)超過(guò)10%。Reads與參 組比對(duì)情況統(tǒng)

表3.1Reads與參 組比對(duì)情況一覽SampleTotalTotal(80(82(80Multiple(14(23(14Uniquely(65(59Read-(33Read-(29Readsmapto(32(29Readsmapto'-(32Non-splice(47(41Splice(18(17(18ReadsmappedinproperProper-pairedreadsmaptodifferent526520582450504632比對(duì)結(jié)果統(tǒng)計(jì)詳細(xì)內(nèi)容如Totalreads：序列經(jīng)過(guò)數(shù)據(jù)過(guò)濾后的數(shù)量統(tǒng)計(jì)(Cleandata)Totlmaped：能定位到組上的序列的數(shù)量的統(tǒng)計(jì)；一般情況下，如果不存在污染并且參考組選擇合適的情況下，這部分?jǐn)?shù)據(jù)的百分比大于70%。Multiplemapped：在參考序列上有多個(gè)比對(duì)位置 序列的數(shù)量統(tǒng)計(jì)；這部分?jǐn)?shù)據(jù)的百分比一般會(huì)小于10%Uniquelymapped：在參考序列上有唯一比對(duì)位置 序列的數(shù)量統(tǒng)計(jì)Readsmapto'+'，Readsmapto'-'：序列比對(duì)到組上正鏈和負(fù)鏈的統(tǒng)計(jì)Splicereads：(2)中，分段比對(duì)到兩個(gè)外顯子上的 序列(也稱為Junctionreads)的統(tǒng)計(jì)，Non-splicereads為整段比對(duì)到外顯子 序列的統(tǒng)計(jì)，Splicereads的百分比取決于 Reads在參考組不同區(qū)域的分布情對(duì)Totalmppdreds的比對(duì)到組上的各個(gè)部分的情況進(jìn)行統(tǒng)計(jì)，定位區(qū)域分為Exon(外顯子)、Intro(內(nèi)含子)和Integeic(間隔區(qū)域)。正常情況下，Exon(外顯子)區(qū)域的序列定位的百分比含量應(yīng)該最高，定位到Intron(內(nèi)含子)區(qū)域的序列可能是由于非成mRNA的污染或者組注釋不完全導(dǎo)致的，而定位到Intergenic(間隔區(qū)域)的序列可能是因?yàn)榻M注釋不完全以及背景噪音圖3.2Reads在參 組不同區(qū)域的分布情Reads在上的密度分布情對(duì)Totalmappedreads的比對(duì)到組上的各個(gè)（分正負(fù)鏈）的密度進(jìn)行統(tǒng)計(jì)，如下圖所示，具體作圖的方法為用滑動(dòng)窗口(windowsize)為1K，計(jì)算窗口內(nèi)部比對(duì)到堿基位置上的reads的中位數(shù)，并轉(zhuǎn)化成log2。正常情況下，整個(gè)長(zhǎng)度越長(zhǎng)，該內(nèi)部定位的reads總數(shù)會(huì)越多(Marquezetal.)。從定位到上的reads數(shù)與長(zhǎng)度的關(guān)系圖中，可以更加直出長(zhǎng)度和reads總數(shù)的關(guān)圖3.3Reads在上的密度分布左圖：圖中最外圈是選擇展示的各條；中間的灰色背景區(qū)是抽取其中了10000reads的分布情況，紅色map到正鏈，藍(lán)色到負(fù)鏈；最里面的圓圈區(qū)是比對(duì)到該上的所有reads，橘黃色為正鏈coverage分布，綠色為負(fù)鏈coverage分布超過(guò)所有coverage集均值+3倍標(biāo)準(zhǔn)差的異點(diǎn)將被舍棄。圖：橫坐標(biāo)為的長(zhǎng)度信息(單位為Mb)，縱坐標(biāo)為mapped到上的reads數(shù)(單位為M)，圖中灰域表示95%的置信區(qū)間Reads比對(duì)結(jié)果IGV我們提供RNA-seqReds在組上比對(duì)結(jié)果的bam格式文件，部分物種還提供相應(yīng)的參考組和注釋文件，并推薦使用IGV(eteGenmisieer)瀏覽器對(duì)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瀏覽器使用方法可參考我們提供的使用說(shuō)明文檔(IGVQuikSar.pf)。圖3.4IGV瀏覽器界4表達(dá)分一個(gè)表達(dá)水平的直接體現(xiàn)就是其轉(zhuǎn)錄本的豐度情況，轉(zhuǎn)錄本豐度程度越高，則表達(dá)水平越高。在RNA-seq分析中，我們可以通過(guò)定位到組區(qū)域或外顯子區(qū)的序列(reads)的計(jì)數(shù)來(lái)估計(jì)的表達(dá)水平。通過(guò)不同Reads計(jì)數(shù)除了與的真實(shí)表達(dá)水平成正比外，還與的長(zhǎng)度和深度成正相關(guān)。為了使不同、不同實(shí)驗(yàn)間估計(jì)的表達(dá)水平具有可比性，人們引入了RPKM的概念，RPKM(RedsPerKilobsesprilionrad)是每百萬(wàn)reads中來(lái)自某一每千堿基長(zhǎng)度的reads數(shù)目。RPKM同時(shí)考慮了深度和長(zhǎng)度對(duì)read影響，是目前最為常用的表達(dá)水平估算方法(Mortzaitl.,208)。已知注釋類型含量分采用HTSeq軟件對(duì)各樣品不同進(jìn)行表達(dá)水平分析，使用的模型為union。根據(jù)表達(dá)量統(tǒng)計(jì)樣品中各類型的的表達(dá)分布，結(jié)果展示下圖4.1各類已 表達(dá)分布已知表達(dá)水平分分別統(tǒng)計(jì)不同表達(dá)水平下的數(shù)以及單個(gè)的表達(dá)水平一般況下，KM數(shù)1或1作為判斷是否達(dá)的閾值，不同的文獻(xiàn)所采用的閾值不同。表4.2.1不同表達(dá)水平區(qū)間的數(shù)量統(tǒng)計(jì)RPKM0-37473(7136140(691-6772(126786(135768(116712(123-6421(126441(126026(1115-2121(42222(4964(1表 表達(dá)水平統(tǒng)計(jì)10330000RNA-seq整體質(zhì)量評(píng)生物學(xué)重復(fù)是任何生物學(xué)實(shí)驗(yàn)所必須的，高通量技術(shù)也不例外(Hanseneta.)。生物學(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ì)劃建議皮爾遜相關(guān)系數(shù)的平方(R2)大于0.9(理想的取樣和實(shí)驗(yàn)條件下)。具體的項(xiàng)目操作中，我們要求R2至少要大于0.8，否則需要對(duì)樣品做出合適的解釋，或者重新進(jìn)行實(shí)驗(yàn)。圖5.1樣品間相關(guān)性檢左圖：樣品間的相關(guān)系數(shù)散點(diǎn)圖，R2:pearson相關(guān)系數(shù)的平方；右圖：樣品間相關(guān)系數(shù)熱樣品間聚類及PCA分當(dāng)樣本數(shù)目較多時(shí)>4），可利用的表達(dá)量進(jìn)行樣間聚類析及C分析，對(duì)樣本間關(guān)系進(jìn)探究或者對(duì)實(shí)驗(yàn)設(shè)計(jì)進(jìn)行證。CA圖5.2樣本間聚類及PCA分左圖：樣品的表達(dá)水平的層次聚類圖；右圖：樣品的表達(dá)水平的PCA聚類理想條件下，對(duì)于RNA-sq技術(shù)來(lái)說(shuō)，序列(reads)之間為獨(dú)立抽樣并且reads在所有表達(dá)的轉(zhuǎn)錄本上的分布應(yīng)該呈現(xiàn)均一化分布。然而很多研究表明，很多偏好型的因素都會(huì)影響這種均一化的分布(Dohmta.,208)。例如，在RNA-sq建庫(kù)過(guò)程中，片段破碎和RNA一樣會(huì)導(dǎo)致RNA-seq最終的數(shù)據(jù)呈現(xiàn)嚴(yán)重的3’偏。其他因素還包括轉(zhuǎn)錄區(qū)域的GC含量不同、隨機(jī)引物等等，并且生物體內(nèi)從5’或者3’的降解過(guò)程同樣會(huì)導(dǎo)致不均一性分布。圖5.3不同表達(dá)水平的轉(zhuǎn)錄本的reads密度分High：高表達(dá)量轉(zhuǎn)錄本；Medium：中度表達(dá)量轉(zhuǎn)錄本；Low：低表達(dá)量轉(zhuǎn)錄本；橫坐標(biāo)為距離轉(zhuǎn)錄本5’端的相對(duì)位置(以百分比表示)，縱坐標(biāo)為覆蓋深度的平均轉(zhuǎn)錄本拼采用Cufflinks(Trapnelletal,2013)和Scripture(Guttmanetal,2010)兩種軟件同時(shí)對(duì)比對(duì)結(jié)果進(jìn)行組裝，在此基礎(chǔ)上進(jìn)行l(wèi)ncRNA的選cufflinks拼uisfmcufflnk拼接結(jié)果展示如下：表6.1cuffliks拼接結(jié)果展示（部分ChrChrStartEnd1+gene_id"XLOC_000001";transcript_id";"1";oId"CUFF5.1";tss_id1+gene_id"XLOC_000001";transcript_id";"2";oId"CUFF5.1";tss_id1+gene_id"XLOC_000002";transcript_id";"1";oId"CUFF8.1";tss_id1+gene_id"XLOC_000002";transcript_id";"2";oId"CUFF8.1";tss_id1+gene_id"XLOC_000003";transcript_id";"1";oId"CUFF9.1";tss_id 表格說(shuō)明如下第1列：序號(hào)；第2列：來(lái)源描述；第3列：類型；第4列：起始坐標(biāo)；第5列：終止坐標(biāo)；第7列：鏈的信息；第9列：id等描述信息scripture拼表6.2scripture拼接結(jié)果展示（部分ChrStartEnd1chr1:14656--21chr1:16727--61chr1:24848--21chr1:135989-+31chr1:135989-+6表格說(shuō)明如下第1列：序號(hào)；第2列：起始坐標(biāo)；第3列：終止坐標(biāo)；第4列：轉(zhuǎn)錄本id；第6列：鏈的信息；第10列：exon個(gè)數(shù)；第11列：exon長(zhǎng)度；第12列：exon起始位置；候選lncRNA篩lncNA為一類長(zhǎng)度>200bp的長(zhǎng)鏈非編碼RNA，根據(jù)與編碼序列的位置關(guān)系可分為intergeiclnRN(簡(jiǎn)稱linRN)，inroiclnRN,ant-snselcRA，seselnRNA，idretinallcRA等類型。其中l(wèi)incRA所占比例最高，這里主要進(jìn)行前3種類型的篩選。我們根據(jù)lncNA的特點(diǎn)設(shè)置一系列嚴(yán)格的篩選條件，基于cuffliks和sritue的拼接結(jié)果同時(shí)進(jìn)行以下步驟的篩選，最終選擇在>=2lncRNA作為最終的候選lncRNA集進(jìn)行后續(xù)分析?；竞Y基本篩選主要由三個(gè)部分組成step1:選擇長(zhǎng)度>=200bp，Exon個(gè)數(shù)>=2的轉(zhuǎn)錄本step2:通過(guò)cufflinks計(jì)算每條轉(zhuǎn)錄本的reads覆蓋度，選擇Reads最小覆蓋度>=3的轉(zhuǎn)錄本step3:通過(guò)與已知非lncRNA比較過(guò)濾掉非lncRNA，并利用 下圖展示的是lncRNA的篩選過(guò)程統(tǒng)計(jì)圖7.1.1lncRNA的篩選統(tǒng)計(jì)縱坐標(biāo)為篩選步驟（C代表cufflinks，S代表scripture，assembly為原始拼接出來(lái)的轉(zhuǎn)錄本條數(shù)），橫坐標(biāo)為對(duì)應(yīng)步驟篩選過(guò)后的轉(zhuǎn)錄本條

圖7.1.2lncRNA的篩選統(tǒng)計(jì)橫坐標(biāo)為各class_code類型，縱坐標(biāo)為對(duì)應(yīng)類型的轉(zhuǎn)錄本條數(shù)；左圖為cufflinks的結(jié)果，右圖為scripture的結(jié)諾禾致源生物信息科技具有編碼潛能與否是判斷轉(zhuǎn)錄本是否為lncRNA的關(guān)鍵條件，我們綜合了目前主流的編碼潛能分析方法進(jìn)行該項(xiàng)篩選，主要包括：CPCCNCI分析、pfam蛋白結(jié)構(gòu)域分析、PhyloCSF分析四種方法CPC分CPC(CodingPotentialCalculator)結(jié)果展示如下-0CNCI分CNCI(Coding-Non-CodingIndex)結(jié)果展示如下score:-0start:stop:start:stop:score:-0start:stop:start:stop:start:stop:start:stop:pfampfam蛋白結(jié)構(gòu)域搜索結(jié)果展示如下hmmhmmE-0-77135e-10-PF1281324.4e-10-2210-332e-10-PF00439298e-1(3)phyloCSF分phyoCF(hyogneiccoonsusttuionfeqeny)進(jìn)化子置換頻率分析，利用多物種間的全組序列比對(duì)文件定義一段組區(qū)域是否有編碼潛能。通過(guò)文獻(xiàn)查詢，我們發(fā)現(xiàn)不同的物種間phyloCSF闕值不盡相同，故首先隨機(jī)選擇本項(xiàng)目研究物種一定數(shù)目的已知og和lncRA進(jìn)行闕值分析，再篩選候選轉(zhuǎn)錄本分析結(jié)果。phyloS結(jié)果展示如下：圖7.2.1lncRNA的篩選統(tǒng)計(jì)橫坐標(biāo)為phyloCSF的分值，縱坐標(biāo)為對(duì)應(yīng)分值的轉(zhuǎn)錄本占所有轉(zhuǎn)錄本條數(shù)的比將4種軟件的結(jié)果取交集圖7.2.24種方法結(jié)果維恩圖展左邊為cuffliks的結(jié)果，右邊為scripture的結(jié)2cAcAnA橫坐標(biāo)為lncRNA類型，縱坐標(biāo)為對(duì)應(yīng)類型的轉(zhuǎn)錄本條候選lncRNA對(duì)篩選得到的lncRNA進(jìn)行長(zhǎng)度，exon個(gè)數(shù)等方面的統(tǒng)計(jì)，有助于進(jìn)一步觀察篩選得到的候選lncRNA的特點(diǎn)，并通過(guò)與本物種已知lncRNA得到已知lncRNA和新預(yù)測(cè)的lncRNA。圖8.1lncRNA長(zhǎng)度分布橫坐標(biāo)為lncRNA長(zhǎng)度(bp)，縱坐標(biāo)為對(duì)應(yīng)長(zhǎng)度的轉(zhuǎn)錄本密橫坐標(biāo)為外顯子個(gè)數(shù)，縱坐標(biāo)為對(duì)應(yīng)轉(zhuǎn)錄本的數(shù)已知和預(yù)測(cè)lncRNA統(tǒng)

圖8.3知和預(yù)測(cè)lncRNA維恩lncRNA保守性分lncNA的序列保守性相對(duì)蛋白編碼要低，采用phastCos( ge.bcbconel.duphst)分別對(duì)蛋白編碼和lncRNA基因進(jìn)行保守性打分，得到如下保守性分值累積分布圖。圖9.3lncRNA和蛋白編 的保守性分值累積分布lncRNA的序列在物種間有一定的位點(diǎn)保守性，通過(guò)UCSC瀏覽器可視化lncRNA在不同物種中的位置圖9.3lncRNA和蛋白編 的保守性分值累積分布lncRNA差異表達(dá)分lncRNA通過(guò)所有l(wèi)ncRNA的RPM的分布圖以及盒形圖對(duì)不同實(shí)驗(yàn)條件下的lncRNA表達(dá)水平進(jìn)行比較。對(duì)于同一實(shí)驗(yàn)條件下的重復(fù)樣品，最終的RPKM為所有重復(fù)數(shù)據(jù)的平均值。圖10.1不同實(shí)驗(yàn)條件下lncRNA表達(dá)水平比對(duì)圖一：RPKM盒形圖，橫坐標(biāo)為樣品名稱，縱坐標(biāo)為log1(RPKMRPKM分布圖，橫坐標(biāo)為log10(RPKM),縱坐標(biāo)為的密度lncRNAlncRNA差異表達(dá)的輸入數(shù)據(jù)為lncRNA表達(dá)水平分析中得到的readcount數(shù)據(jù)。對(duì)于有生物學(xué)重復(fù)的樣品，分析我們采用DESeq（Andersetal,2010）進(jìn)行分析：該分析方法基于的模型是負(fù)二項(xiàng)分布，第i個(gè)在第j個(gè)樣本中的readcount值為Kij，則Kij～對(duì)于無(wú)生物學(xué)重復(fù)的樣品，先采用TMM對(duì)readcount數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，之后用DEGseq進(jìn)行差異分析。差異表達(dá)列表如下表9.2差異列-32.7547e-13326e-1431.0235e-20042e--51.2429e-20042e-0-69.0982e-9.7672e-差異lncRNA列表主要包括的內(nèi)容Gene_id:編readcount_Sample1：校正后樣品組1的readcountreadcount_Sample2：校正后樣品組2的readcountlog2FoldChange:pvalue(pval):統(tǒng)計(jì)學(xué)差異顯著性檢驗(yàn)指qvalue(padj):校正后的pvalue。qvalue越小，表 表達(dá)差異越顯差異表達(dá)lncRNA篩lNAlncRNA進(jìn)行篩選，閾值設(shè)定一般為:|lo2FodCane)|>1且qvlue<0.05。對(duì)于有生物學(xué)重復(fù)的實(shí)驗(yàn)，由于DESeq已經(jīng)進(jìn)行了生物學(xué)變異的消除，我們對(duì)差異lncRNA篩選的標(biāo)準(zhǔn)一般為:padj<0.05。圖10.3差異lncRNA火山有顯著性差異表達(dá)的lncRNA用紅色點(diǎn)表示；橫坐標(biāo)代表lncRNA在不同樣本中表達(dá)倍數(shù)變化；縱坐標(biāo)代表lncRNA表達(dá)量變化差異的統(tǒng)計(jì)學(xué)顯著lncRNAlncRNA與mRNA通常認(rèn)為lincRNA相對(duì)于mRNA有較高的組織表達(dá)特異性，隨機(jī)抽取一定比例的lncRNA和mNA比較兩種類型在不同組織中表達(dá)水平的聚類情況。圖11.llncRNA和mRNA的表達(dá)熱左圖為lncRNA在各樣品中的表達(dá)情況，右圖為mRNA在個(gè)樣品中的表達(dá)情況；橫坐標(biāo)為樣品，縱坐標(biāo)為，顏色越深表示表達(dá)水平越組織或表型特異性分我們基于JSdvegece這一衡量指標(biāo)對(duì)于各轉(zhuǎn)錄本在不同組織樣本中的表達(dá)模式（pattern）進(jìn)行分析。參考文獻(xiàn)（Cabil,MN.tal.，01），預(yù)先設(shè)定每個(gè)轉(zhuǎn)錄本僅在中特異性表達(dá)有N（N為組織個(gè)數(shù)）種模式，定義每?jī)蓚€(gè)轉(zhuǎn)錄本表達(dá)模式之間的距離為S的平方根，則一個(gè)轉(zhuǎn)錄本在N個(gè)組織中的組織特異性定義為:其中et為預(yù)先設(shè)定的轉(zhuǎn)錄本表達(dá)模式。選擇其中最大的值作為該轉(zhuǎn)錄本在各組織中的特異性分值。分值范圍為0-1，分?jǐn)?shù)越接近于1，表示該轉(zhuǎn)錄本的組織特異性越高。圖11.2轉(zhuǎn)錄本的組織特異性分值密度分布橫坐標(biāo)為轉(zhuǎn)錄本的組織特異性分值；縱坐標(biāo)為轉(zhuǎn)錄本的密lncRNA靶預(yù)lncRNA功能主要通過(guò)cis或trans作用于蛋白編碼靶的方式實(shí)現(xiàn)，因此分成兩種情況預(yù)測(cè)lncRNA的靶cis作用靶預(yù)is功預(yù)測(cè)基本原理認(rèn)nA的功能與其坐標(biāo)的編碼蛋白相關(guān)，于是將RA~上下k0)蛋白編碼的找出進(jìn)行功能富集分析，以推測(cè)lncRNA的主要功能。cis作用靶預(yù)測(cè)結(jié)果如下表所示：表12.1cis作用靶統(tǒng)計(jì)123459注Expressed_Samplenum：lncRNA在n個(gè)樣品中表達(dá)lncRNA_num：lncRNA在n個(gè)樣品中表達(dá)的數(shù)目10kb(lncRNA/mRNA)：在上下游10kb范圍內(nèi)(2)中的lncRNA能檢測(cè)到的mRNA的數(shù)目100kb(lncRNA/mRNA)：在上下游100kb范圍內(nèi)(2)中的lncRNA能檢測(cè)到的mRNA的數(shù)目trans作用靶預(yù)trans功能預(yù)測(cè)基本原理認(rèn)為lncRNA的功能與樣品中共表達(dá)的編碼蛋白相關(guān)，可以通過(guò)樣本間lncRNA與蛋白編碼的相關(guān)性分析或共表達(dá)分析來(lái)預(yù)測(cè)。當(dāng)樣本量>=5時(shí)采用Pearson相關(guān)系數(shù)法分析樣本間lncRNA與蛋白編碼的相關(guān)性；當(dāng)樣本數(shù)> etal,2008)將不同的組織、處理或者時(shí)間點(diǎn)間表達(dá)模式相似的聚類，以得到不同的共表達(dá)模塊，根據(jù)模塊內(nèi)已知的蛋白編碼的功能進(jìn)圖9.6共表達(dá)聚類熱 特異lncRNA靶功能富集分特異lncRNA一般指差異表達(dá)的或者組織或表型特異性表達(dá)的lncRNA，對(duì)這些lncRNA對(duì)應(yīng)的靶分別進(jìn)行GO和KEGG功能富集分析GO富集分GeneOntology（簡(jiǎn)稱GO,）是功能國(guó)際標(biāo)準(zhǔn)分類體系。根據(jù)實(shí)驗(yàn)?zāi)康暮Y選特定lncRNA后，研究該lncRNA對(duì)應(yīng)的靶在GeneOntology中的分布狀況將闡明實(shí)驗(yàn)中樣本差異在功能上的體現(xiàn)。GO富集分析方法為GOseq（Youngetal,2010）,此方法基于Walleniusnon-centralhyper-geometricdistribution。相對(duì)于普通的超幾何分布(Hyper-geometricdistribution)，此分布的特點(diǎn)是從某個(gè)類別中抽取的概率與從某個(gè)類別之外抽取一個(gè)的概率是不同的，這種概率的不同是通過(guò)對(duì)長(zhǎng)度的偏進(jìn)行估計(jì)得到的，從而能更為準(zhǔn)確地計(jì)算出GOterm被靶富集的概率。表13.1.1靶GO富集列sequence-specificDNAbindingtranscriptionfactor3.7832e-08nucleicacidbindingtranscriptionfactor3846e-08regulationofgene3.4857e-0nucleobase-containingcompoundbiosynthetic5.7706e-0結(jié)果表格詳細(xì)內(nèi)容如下GO_accession：GeneOntology數(shù)據(jù)庫(kù)中唯一的標(biāo)號(hào)信Description：GeneOntology功能的描述信Term_type：該GO的類別 ponent：細(xì)胞組分；biological_process：生物學(xué)過(guò)程；molecular_function：分子功能Over_represented_pValue：富集分析統(tǒng)計(jì)學(xué)顯著水Corrected_pValue：矯正后的P-Value，一般情況下，P-value0.05功能為富集DEG_item：與該GO相關(guān)的靶的數(shù)DEG_list：GO注釋的靶數(shù)有向無(wú)環(huán)圖(DirectdcylicGaph，DG)為差異GO富集分析結(jié)果的圖形化展示方式，分支代表包含關(guān)系，從上至下所定義的功能范圍越來(lái)越小，一般選取GO富集分析的結(jié)果前10位作為有向無(wú)環(huán)圖的主節(jié)點(diǎn)，并通過(guò)包含關(guān)系，將相關(guān)聯(lián)的GOTerm一起展示，顏色的深淺代表富集程度。我們的項(xiàng)目中分別繪制生物過(guò)程(bioloiclroes)、分子功能(moleclarfncio)和細(xì)胞組分(celluarcopoen)的DA圖13.1.2GO富集有向無(wú)環(huán)每個(gè)節(jié)點(diǎn)代表一個(gè)GO術(shù)語(yǔ)，方框代表的是富集程度為T(mén)OP10GO，顏色的深淺代表富集程度，顏色越深就表示富集程度越高，每個(gè)節(jié)點(diǎn)上展示了該TERM析的p-value靶GO富集柱狀圖，直觀的反映出在生物過(guò)程(biologicalprocess)、細(xì)胞組分(cellularcomponent)和分子功能(molecularfunction)富集的GOterm上靶的個(gè)數(shù)分布情況。我們挑選了富集最顯著的30個(gè)GOterm在圖中展示，如果不足30條，則全部展示。圖13.1.3GO富集柱狀每組兩張圖；左圖：縱坐標(biāo)為富集的GOterm，橫坐標(biāo)為該term中靶個(gè)數(shù)。不同顏色用來(lái)區(qū)分生物過(guò)程、細(xì)胞組分和分子功能，帶“*”為富集的GOterm右圖：對(duì)KEGG富集分

諾禾致源生物信息科技在生物體內(nèi)，不同相互協(xié)調(diào)行使其生物學(xué)功能，通過(guò)Pathway顯著性富集能確定靶參與的最主要生化代謝途徑和信號(hào)轉(zhuǎn)導(dǎo)途徑。KEGG（KyotoEncyclopediaofGenesandGenomes）是有關(guān)Pathway的主要公共數(shù)據(jù)庫(kù)(Kanehisa,2008）。Pathway顯著性富集分析以KEGGPathway為單位，應(yīng)用超幾何檢驗(yàn)，找出與整個(gè)組背景相比，在靶中顯著性富集的Pathway。該分析的計(jì)算公式:在這里N為所有中具有Pathway注釋的數(shù)目；n為N中差異表達(dá)的數(shù)目；M為所有中注釋為某特定Pathway的數(shù)目；m為注釋為某特定Pathway的差異表達(dá)數(shù)目。FDR≤0.05的Pathway定義為在差異表達(dá)中顯著富集的Pathway，我們使用KOBAS（2.0）進(jìn)行表13.2.1差異KEGG富集列SampleBackgroundP-CorrectedP-GlycosphingolipidKEGG10TGF-betasignalingKEGG1AxonKEGG1MetabolicKEGG1結(jié)果表格詳細(xì)內(nèi)容如下#Term：KEGG通路的描述信息Id：KEGG數(shù)據(jù)庫(kù)中通路唯一的編號(hào)信息Samplenumber：該通路下靶的個(gè)數(shù)Backgroundnumber：該通路下的個(gè)數(shù)P-value：富集分析統(tǒng)計(jì)學(xué)顯著水平CorrectedP-value：矯正后的統(tǒng)計(jì)學(xué)顯著水平，一般情況下，P-value0.05功能為富集項(xiàng)散點(diǎn)圖是KEGG富集分析結(jié)果的圖形化展示方式。在此圖中，KEGG富集程度通過(guò)Richfator、vaue和富集到此通的個(gè)數(shù)來(lái)衡量。其中Richacor指差異表達(dá)的中位于該pathway條目的數(shù)目與所有有注釋中位于該pathway條目的總數(shù)的比值。Richfator越大，表示富集的程度越大。Qvalue是做過(guò)多重假設(shè)檢驗(yàn)校正之后的Pvalue，Qale的取值范圍為[0,1]，越接近于零，表示富集越顯著。我們挑選了富集最顯著的20條pathay條目在該圖中進(jìn)行展示，若富集的pathway條目不足20條，則全部展示。圖13.2.2富集的KEGG代謝通路的散點(diǎn)縱軸表示pathway名稱，橫軸表示 factor，點(diǎn)的大小表示此pathway中差異表達(dá)個(gè)數(shù)多少，而點(diǎn)的顏色對(duì)應(yīng)于不同的Qvalue范圍將差異富集出的通路圖展示出，該通路圖中，包含上調(diào)的節(jié)點(diǎn)標(biāo)紅色，包含下調(diào)O節(jié)點(diǎn)標(biāo)綠色，包含下調(diào)的標(biāo)黃色。鼠標(biāo)懸停于標(biāo)記的KO節(jié)點(diǎn)，彈出差異細(xì)節(jié)框，標(biāo)色同上，括號(hào)中數(shù)字為log2(Fldchng)。以上步驟可脫機(jī)實(shí)現(xiàn)，如連接互聯(lián)網(wǎng)，點(diǎn)擊各個(gè)節(jié)點(diǎn)，可以連接KG數(shù)據(jù)庫(kù)中個(gè)O的具體信頁(yè)。圖13.2.3著富集的KEGGpathway代謝通路KEGG代謝通路圖中，包含上調(diào)的KO節(jié)點(diǎn)標(biāo)紅色，包含下調(diào)的KO節(jié)點(diǎn)標(biāo)綠色，包含上下調(diào)的標(biāo)黃色。在網(wǎng)頁(yè)界面上，鼠標(biāo)懸停于標(biāo)記的KO節(jié)點(diǎn)，出差細(xì)節(jié)框，標(biāo)色同上，括號(hào)中數(shù)字為log2(Foldchange)。特異lncRNA與mRNAlncNA與mNA可以通過(guò)靶向關(guān)系進(jìn)行關(guān)聯(lián)，mRNA和mRNA之間可以通過(guò)蛋白質(zhì)互作關(guān)系進(jìn)行關(guān)聯(lián)，從而可以形成lncRNAmRA-roein網(wǎng)絡(luò)關(guān)系。mRNA和RNA之間主要應(yīng)用STRING蛋白質(zhì)互作數(shù)據(jù)庫(kù)()中的互作關(guān)系，針對(duì)數(shù)據(jù)庫(kù)中包含的物種，直接從數(shù)據(jù)庫(kù)中提取出目標(biāo)集(比如差異list)的互作關(guān)系構(gòu)建網(wǎng)絡(luò)。我們提供特異lncRNA與靶，靶蛋白互作網(wǎng)絡(luò)數(shù)據(jù)文件，此文件可以直接導(dǎo)入Cytoscae軟件進(jìn)行可視化編輯。Cytosape用方法可參考我們提供的使用說(shuō)明文檔(CytosapQucktat.df)?？蛻艨梢葬槍?duì)一些網(wǎng)絡(luò)的拓?fù)鋵傩赃M(jìn)行統(tǒng)計(jì)和標(biāo)示作圖，比如：互作網(wǎng)絡(luò)圖中節(jié)點(diǎn)(node)的大小與此節(jié)點(diǎn)的度(degre)成正比，即與此節(jié)點(diǎn)相連的邊越多，它的度越大，節(jié)點(diǎn)也就越大，這些節(jié)點(diǎn)在網(wǎng)絡(luò)中可能處于較為的位置。節(jié)點(diǎn)的顏色與此節(jié)點(diǎn)的系數(shù)(clustrigoefiiet相關(guān)，顏色梯度由綠到紅對(duì)應(yīng)系數(shù)的值由低到高；系數(shù)表示此節(jié)點(diǎn)的鄰接點(diǎn)之間的連通性好壞，系數(shù)值越高表示此節(jié)點(diǎn)的鄰接點(diǎn)之間的連通性越好等等。根據(jù)不同的研究目的和需求，客戶還可以在網(wǎng)絡(luò)圖中進(jìn)行調(diào)整節(jié)點(diǎn)位置和顏色、標(biāo)注表達(dá)量水操作。需要注意的是，通過(guò)blast比對(duì)得到的結(jié)果不能保證較好的準(zhǔn)確性，這部分的工作只是給客戶提供參考，輔助客戶發(fā)現(xiàn)一些可能的重要的。按我們提供的使用說(shuō)明將文件導(dǎo)入Cytosape軟件后的效果圖如下：圖14Cytoscape軟件界Anders,S.(2010).HTSeq:ysinghigh-throughputsequencingdatawithAnders,S.,andHuber,W.(2010).Differentialexpressionysisforsequencecountdata.GenomeBiol.(DESeq)Anders,S.andHuber,W.(2012).DifferentialexpressionofRNA-Seqdataatthegenelevel-theDESeqpackage.(DEseq)AndersS,ReyesA,HuberW.(2012).DetectingdifferentialusageofexonsfromRNA-seqdata.GenomeResearch.(DEXSeq)Kanehisa,M.,M.Araki,etal.(2008).KEGGforlinkinggenomestolifeandtheenvironment.Nucleicacidsresearch.(KEGG)Kim,D.,G.Pertea,etal.(2012).TopHat2:ParallelmapoftranscriptomestodetectInDels,genefusions,andLangmead,B.,Trapnell,C.,Pop,M.&Salz