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

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

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

圖7.1.2lncRNA的篩選統(tǒng)計橫坐標(biāo)為各class_code類型，縱坐標(biāo)為對應(yīng)類型的轉(zhuǎn)錄本條數(shù)；左圖為cufflinks的結(jié)果，右圖為scripture的結(jié)諾禾致源生物信息科技具有編碼潛能與否是判斷轉(zhuǎn)錄本是否為lncRNA的關(guān)鍵條件，我們綜合了目前主流的編碼潛能分析方法進行該項篩選，主要包括：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)進化子置換頻率分析，利用多物種間的全組序列比對文件定義一段組區(qū)域是否有編碼潛能。通過文獻查詢，我們發(fā)現(xiàn)不同的物種間phyloCSF闕值不盡相同，故首先隨機選擇本項目研究物種一定數(shù)目的已知og和lncRA進行闕值分析，再篩選候選轉(zhuǎn)錄本分析結(jié)果。phyloS結(jié)果展示如下：圖7.2.1lncRNA的篩選統(tǒng)計橫坐標(biāo)為phyloCSF的分值，縱坐標(biāo)為對應(yīng)分值的轉(zhuǎn)錄本占所有轉(zhuǎn)錄本條數(shù)的比將4種軟件的結(jié)果取交集圖7.2.24種方法結(jié)果維恩圖展左邊為cuffliks的結(jié)果，右邊為scripture的結(jié)2cAcAnA橫坐標(biāo)為lncRNA類型，縱坐標(biāo)為對應(yīng)類型的轉(zhuǎn)錄本條候選lncRNA對篩選得到的lncRNA進行長度，exon個數(shù)等方面的統(tǒng)計，有助于進一步觀察篩選得到的候選lncRNA的特點，并通過與本物種已知lncRNA得到已知lncRNA和新預(yù)測的lncRNA。圖8.1lncRNA長度分布橫坐標(biāo)為lncRNA長度(bp)，縱坐標(biāo)為對應(yīng)長度的轉(zhuǎn)錄本密橫坐標(biāo)為外顯子個數(shù)，縱坐標(biāo)為對應(yīng)轉(zhuǎn)錄本的數(shù)已知和預(yù)測lncRNA統(tǒng)

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

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