有參考基因組的轉(zhuǎn)錄組生物信息分析模板

1、一、生物信息分析流程獲得原始測序序列(Sequenced Reads)后，在有相關(guān)物種參考序列或參考基因組的情況下，通過如下流程進行生物信息分析：二、項目結(jié)果說明1原始序列數(shù)據(jù)項目結(jié)果文件高通量測序(如illumina HiSeqTM2000/MiSeq等測序平臺)測序得到的原始圖像數(shù)據(jù)文件經(jīng)堿基識別(Base Calling)分析轉(zhuǎn)化為原始測序序列(Sequenced Reads)，我們稱之為Raw Data或Raw Reads，結(jié)果以FASTQ(簡稱為fq)文件格式存儲，其中包含測序序列(reads)的序列信息以及其對應(yīng)的測序質(zhì)量信息。測序樣品真實數(shù)據(jù)隨機截取展示FASTQ格式文件中每個r

2、ead由四行描述，如下：EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG GCTCTTTGCCCTTCTCGTCGAAAATTGTCTCCTCATTCGAAACTTCTCTGT + CFFFDEHHHHFIJJJFHGIIIEHIIJBHHHIJJEGIIJJIGHIGHCCF 其中第一行以“”開頭，隨后為illumina 測序標識符(Sequence Identifiers)和描述文字(選擇性部分)；第二行是堿基序列；第三行以“+”開頭，隨后為illumina 測序標識符(選擇性部分)；第四行是對應(yīng)序列的測序質(zhì)量(Cock et al.

3、)。illumina 測序標識符詳細信息如下：EAS139Unique instrument name136Run IDFC706VJFlowcell ID2Flowcell lane2104Tile number within the flowcell lane15343'x'-coordinate of the cluster within the tile197393'y'-coordinate of the cluster within the tile1Member of a pair, 1 or 2 (paired-end or mate-pair

4、reads only)YY if the read fails filter (read is bad), N otherwise180 when none of the control bits are on, otherwise it is an even numberATCACGIndex sequence第四行中每個字符對應(yīng)的ASCII值減去33，即為對應(yīng)第二行堿基的測序質(zhì)量值。如果測序錯誤率用e表示，illumina HiSeqTM2000/MiSeq的堿基質(zhì)量值用Qphred表示，則有下列關(guān)系：公式一：Qphred=-10log10(e)illumina Casava 1.8版本測

5、序錯誤率與測序質(zhì)量值簡明對應(yīng)關(guān)系如下：測序錯誤率測序質(zhì)量值對應(yīng)字符5%13.1%2050.1%30?0.01%40I2測序數(shù)據(jù)質(zhì)量評估項目結(jié)果文件2.1測序錯誤率分布檢查每個堿基測序錯誤率是通過測序Phred數(shù)值(Phred score, Qphred)通過公式1轉(zhuǎn)化得到，而Phred 數(shù)值是在堿基識別(Base Calling)過程中通過一種預(yù)測堿基判別發(fā)生錯誤概率模型計算得到的，對應(yīng)關(guān)系如下表所顯示：illumina Casava 1.8版本堿基識別與Phred分值之間的簡明對應(yīng)關(guān)系Phred分值不正確的堿基識別堿基正確識別率Q-sorce101/1090%Q10201/10099%Q20

6、301/100099.9%Q30401/1000099.99%Q40測序錯誤率與堿基質(zhì)量有關(guān)，受測序儀本身、測序試劑、樣品等多個因素共同影響。對于RNA-seq技術(shù)，測序錯誤率分布具有兩個特點：(1)測序錯誤率會隨著測序序列(Sequenced Reads)的長度的增加而升高，這是由于測序過程中化學(xué)試劑的消耗而導(dǎo)致的，并且為illumina高通量測序平臺都具有的特征(Erlich and Mitra, 2008; Jiang et al.)。(2)前6個堿基的位置也會發(fā)生較高的測序錯誤率，而這個長度也正好等于在RNA-seq建庫過程中反轉(zhuǎn)錄所需要的隨機引物的長度。所以推測前6個堿基測序錯誤率較

7、高的原因為隨機引物和RNA模版的不完全結(jié)合(Jiang et al.)。測序錯誤率分布檢查用于檢測在測序長度范圍內(nèi)，有無異常的堿基位置存在高錯誤率，比如中間位置的堿基測序錯誤率顯著高于其他位置。一般情況下，每個堿基位置的測序錯誤率都應(yīng)該低于0.5%。圖2.1測序錯誤率分布圖橫坐標為reads的堿基位置，縱坐標為單堿基錯誤率2.2GC含量分布檢查GC含量分布檢查用于檢測有無AT、GC 分離現(xiàn)象，而這種現(xiàn)象可能是測序或者建庫所帶來的，并且會影響后續(xù)的定量分析。在illumina測序平臺的轉(zhuǎn)錄組測序中，反轉(zhuǎn)錄成cDNA時所用的6bp 的隨機引物會引起前幾個位置的核苷酸組成存在一定的偏好性。而這種偏好

8、性與測序的物種和實驗室環(huán)境無關(guān)，但會影響轉(zhuǎn)錄組測序的均一化程度(Hansen et al.)。除此之外，理論上G和C堿基及A和T堿基含量每個測序循環(huán)上應(yīng)分別相等，且整個測序過程穩(wěn)定不變，呈水平線。對于DGE測序來說，由于隨機引物擴增偏差等原因，常常會導(dǎo)致在測序得到的每個read前6-7個堿基有較大的波動，這種波動屬于正常情況。圖2.2GC含量分布圖橫坐標為reads的堿基位置，縱坐標為單堿基所占的比例；不同顏色代表不同的堿基類型2.3測序數(shù)據(jù)過濾測序得到的原始測序序列，里面含有帶接頭的、低質(zhì)量的reads，為了保證信息分析質(zhì)量，必須對raw reads進行過濾，得到clean reads，后續(xù)

9、分析都基于clean reads。數(shù)據(jù)處理的步驟如下：(1) 去除帶接頭(adapter)的reads；(2) 去除N(N表示無法確定堿基信息)的比例大于10%的reads；(3) 去除低質(zhì)量reads。RNA-seq 的接頭(Adapter, Oligonucleotide sequences for TruSeqTM RNA and DNA Sample Prep Kits) 信息：RNA 5 Adapter (RA5), part # 15013205：5-AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT-3RNA 3

10、Adapter (RA3), part # 15013207：5-GATCGGAAGAGCACACGTCTGAACTCCAGTCAC(6位index)ATCTCGTATGCCGTCTTCTGCTTG-3圖2.3原始數(shù)據(jù)過濾結(jié)果 2.4測序數(shù)據(jù)質(zhì)量情況匯總表2.4數(shù)據(jù)產(chǎn)出質(zhì)量情況一覽表Sample nameRaw readsClean readsclean basesError rate(%)Q20(%)Q30(%)GC content(%)HS1_136579608351752053.52G0.0397.8892.8849.39HS1_236579608351752053.52G0.

11、0396.5090.3849.59HS2_136547734351194633.51G0.0397.8592.8149.53數(shù)據(jù)質(zhì)量情況詳細內(nèi)容如下：(1) Raw reads：統(tǒng)計原始序列數(shù)據(jù)，以四行為一個單位，統(tǒng)計每個文件的測序序列的個數(shù)。(2) Clean reads：計算方法同 Raw Reads，只是統(tǒng)計的文件為過濾后的測序數(shù)據(jù)。后續(xù)的生物信息分析都是基于Clean reads。(3) Clean bases：測序序列的個數(shù)乘以測序序列的長度，并轉(zhuǎn)化為以G為單位。(4) Error rate：通過公式1計算得到。(5) Q20、Q30：分別計算 Phred 數(shù)值大于20、30的堿基占

12、總體堿基的百分比。(6) GC content：計算堿基G和C的數(shù)量總和占總的堿基數(shù)量的百分比。3參考序列比對分析項目結(jié)果文件測序序列定位算法：根據(jù)不同的基因組的特征，我們選取相對合適的軟件(動植物用TopHat(Trapnell et al., 2009)、真菌或者基因密度較高的物種用Bowtie)，合適的參數(shù)設(shè)置(如最大的內(nèi)含子長度，會根據(jù)已知的該物種的基因模型來進行統(tǒng)計分析)，將過濾后的測序序列進行基因組定位分析。下圖為TopHat的算法示意圖： Tophat的算法主要分為兩個部分：(1) 將測序序列整段比對到外顯子上。(2) 將測序序列分段比對到兩個外顯子上。我們統(tǒng)計了實驗所產(chǎn)生的測序

13、序列的定位個數(shù)(Total Mapped Reads)及其占clean reads的百分比，其中包括多個定位的測序序列個數(shù)(Multiple Mapped Reads)及其占總體（clean reads）的百分比，以及單個定位的測序序列個數(shù)(Uniquely Mapped Reads)及其占總體（clean reads）的百分比。3.1Reads與參考基因組比對情況統(tǒng)計表3.1Reads與參考基因組比對情況一覽表Sample nameHS1HS2HT1HT2HW1HW2Total reads703504107023892676161678506660844657366240543118Tota

14、l mapped60529821 (86.04%)60232484 (85.75%)63555439 (83.45%)43461327 (85.78%)40246848 (86.42%)34971284 (86.26%)Multiple mapped606556 (0.86%)633575 (0.9%)714678 (0.94%)450156 (0.89%)389470 (0.84%)335509 (0.83%)Uniquely mapped59923265 (85.18%)59598909 (84.85%)62840761 (82.51%)43011171 (84.89%)39857378

15、(85.58%)34635775 (85.37%)Read-130176973 (42.9%)29987004 (42.69%)31592931 (41.48%)21654629 (42.74%)20028779 (43%)17411209 (43.02%)Read-229746292 (42.28%)29611905 (42.16%)31247830 (41.03%)21356542 (42.15%)19828599 (42.57%)17224566 (42.35%)Reads map to '+'29930036 (42.54%)29783311 (42.4%)314099

16、12 (41.24%)21476601 (42.39%)19923501 (42.78%)17289330 (42.61%)Reads map to '-'29993229 (42.63%)29815598 (42.45%)31430849 (41.27%)21534570 (42.5%)19933877 (42.8%)17346445 (42.76%)Non-splice reads42357242 (60.21%)42528691 (60.55%)45227757 (59.38%)31347392 (61.87%)28062847 (60.25%)24725216 (61.

17、1%)Splice reads17566023 (24.97%)17070218 (24.3%)17613004 (23.13%)11663779 (23.02%)11794531 (25.32%)9910559 (24.26%)Reads mapped in proper pairs53795182 (76.47%)54428240 (77.49%)56181352 (73.77%)38524314 (76.04%)36101400 (77.51%)31246362 (77.25%)比對結(jié)果統(tǒng)計詳細內(nèi)容如下：(1) Total reads：測序序列經(jīng)過測序數(shù)據(jù)過濾后的數(shù)量統(tǒng)計(Clean d

18、ata)。(2) Total mapped：能定位到基因組上的測序序列的數(shù)量的統(tǒng)計；一般情況下，如果不存在污染并且參考基因組選擇合適的情況下，這部分數(shù)據(jù)的百分比大于 70%。(3) Multiple mapped：在參考序列上有多個比對位置的測序序列的數(shù)量統(tǒng)計；這部分數(shù)據(jù)的百分比一般會小于10%。(4) Uniquely mapped：在參考序列上有唯一比對位置的測序序列的數(shù)量統(tǒng)計。(5) Reads map to '+'，Reads map to '-'：測序序列比對到基因組上正鏈和負鏈的統(tǒng)計。(6) Splice reads：(2)中，分段比對到兩個外顯子上

19、的測序序列(也稱為Junction reads)的統(tǒng)計，Non-splice reads為整段比對到外顯子的將測序序列的統(tǒng)計，Splice reads的百分比取決于測序片段的長度。3.2Reads在參考基因組不同區(qū)域的分布情況對Total mapped reads的比對到基因組上的各個部分的情況進行統(tǒng)計，定位區(qū)域分為Exon(外顯子)、Intron(內(nèi)含子)和Intergenic(基因間隔區(qū)域)。正常情況下，Exon (外顯子)區(qū)域的測序序列定位的百分比含量應(yīng)該最高，定位到Intron (內(nèi)含子) 區(qū)域的測序序列可能是由于非成熟的mRNA的污染或者基因組注釋不完全導(dǎo)致的，而定位到Interge

20、nic(基因間隔區(qū)域)的測序序列可能是因為基因組注釋不完全以及背景噪音。圖3.2Reads在參考基因組不同區(qū)域的分布情況 3.3Reads在染色體上的密度分布情況對Total mapped reads的比對到基因組上的各個染色體（分正負鏈）的密度進行統(tǒng)計，如下圖所示，具體作圖的方法為用滑動窗口(window size)為1K，計算窗口內(nèi)部比對到堿基位置上的reads的中位數(shù)，并轉(zhuǎn)化成 log2 。正常情況下，整個染色體長度越長，該染色體內(nèi)部定位的reads總數(shù)會越多(Marquez et al.)。從定位到染色體上的reads數(shù)與染色體長度的關(guān)系圖中，可以更加直觀看出染色體長度和re

21、ads總數(shù)的關(guān)系。圖3.3Reads在染色體上的密度分布圖上圖：橫坐標為染色體的長度信息(以百萬堿基為單位)，縱坐標為log2(reads的密度的中位數(shù))，綠色為正鏈，紅色為負鏈 下圖：橫坐標為染色體的長度信息(單位為Mb)，縱坐標為mapped到染色體上的reads數(shù)(單位為M)3.4Reads比對結(jié)果可視化我們提供RNA-seq Reads在基因組上比對結(jié)果的bam格式文件，部分物種還提供相應(yīng)的參考基因組和注釋文件，并推薦使用IGV (Integrative Genomics Viewer) 瀏覽器對bam文件進行可視化瀏覽。IGV瀏覽器具有以下特點：(1)能在不同尺度下顯示單個或多個讀段

22、在基因組上的位置，包括讀段在各個染色體上的分布情況和在注釋的外顯子、內(nèi)含子、剪接接合區(qū)、基因間區(qū)的分布情況等；(2)能在不同尺度下顯示不同區(qū)域的讀段豐度，以反映不同區(qū)域的轉(zhuǎn)錄水平；(3)能顯示基因及其剪接異構(gòu)體的注釋信息；(4)能顯示其他注釋信息；(5)既可以從遠程服務(wù)器端下載各種注釋信息，又可以從本地加載注釋信息。IGV瀏覽器使用方法可參考我們提供的使用說明文檔(IGVQuickStart.pdf)。圖3.4IGV瀏覽器界面4可變剪切分析項目結(jié)果文件用ASprofile軟件對Cufflinks (Trapnell et al.)預(yù)測出的基因模對每個樣品的可變剪切事件分別進行分類和表達量統(tǒng)計。

23、分析流程及ASprofile中的可變剪切事件分類如下圖所示： 12類可變剪切事件定義如下:(1) TSS: Alternative 5' first exon (transcription start site) 第一個外顯子可變剪切(2) TTS: Alternative 3' last exon (transcription terminal site)最后一個外顯子可變剪切(3) SKIP: Skipped exon (SKIP_ON,SKIP_OFF pair)單外顯子跳躍(4) XSKIP: Approximate SKIP (XSKIP_ON,XSKIP_OFF p

24、air)單外顯子跳躍（模糊邊界）(5) MSKIP: Multi-exon SKIP (MSKIP_ON,MSKIP_OFF pair)多外顯子跳躍(6) XMSKIP: Approximate MSKIP (XMSKIP_ON,XMSKIP_OFF pair)多外顯子跳躍（模糊邊界）(7) IR: Intron retention (IR_ON, IR_OFF pair)單內(nèi)含子滯留(8) XIR: Approximate IR (XIR_ON, XIR_OFF pair)單內(nèi)含子滯留（模糊邊界）(9) MIR: Multi-IR (MIR_ON, MIR_OFF pair)多內(nèi)含子滯留(1

25、0) XMIR: Approximate MIR (XMIR_ON, XMIR_OFF pair)多內(nèi)含子滯留（模糊邊界）(11) AE: Alternative exon ends (5', 3', or both)可變 5'或3'端剪切(12) XAE: Approximate AE可變 5'或3'端剪切（模糊邊界） 4.1可變剪切事件分類和數(shù)量統(tǒng)計圖4.1AS分類和數(shù)量統(tǒng)計縱軸為可變剪切事件的分類縮寫，橫軸為該種事件下可變剪切的數(shù)量，不同樣品用不同子圖和顏色區(qū)分4.2可變剪切事件結(jié)構(gòu)和表達量統(tǒng)計表4.2AS結(jié)構(gòu)和表達量統(tǒng)計event_ide

26、vent_typegene_idchromevent_startevent_endevent_patternstrandfpkmref_id1000001TSSCUFF.1001343827734383303438330+1.0000000000ENSGALT000000102251000002TSSCUFF.1001345021834502533450253+3.0000000000ENSGALT000000102251000003TSSCUFF.1001345674434571653457165+2.0000000000ENSGALT000000102251000004TTSCUFF.10

27、01347480634781783474806+5.0000000000ENSGALT00000010225(1) event_id: AS事件編號(2) event_type: AS事件類型 (TSS, TTS, SKIP_ON,OFF, XSKIP_ON,OFF, MSKIP_ON,OFF, XMSKIP_ON,OFF, IR_ON ,OFF, XIR_ON,OFF, AE, XAE)(3) gene_id: cufflink組裝結(jié)果中的基因編號(4) chrom: 染色體編號(5) event_start: AS事件起始位置(6) event_end: AS事件結(jié)束位置(7) event

28、_signature: AS事件特征 (for TSS, TTS - inside boundary of alternative marginal exon; for *SKIP_ON,the coordinates of the skipped exon(s); for *SKIP_OFF, the coordinates of the enclosing introns; for *IR_ON, the end coordinates of the long, intron-containing exon; for *IR_OFF, the listing of coordinates

29、of all the exons along the path containing the retained intron; for *AE, the coordinates of the exon variant)(8) strand: 基因正負鏈信息(9) fpkm: 此AS類型該基因表達量(10) ref_id: 此基因在參考注釋文件中的編號 5新轉(zhuǎn)錄本預(yù)測項目結(jié)果文件將所有測序reads數(shù)據(jù)的基因組定位結(jié)果放到一起，用 Cufflinks 進行組裝，然后用Cuffcompare和已知的基因模型進行比較，可以:(1)發(fā)現(xiàn)新的未知基因（相對于原有基因注釋文件）；(2)發(fā)現(xiàn)已知基因新的外顯

30、子區(qū)域；(3)對已知基因的起始和終止位置進行優(yōu)化。新基因和新外顯子區(qū)域預(yù)測結(jié)果為GTF格式的注釋文件。GTF格式的詳細說明可參考(/GTF22.html)表5.1新轉(zhuǎn)錄本結(jié)構(gòu)注釋結(jié)果seqnamesourcefeaturestartendscorestrandframeattributes1novelGeneexon1853119499.+.gene_id "Novel00001" transcript_id "Novel00001.1" exon_number "1"1novelGeneex

31、on2081321813.+.gene_id "Novel00002" transcript_id "Novel00002.1" exon_number "1"1novelGeneexon2391724402.+.gene_id "Novel00003" transcript_id "Novel00003.1" exon_number "1"1novelGeneexon2518926100.+.gene_id "Novel00004" transcript

32、_id "Novel00004.1" exon_number "1"(1) seqname：染色體編號(2) source：來源標簽，這里的novelGene指新基因(3) feature：區(qū)域類型，目前我們預(yù)測外顯子區(qū)域(4) start：起始坐標(5) end：終止坐標(6) score：不必關(guān)注(7) strand：正負鏈信息(8) frame：不必關(guān)注(9) attributes：屬性，包括基因編號、轉(zhuǎn)錄本編號等信息表5.2已知基因結(jié)構(gòu)優(yōu)化Gene_idChromosomeStrandOriginal_spanAssembled_spanENSG

33、ALG000000000031+19947199199676621994719919967662ENSGALG00000000004Z-34293201342995233429320134299554ENSGALG000000000116-30928671309326163092867130932616ENSGALG0000000001322+2783575278733727835752787453(1) Gene_id：原注釋文件中基因命名編號(2) Chromosome：染色體編號(3) Strand：正負鏈信息(4) Original_span：原注釋文件中基因起始位置終止位置(5) A

34、ssembled_span：轉(zhuǎn)錄組拼接結(jié)果中基因起始位置終止位置6SNP和Indel分析項目結(jié)果文件SNP全稱Single Nucleotide Polymorphisms，是指在基因組上由單個核苷酸變異形成的遺傳標記，其數(shù)量很多，多態(tài)性豐富。從理論上來看每一個SNP 位點都可以有4 種不同的變異形式，但實際上發(fā)生的只有兩種，即轉(zhuǎn)換和顛換，二者之比為1:2。SNP在CG序列上出現(xiàn)最為頻繁，而且多是C轉(zhuǎn)換為T，原因是CG中的C常為甲基化的，自發(fā)地脫氨后即成為胸腺嘧啶。一般而言，SNP是指變異頻率大于1%的單核苷酸變異。Indel(insertion-deletion)是指相對于參考基因組，樣本中

35、發(fā)生的小片段的插入缺失，該插入缺失可能含一個或多個堿基。我們通過samtools和picard-tools等工具對比對結(jié)果進行染色體坐標排序、去掉重復(fù)的reads等處理，最后通過變異檢測軟件GATK(McKenna et al., 2010)分別進行SNP Calling和Indel Calling，并對原始結(jié)果進行過濾，得到如下表形式的分析結(jié)果。其中Indel分析結(jié)果每列的含義和SNP結(jié)果是一致的。 表6SNP分析結(jié)果#CHROMPOSREFALTHS1HS2HT1HT2HW1HW21502AG.11.00.1563CA.111100.11213AG.11.11316GA000001.000

36、0#CHROM：SNP位點所在染色體POS：SNP位點坐標REF：參考序列在該位點的基因型ALT：該位點的其它基因型other coloums：每個個體該位點的基因型(0 與REF一致；1 與ALT一致；. 缺少數(shù)據(jù)支持) 7基因表達水平分析項目結(jié)果文件一個基因表達水平的直接體現(xiàn)就是其轉(zhuǎn)錄本的豐度情況，轉(zhuǎn)錄本豐度程度越高，則基因表達水平越高。在RNA-seq分析中，我們可以通過定位到基因組區(qū)域或基因外顯子區(qū)的測序序列(reads)的計數(shù)來估計基因的表達水平。Reads計數(shù)除了與基因的真實表達水平成正比外，還與基因的長度和測序深度成正相關(guān)。為了使不同基因、不同實驗間估計的基因表達水平具有可比性，

37、人們引入了RPKM的概念，RPKM(Reads Per Kilo bases per Million reads)是每百萬reads中來自某一基因每千堿基長度的reads數(shù)目。RPKM同時考慮了測序深度和基因長度對reads計數(shù)的影響，是目前最為常用的基因表達水平估算方法 (Mortazavi et al., 2008)。結(jié)果文件分別統(tǒng)計了不同表達水平下基因的數(shù)量以及單個基因的表達水平。一般情況下，RPKM數(shù)值0.1或者1作為判斷基因是否表達的閾值，不同的文獻所采用的閾值不同。表7.1不同表達水平區(qū)間的基因數(shù)量統(tǒng)計表RPKM IntervalHS1HS2HT1HT2HW1HW20-111678

38、(44.62%)11157(42.63%)11644(44.49%)11552(44.14%)11663(44.56%)11652(44.52%)1-33416(13.05%)3829(14.63%)3497(13.36%)3622(13.84%)3359(12.83%)3503(13.39%)3-156586(25.17%)6741(25.76%)6719(25.67%)6731(25.72%)6441(24.61%)6522(24.92%)15-603436(13.13%)3421(13.07%)3278(12.53%)3277(12.52%)3612(13.80%)3442(13.15%

39、)>601055(4.03%)1023(3.91%)1033(3.95%)989(3.78%)1096(4.19%)1052(4.02%)表7.2基因表達水平統(tǒng)計表geneIDHS1HS2HT1HT2HW1HW2ENSGALG000000000030.177462748016770.5750992428404860.1820767701852560.3287502168704620.2242685548196920.283162389409366ENSGALG000000000049.444845344415057.218822649746968.93638040283769.69400

40、53828110310.40638440789559.56312902348064ENSGALG0000000001114.702355296083514.235484571846411.320256695686512.294365722600213.35758690541310.060029486986ENSGALG000000000136.748868841182057.488732937378276.668261449508977.413383587875737.536370222120317.508286101238598RNA-seq整體質(zhì)量評估項目結(jié)果文件8.1表達水平的飽和曲線檢

41、查定量飽和曲線檢查反映了基因表達水平定量對數(shù)據(jù)量的要求。表達量越高的基因，就越容易被準確定量；反之，表達量低的基因，需要較大的測序數(shù)據(jù)量才能被準確定量。表達水平的飽和曲線的具體算法描述如下：分別對10%、20%、30%90%的總體測序數(shù)據(jù)單獨進行基因定量分析，并把所有數(shù)據(jù)條件下得到的基因的表達水平作為最終的數(shù)值。用每個百分比條件下求出的單個基因的RPKM數(shù)值和最終對應(yīng)基因的表達水平數(shù)值進行比較，如果差異小于15%，則認為這個基因在這個條件下定量是準確的。圖8.1定量飽和曲線檢查分布圖橫坐標代表定位到基因組上的reads數(shù)占總reads數(shù)的百分比，縱坐標代表定量誤差在15%以內(nèi)的基因的比例8.2

42、RNA-Seq相關(guān)性檢查生物學(xué)重復(fù)是任何生物學(xué)實驗所必須的，高通量測序技術(shù)也不例外(Hansen et al.)。生物學(xué)重復(fù)主要有兩個用途：一個是證明所涉及的生物學(xué)實驗操作是可以重復(fù)的且變異不大，另一個為后續(xù)的差異基因分析所需要的。樣品間基因表達水平相關(guān)性是檢驗實驗可靠性和樣本選擇是否合理性的重要指標。相關(guān)系數(shù)越接近1，表明樣品之間表達模式的相似度越高。Encode計劃建議皮爾遜相關(guān)系數(shù)的平方(R2)大于0.92(理想的取樣和實驗條件下)。具體的項目操作中，我們要求R2至少要大于0.8，否則需要對樣品做出合適的解釋，或者重新進行實驗。此部分，我們同時計算了spearman相關(guān)系數(shù)和kendal

43、l-tau相關(guān)系數(shù)作為參考，這兩個主要是針對順序變量的相關(guān)系數(shù)，即秩相關(guān)。圖8.2RNA-Seq相關(guān)性檢查R2:pearson相關(guān)系數(shù)的平方; rho:spearman相關(guān)系數(shù); tau:kendall-tau相關(guān)系數(shù)8.3均一性分布檢查理想條件下，對于RNA-seq技術(shù)來說，測序序列(reads)之間為獨立抽樣并且reads在所有表達的轉(zhuǎn)錄本上的分布應(yīng)該呈現(xiàn)均一化分布。然而很多研究表明，很多偏好型的因素都會影響這種均一化的分布(Dohm et al., 2008)。例如，在RNA-seq建庫過程中，片段破碎和RNA反轉(zhuǎn)錄的順序不一樣會導(dǎo)致RNA-seq最終的數(shù)據(jù)呈現(xiàn)嚴重的3偏好性。其他因素還

44、包括轉(zhuǎn)錄區(qū)域的GC含量不同、隨機引物等等，并且生物體內(nèi)從5或者3的降解過程同樣會導(dǎo)致不均一性分布。圖8.3不同表達水平的轉(zhuǎn)錄本的reads密度分布圖High：高表達量轉(zhuǎn)錄本；Medium：中度表達量轉(zhuǎn)錄本；Low：低表達量轉(zhuǎn)錄本；橫坐標為距離轉(zhuǎn)錄本5端的相對位置(以百分比表示)，縱坐標為覆蓋深度的平均值9基因差異表達分析項目結(jié)果文件9.1基因表達水平對比通過所有基因的RPKM的分布圖以及盒形圖對不同實驗條件下的基因表達水平進行比較。對于同一實驗條件下的重復(fù)樣品，最終的RPKM為所有重復(fù)數(shù)據(jù)的平均值。圖9.1不同實驗條件下基因表達水平比對圖RPKM分布圖(左圖)的橫坐標為log10(RPKM), 縱坐標為基因的密度。RPKM盒形圖(右圖)的橫坐標為樣品名稱，縱坐標為log10(RPKM)，每個區(qū)域的盒形圖對五個統(tǒng)計量(至上而下分別為最大值，上四分位數(shù)，中值，下四分位數(shù)和最小值)9.2差異表達基因列表基因差異表達的輸入數(shù)據(jù)為基因表達水平分析中得到的readcount數(shù)據(jù)。對于有生物學(xué)重復(fù)的樣品，分析我們采用DESeq（Anders et al, 2010）進行分析：該分析方法基于的模型是負二項分布，第 i 個基因在第 j 個樣本中的 read count 值為Kij，則有Kij NB(i