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

一、生物信息剖析流程获取原始测序序列(SequencedReads)后，在有有关物种参照序列或参照基因组的状况下，经过以下贱程进行生物信息剖析：二、项目结果说明原始序列数据高通量测序(如illuminaHiSeqTM2000/MiSeq等测序平台)测序获取的原始图像数据文件经碱基辨别(BaseCalling)剖析转变为原始测序序列(SequencedReads)，我们称之为RawData或RawReads，结果以FASTQ(简称为fq)文件格式储存，此中包含测序序列(reads)的序列信息以及其对应的测序质量信息。FASTQ格式文件中每个read由四行描绘，以下：@EAS139:136:FC706VJ:2:2104:15343:1973931:Y:18:ATCACGGCTCTTTGCCCTTCTCGTCGAAAATTGTCTCCTCATTCGAAACTTCTCTGT+@@CFFFDEHHHHFIJJJ@FHGIIIEHIIJBHHHIJJEGIIJJIGHIGHCCF此中第一行以“@”开头，随后为illumina测序表记符(SequenceIdentifiers)和描绘文字(选择性部分)；第二行是碱基序列；第三行以“+”开头，随后为illumina测序表记符(选择性部分)；第四行是对应序列的测序质量(Cocketal.)

。illumina

测序表记符详尽信息以下：EAS139136FC706VJ22104153431973931Y18ATCACG

UniqueinstrumentnameRunIDFlowcellIDFlowcelllaneTilenumberwithintheflowcelllane'x'-coordinateoftheclusterwithinthetile'y'-coordinateoftheclusterwithinthetileMemberofapair,1or2(paired-endormate-pairreadsonly)Yifthereadfailsfilter(readisbad),Notherwise0whennoneofthecontrolbitsareon,otherwiseitisanevennumberIndexsequence第四行中每个字符对应的ASCII值减去量值。假如测序错误率用e表示，illuminaQphred表示，则有以下关系：

33，即为对应第二行碱基的测序质HiSeqTM2000/MiSeq的碱基质量值用公式一：

Qphred

=-10log

10(e)illuminaCasava版本测序错误率与测序质量值简洁对应关系以下：测序错误率测序质量值对应字符5%13.1%205%30?%40I测序数据质量评估测序错误率散布检查每个碱基测序错误率是经过测序Phred数值(Phredscore,Qphred)经过公式1转变获取，而Phred数值是在碱基辨别(BaseCalling)过程中经过一种展望碱基鉴别发生错误概率模型计算获取的，对应关系以下表所显示：illuminaCasava版本碱基辨别与Phred分值之间的简洁对应关系Phred分值不正确的碱基辨别碱基正确辨别率Q-sorce101/1090%Q10201/10099%Q20301/1000%Q30401/10000%Q40测序错误率与碱基质量有关，受测序仪自己、测序试剂、样品等多个要素共同影响。关于RNA-seq技术，测序错误率散布拥有两个特色：测序错误率会跟着测序序列(SequencedReads)的长度的增添而高升，这是因为测序过程中化学试剂的耗费而致使的，并且为illumina高通量测序平台都拥有的特色(ErlichandMitra,2008;Jiangetal.)。前6个碱基的地点也会发生较高的测序错误率，而这个长度也正好等于在RNA-seq建库过程中反转录所需要的随机引物的长度。因此推测前6个碱基测序错误率较高的原由为随机引物和RNA模版的不完整联合(Jiangetal.)。测序错误率散布检查用于检测在测序长度范围内，有无异样的碱基地点存在高错误率，比方中间地点的碱基测序错误率明显高于其余地点。一般状况下，每个碱基地点的测序错误率都应当低于%。图测序错误率散布图横坐标为reads的碱基地点，纵坐标为单碱基错误率GC含量散布检查GC含量散布检查用于检测有无AT、GC分别现象，而这类现象可能是测序或许建库所带来的，并且会影响后续的定量剖析。在illumina测序平台的转录组测序中，反转录成cDNA时所用的6bp的随机引物会惹起前几个地点的核苷酸构成存在必定的偏好性。而这类偏好性与测序的物种和实验室环境没关，但会影响转录组测序的均一化程度(Hansenetal.)。除此以外，理论上G和C碱基及A和T碱基含量每个测序循环上应分别相等，且整个测序过程稳固不变，呈水平线。关于DGE测序来说，因为随机引物扩增偏差等原由，经常会致使在测序获取的每个read前6-7个碱基有较大的颠簸，这类颠簸属于正常状况。图GC含量散布图横坐标为reads的碱基地点，纵坐标为单碱基所占的比率；不同颜色代表不同的碱基种类测序数据过滤测序获取的原始测序序列，里面含有带接头的、低质量的reads，为了保证信息剖析质量，一定对rawreads进行过滤，获取cleanreads，后续剖析都基于cleanreads。数据办理的步骤以下：去除带接头(adapter)的reads；去除N(N表示没法确立碱基信息)的比率大于10%的reads；去除低质量reads。RNA-seq的接头(Adapter,OligonucleotidesequencesforTruSeqandDNASamplePrepKits)信息：

TMRNARNA5’Adapter(RA5),part#

：5’-AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT-3’RNA3’Adapter(RA3),part#：5’-GATCGGAAGAGCACACGTCTGAACTCCAGTCAC(6位index)ATCTCGTATGCCGTCTTCTGCTTG3’-图原始数据过滤结果测序数据质量状况汇总表数据产出质量状况一览表SamplenameRawreadsCleanreadscleanbasesErrorrate(%)Q20(%)Q30(%)GCcontent(%)HS1_1HS1_2HS2_1数据质量状况详尽内容以下：Rawreads：统计原始序列数据，以四行为一个单位，统计每个文件的测序序列的个数。Cleanreads：计算方法同RawReads，不过统计的文件为过滤后的测序数据。后续的生物信息剖析都是鉴于Cleanreads。(3)Cleanbases：测序序列的个数乘以测序序列的长度，并转变为以G为单位。Errorrate：经过公式1计算获取。Q20、Q30：分别计算Phred数值大于20、30的碱基占整体碱基的百分比。GCcontent：计算碱基G和C的数目总和占总的碱基数目的百分比。参照序列比对剖析测序序列定位算法：依据不同的基因组的特色，我们选用相对适合的软件(动植物用TopHat(Trapnelletal.,2009)、真菌或许基因密度较高的物种用Bowtie)，适合的参数设置(如最大的内含子长度，会依据已知的该物种的基因模型来进行统计剖析)，将过滤后的测序序列进行基因组定位剖析。以下图为TopHat的算法表示图：Tophat的算法主要分为两个部分：将测序序列整段比对到外显子上。将测序序列分段比对到两个外显子上。我们统计了实验所产生的测序序列的定位个数(TotalMappedReads)及其占cleanreads的百分比，此中包含多个定位的测序序列个数(MultipleMappedReads)及其占整体（cleanreads）的百分比，以及单个定位的测序序列个数(UniquelyMappedReads)及其占整体（cleanreads）的百分比。Reads与参照基因组比对状况统计表Reads与参照基因组比对状况一览表SamplenameHS1HS2HT1HT2HW1HW2TotalreadsTotalmapped%)%)%)%)%)%)Multiplemapped606556%)633575%)714678%)450156%)389470%)335509%)Uniquelymapped%)%)%)%)%)%)Read-1%)%)%)%)(43%)%)Read-2%)%)%)%)%)%)Readsmapto'+'%)%)%)%)%)%)Readsmapto'-'%)%)%)%)%)%)Non-splicereads%)%)%)%)%)%)Splicereads%)%)%)%)%)9910559%)Readsmappedinproperpairs%)%)%)%)%)%)比对结果统计详尽内容以下：Totalreads：测序序列经过测序数据过滤后的数目统计(Cleandata)。Totalmapped：能定位到基因组上的测序序列的数目的统计；一般状况下，假如不存在污染并且参照基因组选择适合的状况下，这部分数据的百分比大于70%。Multiplemapped：在参照序列上有多个比对地点的测序序列的数目统计；这部分数据的百分比一般会小于10%。Uniquelymapped：在参照序列上有独一比对地点的测序序列的数目统计。(5)Readsmapto'+'，Readsmapto'-'：测序序列比对到基因组上正链和负链的统计。(6)Splicereads：(2)中，分段比对到两个外显子上的测序序列(也称为Junctionreads)的统计，Non-splicereads为整段比对到外显子的将测序序列的统计，Splicereads的百分比取决于测序片段的长度。Reads在参照基因组不同地区的散布状况对Totalmappedreads的比对到基因组上的各个部分的状况进行统计，定位地划分为Exon(外显子)、Intron(内含子)和Intergenic(基因间隔地区)。正常状况下，Exon(外显子)地区的测序序列定位的百分比含量应当最高，定位到Intron(内含子)地区的测序序列可能是因为非成熟的mRNA的污染或许基因组说明不完整致使的，而定位到Intergenic(基因间隔地区)的测序序列可能是因为基因组说明不完整以及背景噪音。图Reads在参照基因组不同地区的散布状况Reads在染色体上的密度散布状况对Totalmappedreads的比对到基因组上的各个染色体（分正负链）的密度进行统计，以以下图所示，详细作图的方法为用滑动窗口(windowsize)为1K，计算窗口内部比对到碱基地点上的reads的中位数，并转变成log2。正常状况下，整个染色体长度越长，该染色体内部定位的reads总数会越多(Marquezetal.)。从定位到染色体上的reads数与染色体长度的关系图中，能够更为直观看出染色体长度和reads总数的关系。图Reads在染色体上的密度散布图上图：横坐标为染色体的长度信息(以百万碱基为单位)，纵坐标为log2(reads的密度的中位数)，绿色为正链，红色为负链以下图：横坐标为染色体的长度信息(单位为Mb)，纵坐标为mapped到染色体上的reads数(单位为M)Reads比对结果可视化我们供给RNA-seqReads在基因组上比对结果的bam格式文件，部分物种还供给相应的参照基因组和说明文件，并介绍使用IGV(IntegrativeGenomicsViewer)阅读器对bam文件进行可视化阅读。IGV阅读器拥有以下特色：(1)能在不同尺度下显示单个或多个读段在基因组上的地点，包含读段在各个染色体上的散布状况和在说明的外显子、内含子、剪接接合区、基因间区的散布状况等；能在不同尺度下显示不同地区的读段丰度，以反应不同地区的转录水平；(3)能显示基因及其剪接异构体的说明信息；(4)能显示其余说明信息；(5)既能够从远程服务器端下载各样说明信息，又能够从当地加载说明信息。IGV阅读器使用方法可参照我们供给的使用说明文档(。图IGV阅读器界面可变剪切剖析用ASprofile软件对Cufflinks(Trapnelletal.)展望出的基因模对每个样品的可变剪切事件分别进行分类和表达量统计。剖析流程及ASprofile中的可变剪切事件分类以以下图所示：12类可变剪切事件定义以下:(1)TSS:Alternative5'firstexon(transcriptionstartsite)第一个外显子可变剪切(2)TTS:Alternative3'lastexon(transcriptionterminalsite)最后一个外显子可变剪切(3)SKIP:Skippedexon(SKIP_ON,SKIP_OFFpair)单外显子跳跃(4)XSKIP:ApproximateSKIP(XSKIP_ON,XSKIP_OFFpair)单外显子跳跃（模糊界限）(5)MSKIP:Multi-exonSKIP(MSKIP_ON,MSKIP_OFFpair)多外显子跳跃(6)XMSKIP:ApproximateMSKIP(XMSKIP_ON,XMSKIP_OFFpair)多外显子跳跃（模糊界限）(7)IR:Intronretention(IR_ON,IR_OFFpair)单内含子滞留(8)XIR:ApproximateIR(XIR_ON,XIR_OFFpair)单内含子滞留（模糊界限）(9)MIR:Multi-IR(MIR_ON,MIR_OFFpair)多内含子滞留(10)XMIR:ApproximateMIR(XMIR_ON,XMIR_OFFpair)多内含子滞留（模糊界限）(11)AE:Alternativeexonends(5',3',orboth)可变5'或3'端剪切XAE:ApproximateAE可变5'或3'端剪切（模糊界限）可变剪切事件分类和数目统计图AS分类和数目统计纵轴为可变剪切事件的分类缩写，横轴为该种事件下可变剪切的数目，不相同品用不同子图和颜色划分可变剪切事件构造和表达量统计表AS构造和表达量统计event_ievent_typgene_ichroevent_starevent_enevent_patterstranfpkdedmtdndref_idm1000001TSS1343827734383303438330+ENSGALT000000102251000002TSS1345021834502533450253+ENSGALT00000010225ENSGALT000000102+25ENSGALT000000102+25(1)event_id:AS事件编号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组装结果中的基因编号(4)chrom:染色体编号(5)event_start:AS事件开端地点event_end:AS事件结束地点event_signature:AS事件特色(forTSS,TTS-insideboundaryofalternativemarginalexon;for*SKIP_ON,thecoordinatesoftheskippedexon(s);for*SKIP_OFF,thecoordinatesoftheenclosingintrons;for*IR_ON,theendcoordinatesofthelong,intron-containingexon;for*IR_OFF,thelistingofcoordinatesofalltheexonsalongthepathcontainingtheretainedintron;for*AE,thecoordinatesoftheexonvariant)(8)strand:基因正负链信息fpkm:此AS种类该基因表达量ref_id:此基因在参照说明文件中的编号新转录本展望将全部测序reads数据的基因组定位结果放到一同，用Cufflinks进行组装，而后用Cuffcompare和已知的基因模型进行比较，能够:(1)发现新的未知基因（有关于原有基因说明文件）；(2)发现已知基因新的外显子地区；(3)对已知基因的开端和停止地点进行优化。新基因和新外显子地区展望结果为GTF格式的说明文件。GTF格式的详尽说明可参照(表新转录本构造说明结果seqnamesourcefeaturestartendscorestrandframeattributes1novelGeneexon1853119499.gene_id"Novel00001";transcript_id+."";exon_number"1";1novelGeneexon2081321813.gene_id"Novel00002";transcript_id+."";exon_number"1";1novelGeneexon2391724402.gene_id"Novel00003";transcript_id+."";exon_number"1";1novelGeneexon2518926100.gene_id"Novel00004";transcript_id+."";exon_number"1";seqname：染色体编号source：根源标签，这里的novelGene指新基因feature：地区种类，当前我们展望外显子地区start：开端坐标end：停止坐标score：不用关注strand：正负链信息frame：不用关注attributes：属性，包含基因编号、转录本编号等信息表已知基因构造优化Gene_idChromosomeStrandOriginal_spanAssembled_spanENSGALG000000000031+~~ENSGALG00000000004Z-~~ENSGALG000000000116-~~ENSGALG0000000001322+2783575~27873372783575~2787453Gene_id：原说明文件中基因命名编号Chromosome：染色体编号Strand：正负链信息Original_span：原说明文件中基因开端地点~停止地点Assembled_span：转录组拼接结果中基因开端地点~停止地点6SNP和Indel剖析SNP全称SingleNucleotidePolymorphisms，是指在基因组上由单个核苷酸变异形成的遗传标志，其数目好多，多态性丰富。从理论上来看每一个SNP位点都能够有4种不同的变异形式，但实质上发生的只有两种，即变换和颠换，两者之比为1:2。SNP在CG序列上出现最为屡次，并且多是C变换为T，原由是CG中的C常为甲基化的，自觉地脱氨后即成为胸腺嘧啶。一般而言，SNP是指变异频次大于1%的单核苷酸变异。Indel(insertion-deletion)是指有关于参照基因组，样本中发生的小片段的插入缺失，该插入缺失可能含一个或多个碱基。我们经过samtools和picard-tools等工具对照对结果进行染色体坐标排序、去掉重复的reads等办理，最后经过变异检测软件GATK(McKennaetal.,2010)分别进行SNPCalling和IndelCalling，并对原始结果进行过滤，获取以下表形式的剖析结果。此中Indel剖析结果每列的含义和SNP结果是一致的。表

6

SNP剖析结果#CHROM

POS

REF

ALT

HS1

HS2

HT1

HT2

HW1

HW21

502A

G

..

..

11..

00..1563CA....111100..11213AG....11......11316GA000001..0000#CHROM：SNP位点所在染色体POS：SNP位点坐标REF：参照序列在该位点的基因型ALT：该位点的其余基因型othercoloums：每个个体该位点的基因型(0与REF一致；1与ALT一致；.缺乏数据支持)基因表达水平剖析一个基因表达水平的直接表现就是其转录本的丰度状况，转录本丰度程度越高，则基因表达水平越高。在RNA-seq剖析中，我们能够经过定位到基因组地区或基因外显子区的测序序列(reads)的计数来预计基因的表达水平。Reads计数除了与基因的真切表达水平成正比外，还与基因的长度和测序深度成正有关。为了使不同基因、不同实验间预计的基因表达水平拥有可比性，人们引入了RPKM的观点，RPKM(ReadsPerKilobasesperMillionreads)是每百万reads中来自某一基因每千碱基长度的reads数目。RPKM同时考虑了测序深度和基因长度对reads计数的影响，是当前最为常用的基因表达水平估量方法(Mortazavial.,2008)。

et结果文件分别统计了不同表达水平下基因的数目以及单个基因的表达水平。一般状况下，RPKM数值或许1作为判断基因能否表达的阈值，不同的文件所采纳的阈值不同。表不同表达水平区间的基因数目统计表RPKMIntervalHS1HS2HT1HT2HW1HW20-111678%)11157%)11644%)11552%)11663%)11652%)1-33416%)3829%)3497%)3622%)3359%)3503%)3-156586%)6741%)6719%)6731%)6441%)6522%)15-603436%)3421%)3278%)3277%)3612%)3442%)>601055%)1023%)1033%)989%)1096%)1052%)表基因表达水平统计表geneID

HS1

HS2

HT1

HT2

HW1

HW2ENSGALG00000000003ENSGALG00000000004ENSGALG00000000011ENSGALG000000000138RNA-seq整体质量评估表达水平的饱和曲线检查定量饱和曲线检查反应了基因表达水平定量对数据量的要求。表达量越高的基因，就越简单被正确立量；反之，表达量低的基因，需要较大的测序数据量才能被正确立量。表达水平的饱和曲线的详细算法描绘以下：分别对10%、20%、30%90%的整体测序数据独自进行基因定量剖析，并把全部数据条件下获取的基因的表达水平作为最后的数值。用每个百分比条件下求出的单个基因的RPKM数值和最后对应基因的表达水平数值进行比较，假如差别小于15%，则以为这个基因在这个条件下定量是正确的。图定量饱和曲线检查散布图横坐标代表定位到基因组上的reads数占总reads数的百分比，纵坐标代表定量偏差在15%之内的基因的比率RNA-Seq有关性检查生物学重复是任何生物学实验所一定的，高通量测序技术也不例外(Hansenetal.)。生物学重复主要有两个用途：一个是证明所波及的生物学实验操作是能够重复的且变异不大，另一个为后续的差别基因剖析所需要的。样品间基因表达水平有关性是查验实验靠谱性和样本选择能否合理性的重要指标。有关系数越靠近1，表示样品之间表达模式的相像度越高。Encode计划建议皮尔逊有关系数的平方(R2)大于(理想的取样和实验条件下)。详细的项目操作中，我们要求R2起码要大于，不然需要对样品做出适合的解说，或许从头进行实验。此部分，我们同时计算了spearman有关系数和kendall-tau有关系数作为参照，这两个主假如针对次序变量的有关系数，即秩有关。图RNA-Seq有关性检查R^2:pearson有关系数的平方;rho:spearman有关系数;tau:kendall-tau有关系数均一性散布检查理想条件下，关于RNA-seq技术来说，测序序列(reads)之间为独立抽样并且reads在全部表达的转录本上的散布应当表现均一化散布。但是好多研究表明，好多偏好型的要素都会影响这类均一化的散布(Dohmetal.,2008)。比如，在RNA-seq建库过程中，片段破裂和RNA反转录的次序不相同会致使RNA-seq最后的数据表现严重的3’偏好性。其余要素还包含转录地区的GC含量不同、随机引物等等，并且生物体内从5’或许3’的降解过程相同会致使不均一性散布。图不同表达水平的转录本的reads密度散布图High：高表达量转录本；Medium：中度表达量转录本；Low：低表达量转录本；横坐标为距离转录本5’端的相对地点(以百分比表示)，纵坐标为覆盖深度的均匀值基因差别表达剖析基因表达水平对照经过全部基因的RPKM的散布图以及盒形图对不同实验条件下的基因表达水平进行比较。关于同一实验条件下的重复样品，最后的RPKM为全部重复数据的均匀值。图不同实验条件下基因表达水平比对图RPKM散布图(左图)的横坐标为log10(RPKM),纵坐标为基因的密度。RPKM盒形图(右图)的横坐标为样品名称，纵坐标为log10(RPKM)，每个地区的盒形图对五个统计量(至上而下分别为最大值，上四分位数，中值，下四分位数和最小值)差别表达基因列表基因差别表达的输入数据为基因表达水平剖析中获取的readcount数据。关于有生物学重复的样品，剖析我们采纳DESeq（Andersetal,2010）进行剖析：该剖析方法鉴于的模型是负二项散布，第i个基因在第j个样本中的readcount值为Kij，则有K～NB(μ,σij2)ijij关于无生物学重复的样品，先采纳TMM对re