有参考基因组的转录组生物信息分析模板.docx

1、有参考基因组的转录组生物信息分析模板一、生物信息分析流程获得原始测序序列(Sequenced Reads)后，在有相关物种参考序列或参考基因组的情况下，通过如下流程进行生物信息分析：二、项目结果说明1原始序列数据项目结果文件高通量测序(如illumina HiSeqTM2000/MiSeq等测序平台)测序得到的原始图像数据文件经碱基识别(Base Calling)分析转化为原始测序序列(Sequenced Reads)，我们称之为Raw Data或Raw Reads，结果以FASTQ(简称为fq)文件格式存储，其中包含测序序列(reads)的序列信息以及其对应的测序质量信息。测序样品真实数据随

2、机截取展示FASTQ格式文件中每个read由四行描述，如下：EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG GCTCTTTGCCCTTCTCGTCGAAAATTGTCTCCTCATTCGAAACTTCTCTGT + CFFFDEHHHHFIJJJFHGIIIEHIIJBHHHIJJEGIIJJIGHIGHCCF 其中第一行以“”开头，随后为illumina 测序标识符(Sequence Identifiers)和描述文字(选择性部分)；第二行是碱基序列；第三行以“+”开头，随后为illumina 测序标识符(选择性部分)；第四行是对应序

3、列的测序质量(Cock et al.)。illumina 测序标识符详细信息如下：EAS139Unique instrument name136Run IDFC706VJFlowcell ID2Flowcell lane2104Tile number within the flowcell lane15343x-coordinate of the cluster within the tile197393y-coordinate of the cluster within the tile1Member of a pair, 1 or 2 (paired-end or mate-pair re

4、ads 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第四行中每个字符对应的ASCII值减去33，即为对应第二行碱基的测序质量值。如果测序错误率用e表示，illumina HiSeqTM2000/MiSeq的碱基质量值用Qphred表示，则有下列关系：公式一：Qphred=-10log10(e)illumina Casava 1.8版本测序错

5、误率与测序质量值简明对应关系如下：测序错误率测序质量值对应字符5%13.1%2050.1%30?0.01%40I2测序数据质量评估项目结果文件2.1测序错误率分布检查每个碱基测序错误率是通过测序Phred数值(Phred score, Qphred)通过公式1转化得到，而Phred 数值是在碱基识别(Base Calling)过程中通过一种预测碱基判别发生错误概率模型计算得到的，对应关系如下表所显示：illumina Casava 1.8版本碱基识别与Phred分值之间的简明对应关系Phred分值不正确的碱基识别碱基正确识别率Q-sorce101/1090%Q10201/10099%Q2030

6、1/100099.9%Q30401/1000099.99%Q40测序错误率与碱基质量有关，受测序仪本身、测序试剂、样品等多个因素共同影响。对于RNA-seq技术，测序错误率分布具有两个特点：(1)测序错误率会随着测序序列(Sequenced Reads)的长度的增加而升高，这是由于测序过程中化学试剂的消耗而导致的，并且为illumina高通量测序平台都具有的特征(Erlich and Mitra, 2008; Jiang et al.)。(2)前6个碱基的位置也会发生较高的测序错误率，而这个长度也正好等于在RNA-seq建库过程中反转录所需要的随机引物的长度。所以推测前6个碱基测序错误率较高的

7、原因为随机引物和RNA模版的不完全结合(Jiang et al.)。测序错误率分布检查用于检测在测序长度范围内，有无异常的碱基位置存在高错误率，比如中间位置的碱基测序错误率显著高于其他位置。一般情况下，每个碱基位置的测序错误率都应该低于0.5%。图2.1测序错误率分布图横坐标为reads的碱基位置，纵坐标为单碱基错误率2.2GC含量分布检查GC含量分布检查用于检测有无AT、GC 分离现象，而这种现象可能是测序或者建库所带来的，并且会影响后续的定量分析。在illumina测序平台的转录组测序中，反转录成cDNA时所用的6bp 的随机引物会引起前几个位置的核苷酸组成存在一定的偏好性。而这种偏好性与

8、测序的物种和实验室环境无关，但会影响转录组测序的均一化程度(Hansen et al.)。除此之外，理论上G和C碱基及A和T碱基含量每个测序循环上应分别相等，且整个测序过程稳定不变，呈水平线。对于DGE测序来说，由于随机引物扩增偏差等原因，常常会导致在测序得到的每个read前6-7个碱基有较大的波动，这种波动属于正常情况。图2.2GC含量分布图横坐标为reads的碱基位置，纵坐标为单碱基所占的比例；不同颜色代表不同的碱基类型2.3测序数据过滤测序得到的原始测序序列，里面含有带接头的、低质量的reads，为了保证信息分析质量，必须对raw reads进行过滤，得到clean reads，后续分析

9、都基于clean reads。数据处理的步骤如下：(1) 去除带接头(adapter)的reads；(2) 去除N(N表示无法确定碱基信息)的比例大于10%的reads；(3) 去除低质量reads。RNA-seq 的接头(Adapter, Oligonucleotide sequences for TruSeqTM RNA and DNA Sample Prep Kits) 信息：RNA 5 Adapter (RA5), part # 15013205：5-AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT-3RNA 3 Ad

10、apter (RA3), part # 15013207：5-GATCGGAAGAGCACACGTCTGAACTCCAGTCAC(6位index)ATCTCGTATGCCGTCTTCTGCTTG-3图2.3原始数据过滤结果2.4测序数据质量情况汇总表2.4数据产出质量情况一览表Sample nameRaw readsClean readsclean basesError rate(%)Q20(%)Q30(%)GC content(%)HS1_136579608351752053.52G0.0397.8892.8849.39HS1_236579608351752053.52G0.0396.509

11、0.3849.59HS2_136547734351194633.51G0.0397.8592.8149.53数据质量情况详细内容如下：(1) Raw reads：统计原始序列数据，以四行为一个单位，统计每个文件的测序序列的个数。(2) Clean reads：计算方法同 Raw Reads，只是统计的文件为过滤后的测序数据。后续的生物信息分析都是基于Clean reads。(3) Clean bases：测序序列的个数乘以测序序列的长度，并转化为以G为单位。(4) Error rate：通过公式1计算得到。(5) Q20、Q30：分别计算 Phred 数值大于20、30的碱基占总体碱基的百分比

12、。(6) GC content：计算碱基G和C的数量总和占总的碱基数量的百分比。3参考序列比对分析项目结果文件测序序列定位算法：根据不同的基因组的特征，我们选取相对合适的软件(动植物用TopHat(Trapnell et al., 2009)、真菌或者基因密度较高的物种用Bowtie)，合适的参数设置(如最大的内含子长度，会根据已知的该物种的基因模型来进行统计分析)，将过滤后的测序序列进行基因组定位分析。下图为TopHat的算法示意图： Tophat的算法主要分为两个部分：(1) 将测序序列整段比对到外显子上。(2) 将测序序列分段比对到两个外显子上。我们统计了实验所产生的测序序列的定位个数(

13、Total Mapped Reads)及其占clean reads的百分比，其中包括多个定位的测序序列个数(Multiple Mapped Reads)及其占总体（clean reads）的百分比，以及单个定位的测序序列个数(Uniquely Mapped Reads)及其占总体（clean reads）的百分比。3.1Reads与参考基因组比对情况统计表3.1Reads与参考基因组比对情况一览表Sample nameHS1HS2HT1HT2HW1HW2Total reads703504107023892676161678506660844657366240543118Total mapped

14、60529821 (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 (85.58%)

15、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%)31409912 (41.24%)2147660

16、1 (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.1%)Splice reads17566023 (24.

17、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%)比对结果统计详细内容如下：(1) Total reads：测序序列经过测序数据过滤后的数量统计(Clean data)。(2) Total mapped：能定位到基因

18、组上的测序序列的数量的统计；一般情况下，如果不存在污染并且参考基因组选择合适的情况下，这部分数据的百分比大于 70%。(3) Multiple mapped：在参考序列上有多个比对位置的测序序列的数量统计；这部分数据的百分比一般会小于10%。(4) Uniquely mapped：在参考序列上有唯一比对位置的测序序列的数量统计。(5) Reads map to +，Reads map to -：测序序列比对到基因组上正链和负链的统计。(6) Splice reads：(2)中，分段比对到两个外显子上的测序序列(也称为Junction reads)的统计，Non-splice reads为整段比

19、对到外显子的将测序序列的统计，Splice reads的百分比取决于测序片段的长度。3.2Reads在参考基因组不同区域的分布情况对Total mapped reads的比对到基因组上的各个部分的情况进行统计，定位区域分为Exon(外显子)、Intron(内含子)和Intergenic(基因间隔区域)。正常情况下，Exon (外显子)区域的测序序列定位的百分比含量应该最高，定位到Intron (内含子) 区域的测序序列可能是由于非成熟的mRNA的污染或者基因组注释不完全导致的，而定位到Intergenic(基因间隔区域)的测序序列可能是因为基因组注释不完全以及背景噪音。图3.2Reads在参考

20、基因组不同区域的分布情况3.3Reads在染色体上的密度分布情况对Total mapped reads的比对到基因组上的各个染色体（分正负链）的密度进行统计，如下图所示，具体作图的方法为用滑动窗口(window size)为1K，计算窗口内部比对到碱基位置上的reads的中位数，并转化成 log2 。正常情况下，整个染色体长度越长，该染色体内部定位的reads总数会越多(Marquez et al.)。从定位到染色体上的reads数与染色体长度的关系图中，可以更加直观看出染色体长度和reads总数的关系。图3.3Reads在染色体上的密度分布图上图：横坐标为染色体的长度信息(以百万碱基为单位)

21、，纵坐标为log2(reads的密度的中位数)，绿色为正链，红色为负链 下图：横坐标为染色体的长度信息(单位为Mb)，纵坐标为mapped到染色体上的reads数(单位为M)3.4Reads比对结果可视化我们提供RNA-seq Reads在基因组上比对结果的bam格式文件，部分物种还提供相应的参考基因组和注释文件，并推荐使用IGV (Integrative Genomics Viewer) 浏览器对bam文件进行可视化浏览。IGV浏览器具有以下特点：(1)能在不同尺度下显示单个或多个读段在基因组上的位置，包括读段在各个染色体上的分布情况和在注释的外显子、内含子、剪接接合区、基因间区的分布情况等

22、；(2)能在不同尺度下显示不同区域的读段丰度，以反映不同区域的转录水平；(3)能显示基因及其剪接异构体的注释信息；(4)能显示其他注释信息；(5)既可以从远程服务器端下载各种注释信息，又可以从本地加载注释信息。IGV浏览器使用方法可参考我们提供的使用说明文档(IGVQuickStart.pdf)。图3.4IGV浏览器界面4可变剪切分析项目结果文件用ASprofile软件对Cufflinks (Trapnell et al.)预测出的基因模对每个样品的可变剪切事件分别进行分类和表达量统计。分析流程及ASprofile中的可变剪切事件分类如下图所示： 12类可变剪切事件定义如下:(1) TSS:

23、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 pair)单外显子跳跃（模糊边界）(5) MSKIP: Multi-exon SKIP (MSKIP_ON,MSKIP_OFF p

24、air)多外显子跳跃(6) XMSKIP: Approximate MSKIP (XMSKIP_ON,XMSKIP_OFF pair)多外显子跳跃（模糊边界）(7) IR: Intron retention (IR_ON, IR_OFF pair)单内含子滞留(8) XIR: Approximate IR (XIR_ON, XIR_OFF pair)单内含子滞留（模糊边界）(9) MIR: Multi-IR (MIR_ON, MIR_OFF pair)多内含子滞留(10) XMIR: Approximate MIR (XMIR_ON, XMIR_OFF pair)多内含子滞留（模糊边界）(11

25、) AE: Alternative exon ends (5, 3, or both)可变 5或3端剪切(12) XAE: Approximate AE可变 5或3端剪切（模糊边界） 4.1可变剪切事件分类和数量统计图4.1AS分类和数量统计纵轴为可变剪切事件的分类缩写，横轴为该种事件下可变剪切的数量，不同样品用不同子图和颜色区分4.2可变剪切事件结构和表达量统计表4.2AS结构和表达量统计event_idevent_typegene_idchromevent_startevent_endevent_patternstrandfpkmref_id1000001TSSCUFF.100134382

26、7734383303438330+1.0000000000ENSGALT000000102251000002TSSCUFF.1001345021834502533450253+3.0000000000ENSGALT000000102251000003TSSCUFF.1001345674434571653457165+2.0000000000ENSGALT000000102251000004TTSCUFF.1001347480634781783474806+5.0000000000ENSGALT00000010225(1) event_id: AS事件编号(2) event_type: AS事件

27、类型 (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组装结果中的基因编号(4) chrom: 染色体编号(5) event_start: AS事件起始位置(6) event_end: AS事件结束位置(7) event_signature: AS事件特征 (for TSS, TTS - inside boundary of alternative marginal exon; for *SKIP_ON,

28、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 of all the exons along the path containing the retained intron; for *AE, the coordinates of th

29、e exon variant)(8) strand: 基因正负链信息(9) fpkm: 此AS类型该基因表达量(10) ref_id: 此基因在参考注释文件中的编号 5新转录本预测项目结果文件将所有测序reads数据的基因组定位结果放到一起，用 Cufflinks 进行组装，然后用Cuffcompare和已知的基因模型进行比较，可以:(1)发现新的未知基因（相对于原有基因注释文件）；(2)发现已知基因新的外显子区域；(3)对已知基因的起始和终止位置进行优化。新基因和新外显子区域预测结果为GTF格式的注释文件。GTF格式的详细说明可参考(http:/mblab.wustl.edu/GTF22.h

30、tml)表5.1新转录本结构注释结果seqnamesourcefeaturestartendscorestrandframeattributes1novelGeneexon1853119499.+.gene_id Novel00001; transcript_id Novel00001.1; exon_number 1;1novelGeneexon2081321813.+.gene_id Novel00002; transcript_id Novel00002.1; exon_number 1;1novelGeneexon2391724402.+.gene_id Novel00003; transcript_id Novel000