RNASeq测序数据分析服务流程试运行精.doc

北京大学生科院/CLS生物信息平台RNA-Seq测序数据分析服务流程(试运行2015.3平台联系人:李程(lch3000文档撰写:张超Table of Contents1. 测序质量评估 (31.1 测序数据过滤 (31.2 质量值分布 (31.3 GC含量分布 (42. 参考序列比对 (43. 基因表达水平 (63.1 基因表达水平定量 (63.2 基因表达水平分步 (63.3 生物学重复相关性分析 (63.4 样本间层次聚类及PCA分析 (74. 差异基因分析 (74.1 基因表达标准化 (74.2 差异基因列表 (84.3 差异基因可视化 (84.4 差异基因聚类 (95. 差异表达基因功能分析 (105.1 GO富集分析 (105.2 信号通路富集分析 (105.3 癌基因功能注释 (116.基因结构差异分析 (116.1 可变剪切分析 (117. SNP分析 (127.1 SNP检测 (127.2 SNP 筛选 (127.3 GO/KEGG富集 (121. 测序质量评估通过测序的数据进行进行质控,保证数据质量适合下游分析。这里我们使用fastqc和RNA-SeQC来对数据进行质量评定。1.1 测序数据过滤测序得到的原始下机数据往往有许多问题,不能直接使用,通常会经过以下过滤,尽量保证测序数据的质量。a.去除带测序接头的测序序列(reads;b.去除低质量的reads1.2 质量值分布按照现有的测序技术(illumina平台单碱基的错误率应控制在1%以下,即质量值在20以上。 横坐标为reads的碱基位置,纵坐标为单碱基质量值质量值与错误率的关系:Q=-10log10(e;其中Q phred为测序碱基质量值,e为测phred序错误率。1.3 GC含量分布对于RNA测序,鉴于序列通过超声随机打断,所以理论上每个测序循环上的C、G及A、T含量应分布相等,并且CG-content对于每个物种应大致相同。 横坐标为reads的碱基位置,纵坐标为各种碱基的不同比例2. 参考序列比对对于通过质量控制的数据,可以进行后续分析。首先需要将clean reads比对到参考基因组上。由于测序时reads是随机的,只有这些reads的碱基信息和质量信息,没有其在基因组上的位置信息,比对这一步就是给所有reads一个在基因组上位置的信息。在RNA测序中,其实测的是cDNA的序列,由于内含子的存在,所以会较常出现一条read跨内含子的情况,tophat2可以较好的处理这种情况,所以我们选用tophat2来做比对。比对率间接反应了测序的质量和建库的质量,若比对率低,很可能建库时混入了其他物种的序列,导致无法比对到研究的物种参考基因组上。 reads比对到基因上的位置统计:Sample IntragenicRateExonicRateIntronicRateIntergenicRateSplitReadsExpressionProfilingEfficiencyTranscriptsDetectedGenesDetected1BJ 0.885 0.738 0.147 0.114 9,910,010 0.738 32,796 15,434(1Sample:样本名(2IntragenicRate:比对到基因内的reads比例(3ExonicRate:比对到外显子的reads比例(4IntronicRate:比对到内含子的reads比例(5IntergenicRate:比对到基因间区的reads比例(6SplitReads:比对到两外显子交接处的reads数(7ExpressionProfilingEfficiency:比对到外显子上的reads占总体的比例(8TranscriptsDetected:比对上reads数大于5的转录本数(9GenesDetected:比对上reads数大于5的基因数 3. 基因表达水平3.1 基因表达水平定量在RNA-seq分析中,我们可以通过定位到基因组区域或基因外显子区的reads的计数来估计基因的表达水平。Reads计数除了与基因的真实表达水平成正比外,还与基因的长度和测序深度成正相关。为了使不同基因、不同实验间估计的基因表达水平具有可比性,人们引入了RPKM的概念,RPKM(Reads Per Kilo bases per Million reads是每百万reads中来自某一基因每千碱基长度的reads数目。RPKM同时考虑了测序深度和基因长度对reads计数的影响,是目前最为常用的基因表达水平估算方法 (Mortazavi et al., 2008。Gene_ID Sample1 Sample2 Sample3 Sample4 Sample5 Sample6 ENSG00000000003 49.32 46.94 48.91 22.51 20.60 22.95 ENSG00000000419 35.92 34.58 33.69 32.80 35.65 32.73 ENSG00000000457 1.34 0.94 1.19 2.06 2.13 2.26 ENSG00000000460 1.19 1.20 1.22 3.00 3.33 3.06(1 Gene_ID:Ensembl基因ID(2 Other columns:各样本中该基因的表达水平(RPKM3.2 基因表达水平分步每个样本所有基因的RPKM盒形图可以展示出不同实验条件下基因表达水平的分布情况。 图3.2.1 不同条件下的基因表达水平分布图3.3 生物学重复相关性分析生物学重复主要有两个用途:一个是证明所涉及的生物学实验可重复性强、差异小,另一个用于估计生物学变异进行差异基因检测。样品间基因表达水平相关性是检验实验可靠性和样本选择是否合理的重要指标。相关系数越接近1,表明样品之间表达模式的相似度越高。 图3.3.1 生物学重复散点图3.4 样本间层次聚类及PCA分析当样本数目较多时,可以利用基因的表达量进行样本间聚类分析及PCA 分析,对样本间关系进行探究或者对实验设计进行验证。样本聚类距离或者PCA距离越近,说明样本越相似。 4. 差异基因分析4.1 基因表达标准化对于有生物学重复的样品,我们采用DESeq2提出的scaling factor的方法对原始的readcount进行标准化(normalization。以消除非生物学引起的readcount的差异(最主要消除各个文库测序数据量不同带来的差异。对于标准化的结果,我们采用MA-plot或box-plot来评价。图4.1.1 MA-plot 横坐标为表达量,纵坐标为log后的表达差异倍数基于大部分基因都是非差异表达的,所以大多点应在log fold change=0左右,并且不随表达量的变化而变化。4.2 差异基因列表对于有生物学重复的的样品,我们采用DESeq2来分析差异表达基因。该方法基于负二项分布模型(K ij NB(ij,ij2来检测差异表达基因。Gene baseMean log2FoldChange pvalue padjFBgn0000370 31324.379200 -1.3665378519 5.6393206e-176 2.9843284e- 172FBgn0033913 17544.483454 -1.1571536021 6.3177309e-90 1.3373372e-87(1Gene: 基因ID(2baseMean:所有样本矫正后的平均reads数(3log2FoldChange:log2后的表达量差异(4pvalue:统计学差异显著性检验指标(5padj:校正后的pvalue。padj越小,表示基因表达差异越显著4.3 差异基因可视化用火山图可以推断差异基因的整体分布情况。 图 4.3.1显著性差异表达基因用红色点表示;横坐标表示基因在不同样本中的表达倍数变化;纵坐标表示统计学上基因表达量变化差异的显著性对于特定基因在不同实验中的表达情况,和此基因的不同转录本在不同实验中的表达情况。 图 4.3.2 左图为regucalcin基因在两个样本中的表达差异情况;右图为此基因在不同转录本中的表达差异情况4.4 差异基因聚类聚类分析用于判断差异基因在不同实验条件下的表达模式。通过将表达模式相同或相近的基因聚集成类,从而识别未知基因的功能或已知基因的未知功能。 5. 差异表达基因功能分析5.1 GO富集分析Gene Ontology(简称 GO, http:/www.geneontology.org/是基因功能标准分类体系。研究差异基因在 Gene Ontology 中的分布情况将阐明差异基因富集的生物学功能。 5.2 信号通路富集分析在生物体内,不同基因相互协调实现其生物学功能,通过Pathway显著性富集能确定差异表达基因参与的最主要信号通路。KEGG(Kyoto Encyclopedia of Genes and Genomes,http:/www.kegg.jp/是有关Pathway的主要公共数据库(Kanehisa,2008。Pathway显著性富集分析以KEGG Pathway为单位,应用统计检验找出差异表达基因中显著性富集的Pathway。 5.3 癌基因功能注释 原癌基因(Proto-oncogene是参与细胞生长、细胞分裂和细胞分化的正 常基因，当其发生突变后(如基因序列被改变就会变成致癌基因(Oncogene。 通常在肿瘤或恶性细胞系中某些特异性癌基因会上调表达，通过了解癌基因在 实验不同组的表达情况有助于深入认识疾病的发病机理。 Cosmic( https:/cancer.sanger.ac.uk/cosmic 是英国 Sanger 实验室开发并维 护的癌基因及相关注释数据库，有较高的权威性及可信度，通过与数据库进行 比对，可对差异表达基因中的癌基因部分进行鉴别和注释。 6.基因结构差异分析 6.1 可变剪切分析 对于 RNA-seq，除了 gene 水平的差异分析外，还可以进行 exon 水平 的差异分析。不用的 exon 表达，表明了有着不同的剪切方式。这时可以使用 Bioconductor 的 DEXSeq 软件包。 该分析可以给出每个基因在不同的实验条件下，外显子的使用情况。比 如上图的 10 号外显子在 control 和 knockdown 两组中的表达差别较大，此外 显子的表达量情况，也反映到了在两组中此基因的剪切形式有差异。 7. SNP 分析 7.1 SNP 检测 SNP 全称 Single Nucleotide Polymorphisms，是指在基因组上由单个核 苷酸变异形成的遗传标记，其数量很多，多态性丰富。一般而言，SNP 是指变 异频率大于 1%的单核苷酸变异。对 RNA-seq 的 SNP 分析可能得到基因在上 的 SNP 位点和 RNA 编辑位点。 Chr Chr1 Chr1 Chr1 Pos 14653 14907 14930 Ref C A A Alt T G C 7.2 SNP 筛选 对 SNP 位点进行注释、过滤、筛选，旨在找出跟表型相关性高的位点。 过滤 dbSNP 中存在的多态性位点，过滤掉同义突变的位点。 与现有 GWAS 位点数据库比对，与 OMIM 数据库，HGMD 数据库等比 对。 7.3 GO/KEGG 富集 对高可信的 SNP 所在的基因进行 GO/KEGG 富集。