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

一、生物信息分析流程获得原始测序序列(Sequenced Reads)后，在有相关物种参考序列或参考基因组的情况下，通过如下流程进行生物信息分析：二、项目结果说明1　原始序列数据项目结果文件高通量测序(如illumina HiSeqTM2000/MiSeq等测序平台)测序得到的原始图像数据文件经碱基识别(Base Calling)分析转化为原始测序序列(Sequenced Reads)，我们称之为Raw Data或Raw Reads，结果以FASTQ(简称为fq)文件格式存储，其中包含测序序列(reads)的序列信息以及其对应的测序质量信息测序样品真实数据随机截取展示FASTQ格式文件中每个read由四行描述，如下：@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG GCTCTTTGCCCTTCTCGTCGAAAATTGTCTCCTCATTCGAAACTTCTCTGT + @@CFFFDEHHHHFIJJJ@FHGIIIEHIIJBHHHIJJEGIIJJIGHIGHCCF 其中第一行以“@”开头，随后为illumina 测序标识符(Sequence Identifiers)和描述文字(选择性部分)；第二行是碱基序列；第三行以“+”开头，随后为illumina 测序标识符(选择性部分)；第四行是对应序列的测序质量(Cock et al.)。

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 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第四行中每个字符对应的ASCII值减去33，即为对应第二行碱基的测序质量值如果测序错误率用e表示，illumina HiSeqTM2000/MiSeq的碱基质量值用Qphred表示，则有下列关系：公式一：　Qphred　=　-10log10(e)illumina Casava 1.8版本测序错误率与测序质量值简明对应关系如下：测序错误率测序质量值对应字符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%Q20301/100099.9%Q30401/1000099.99%Q40测序错误率与碱基质量有关，受测序仪本身、测序试剂、样品等多个因素共同影响。

对于RNA-seq技术，测序错误率分布具有两个特点：(1)测序错误率会随着测序序列(Sequenced Reads)的长度的增加而升高，这是由于测序过程中化学试剂的消耗而导致的，并且为illumina高通量测序平台都具有的特征(Erlich and Mitra, 2008; Jiang et al.)2)前6个碱基的位置也会发生较高的测序错误率，而这个长度也正好等于在RNA-seq建库过程中反转录所需要的随机引物的长度所以推测前6个碱基测序错误率较高的原因为随机引物和RNA模版的不完全结合(Jiang et al.)测序错误率分布检查用于检测在测序长度范围内，有无异常的碱基位置存在高错误率，比如中间位置的碱基测序错误率显著高于其他位置一般情况下，每个碱基位置的测序错误率都应该低于0.5%图2.1　测序错误率分布图横坐标为reads的碱基位置，纵坐标为单碱基错误率2.2　GC含量分布检查GC含量分布检查用于检测有无AT、GC 分离现象，而这种现象可能是测序或者建库所带来的，并且会影响后续的定量分析在illumina测序平台的转录组测序中，反转录成cDNA时所用的6bp 的随机引物会引起前几个位置的核苷酸组成存在一定的偏好性。

而这种偏好性与测序的物种和实验室环境无关，但会影响转录组测序的均一化程度(Hansen et al.)除此之外，理论上G和C碱基及A和T碱基含量每个测序循环上应分别相等，且整个测序过程稳定不变，呈水平线对于DGE测序来说，由于随机引物扩增偏差等原因，常常会导致在测序得到的每个read前6-7个碱基有较大的波动，这种波动属于正常情况图2.2　GC含量分布图横坐标为reads的碱基位置，纵坐标为单碱基所占的比例；不同颜色代表不同的碱基类型2.3　测序数据过滤测序得到的原始测序序列，里面含有带接头的、低质量的reads，为了保证信息分析质量，必须对raw reads进行过滤，得到clean reads，后续分析都基于clean reads数据处理的步骤如下：(1) 去除带接头(adapter)的reads；(2) 去除N(N表示无法确定碱基信息)的比例大于10%的reads；(3) 去除低质量readsRNA-seq 的接头(Adapter, Oligonucleotide sequences for TruSeqTM RNA and DNA Sample Prep Kits) 信息：RNA 5’ Adapter (RA5), part # 15013205：5’-AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT-3’RNA 3’ Adapter (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.5090.3849.59HS2_136547734351194633.51G0.0397.8592.8149.53数据质量情况详细内容如下：(1) Raw reads：统计原始序列数据，以四行为一个单位，统计每个文件的测序序列的个数。

2) Clean reads：计算方法同 Raw Reads，只是统计的文件为过滤后的测序数据后续的生物信息分析都是基于Clean reads3) Clean bases：测序序列的个数乘以测序序列的长度，并转化为以G为单位4) Error rate：通过公式1计算得到5) Q20、Q30：分别计算 Phred 数值大于20、30的碱基占总体碱基的百分比6) GC content：计算碱基G和C的数量总和占总的碱基数量的百分比3　参考序列比对分析项目结果文件测序序列定位算法：根据不同的基因组的特征，我们选取相对合适的软件(动植物用TopHat(Trapnell et al., 2009)、真菌或者基因密度较高的物种用Bowtie)，合适的参数设置(如最大的内含子长度，会根据已知的该物种的基因模型来进行统计分析)，将过滤后的测序序列进行基因组定位分析下图为TopHat的算法示意图： Tophat的算法主要分为两个部分：(1) 将测序序列整段比对到外显子上2) 将测序序列分段比对到两个外显子上我们统计了实验所产生的测序序列的定位个数(Total Mapped Reads)及其占clean reads的百分比，其中包括多个定位的测序序列个数(Multiple Mapped Reads)及其占总体（clean reads）的百分比，以及单个定位的测序序列个数(Uniquely Mapped Reads)及其占总体（clean reads）的百分比。

3.1　Reads与参考基因组比对情况统计表3.1　Reads与参考基因组比对情况一览表Sample nameHS1HS2HT1HT2HW1HW2Total reads703504107023892676161678506660844657366240543118Total 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 (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%)31409912 (41.24%)21476601 (42.39%)19923501 (42。