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

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

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

3、)。illumina 測(cè)序標(biāo)識(shí)符詳細(xì)信息如下：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第四行中每個(gè)字符對(duì)應(yīng)的ASCII值減去33，即為對(duì)應(yīng)第二行堿基的測(cè)序質(zhì)量值。如果測(cè)序錯(cuò)誤率用e表示，illumina HiSeqTM2000/MiSeq的堿基質(zhì)量值用Qphred表示，則有下列關(guān)系：公式一：Qphred=-10log10(e)illumina Casava 1.8版本測(cè)

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

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

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

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

9、分析都基于clean reads。數(shù)據(jù)處理的步驟如下：(1) 去除帶接頭(adapter)的reads；(2) 去除N(N表示無(wú)法確定堿基信息)的比例大于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ù)過(guò)濾結(jié)果 2.4測(cè)序數(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ì)量情況詳細(xì)內(nèi)容如下：(1) Raw reads：統(tǒng)計(jì)原始序列數(shù)據(jù)，以四行為一個(gè)單位，統(tǒng)計(jì)每個(gè)文件的測(cè)序序列的個(gè)數(shù)。(2) Clean reads：計(jì)算方法同 Raw Reads，只是統(tǒng)計(jì)的文件為過(guò)濾后的測(cè)序數(shù)據(jù)。后續(xù)的生物信息分析都是基于Clean reads。(3) Clean bases：測(cè)序序列的個(gè)數(shù)乘以測(cè)序序列的長(zhǎng)度，并轉(zhuǎn)化為以G為單位。(4) Error rate：通過(guò)公式1計(jì)算得到。(5) Q20、Q30：分別計(jì)算 Phred 數(shù)值大于20、30的堿基占

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

13、序列的定位個(gè)數(shù)(Total Mapped Reads)及其占clean reads的百分比，其中包括多個(gè)定位的測(cè)序序列個(gè)數(shù)(Multiple Mapped Reads)及其占總體（clean reads）的百分比，以及單個(gè)定位的測(cè)序序列個(gè)數(shù)(Uniquely Mapped Reads)及其占總體（clean reads）的百分比。3.1Reads與參考基因組比對(duì)情況統(tǒng)計(jì)表3.1Reads與參考基因組比對(duì)情況一覽表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%)比對(duì)結(jié)果統(tǒng)計(jì)詳細(xì)內(nèi)容如下：(1) Total reads：測(cè)序序列經(jīng)過(guò)測(cè)序數(shù)據(jù)過(guò)濾后的數(shù)量統(tǒng)計(jì)(Clean d

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

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

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

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

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

23、分析流程及ASprofile中的可變剪切事件分類如下圖所示： 12類可變剪切事件定義如下:(1) TSS: Alternative 5' first exon (transcription start site) 第一個(gè)外顯子可變剪切(2) TTS: Alternative 3' last exon (transcription terminal site)最后一個(gè)外顯子可變剪切(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)計(jì)圖4.1AS分類和數(shù)量統(tǒng)計(jì)縱軸為可變剪切事件的分類縮寫(xiě)，橫軸為該種事件下可變剪切的數(shù)量，不同樣品用不同子圖和顏色區(qū)分4.2可變剪切事件結(jié)構(gòu)和表達(dá)量統(tǒng)計(jì)表4.2AS結(jié)構(gòu)和表達(dá)量統(tǒng)計(jì)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事件編號(hào)(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é)果中的基因編號(hào)(4) chrom: 染色體編號(hào)(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: 基因正負(fù)鏈信息(9) fpkm: 此AS類型該基因表達(dá)量(10) ref_id: 此基因在參考注釋文件中的編號(hào) 5新轉(zhuǎn)錄本預(yù)測(cè)項(xiàng)目結(jié)果文件將所有測(cè)序reads數(shù)據(jù)的基因組定位結(jié)果放到一起，用 Cufflinks 進(jìn)行組裝，然后用Cuffcompare和已知的基因模型進(jìn)行比較，可以:(1)發(fā)現(xiàn)新的未知基因（相對(duì)于原有基因注釋文件）；(2)發(fā)現(xiàn)已知基因新的外顯

30、子區(qū)域；(3)對(duì)已知基因的起始和終止位置進(jìn)行優(yōu)化。新基因和新外顯子區(qū)域預(yù)測(cè)結(jié)果為GTF格式的注釋文件。GTF格式的詳細(xì)說(shuō)明可參考(/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：染色體編號(hào)(2) source：來(lái)源標(biāo)簽，這里的novelGene指新基因(3) feature：區(qū)域類型，目前我們預(yù)測(cè)外顯子區(qū)域(4) start：起始坐標(biāo)(5) end：終止坐標(biāo)(6) score：不必關(guān)注(7) strand：正負(fù)鏈信息(8) frame：不必關(guān)注(9) attributes：屬性，包括基因編號(hào)、轉(zhuǎn)錄本編號(hào)等信息表5.2已知基因結(jié)構(gòu)優(yōu)化Gene_idChromosomeStrandOriginal_spanAssembled_spanENSG

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

34、ssembled_span：轉(zhuǎn)錄組拼接結(jié)果中基因起始位置終止位置6SNP和Indel分析項(xiàng)目結(jié)果文件SNP全稱Single Nucleotide Polymorphisms，是指在基因組上由單個(gè)核苷酸變異形成的遺傳標(biāo)記，其數(shù)量很多，多態(tài)性豐富。從理論上來(lá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(insertion-deletion)是指相對(duì)于參考基因組，樣本中

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

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

37、人們引入了RPKM的概念，RPKM(Reads Per Kilo bases per Million reads)是每百萬(wàn)reads中來(lái)自某一基因每千堿基長(zhǎng)度的reads數(shù)目。RPKM同時(shí)考慮了測(cè)序深度和基因長(zhǎng)度對(duì)reads計(jì)數(shù)的影響，是目前最為常用的基因表達(dá)水平估算方法 (Mortazavi et al., 2008)。結(jié)果文件分別統(tǒng)計(jì)了不同表達(dá)水平下基因的數(shù)量以及單個(gè)基因的表達(dá)水平。一般情況下，RPKM數(shù)值0.1或者1作為判斷基因是否表達(dá)的閾值，不同的文獻(xiàn)所采用的閾值不同。表7.1不同表達(dá)水平區(qū)間的基因數(shù)量統(tǒng)計(jì)表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基因表達(dá)水平統(tǒng)計(jì)表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ì)量評(píng)估項(xiàng)目結(jié)果文件8.1表達(dá)水平的飽和曲線檢

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

42、RNA-Seq相關(guān)性檢查生物學(xué)重復(fù)是任何生物學(xué)實(shí)驗(yàn)所必須的，高通量測(cè)序技術(shù)也不例外(Hansen et al.)。生物學(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.92(理想的取樣和實(shí)驗(yàn)條件下)。具體的項(xiàng)目操作中，我們要求R2至少要大于0.8，否則需要對(duì)樣品做出合適的解釋，或者重新進(jìn)行實(shí)驗(yàn)。此部分，我們同時(shí)計(jì)算了spearman相關(guān)系數(shù)和kendal

43、l-tau相關(guān)系數(shù)作為參考，這兩個(gè)主要是針對(duì)順序變量的相關(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均一性分布檢查理想條件下，對(duì)于RNA-seq技術(shù)來(lái)說(shuō)，測(cè)序序列(reads)之間為獨(dú)立抽樣并且reads在所有表達(dá)的轉(zhuǎn)錄本上的分布應(yīng)該呈現(xiàn)均一化分布。然而很多研究表明，很多偏好型的因素都會(huì)影響這種均一化的分布(Dohm et al., 2008)。例如，在RNA-seq建庫(kù)過(guò)程中，片段破碎和RNA反轉(zhuǎn)錄的順序不一樣會(huì)導(dǎo)致RNA-seq最終的數(shù)據(jù)呈現(xiàn)嚴(yán)重的3偏好性。其他因素還

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