RNA-seq数据分析
我们对工具和技术的评估落脚点在它们鉴定出的差异基因的准确性。为了完成这个评估,至少需要四个不同的分析阶段(图2;表2)。第一阶段把测序平台生成的原始测序数据比对到转录组。第二阶段量化与每个基因或转录本来源的reads数量,构建表达矩阵。该过程可能包括1个或多个子过程如比对,组装和定量,或者它也可以一个从读取计数生成表达矩阵。通常有一个第三阶段,包括过滤低表达的基因和至关重要的移除样品间技术差异的标准化过程。DGE的最后阶段是构建样本分组和其它协变量的统计模型,计算差异表达置信度。
第1阶段-测序reads的比对和组装
测序完成后,分析的起点是包含测序碱基的FASTQ文件。最常见的第一步是将测序reads比对到已知的转录组(或注释的基因组),将每个测序reads转换为一个或多个基因组坐标。传统上,该过程是通过几个不同的比对工具(如TopHat,STAR或HISAT)完成的,其都依赖参考基因组的存在。由于测序的cDNA来自RNA,可能跨越外显子边界,因此与参考基因组(包含内含子和外显子)比对时需要进行剪接比对,即允许reads中出现大片段gap。
如果没有可用的包含已知外显子边界的高质量基因组注释,或者如果希望将reads与转录本(而不是基因)相关联,则需要在比对后执行转录组组装步骤。诸如StringTie和SOAPdenovo-Trans之类的组装工具使用比对reads的gap来推测外显子边界和可能的剪接位点。转录本重头组装特别适用于参考基因组注释缺失或不完整的物种,或者对异常转录本感兴趣(例如在肿瘤组织中)的研究。转录组组装方法受益于双端测序和/或更长的reads的使用,增加跨越splice junctions的可能性。但是,通常不需要从RNA-seq数据中从头做转录组组装来确定DGE (生信宝典注:无参分析组装是必须的)。
最近,涌现了一些计算效率高的“alignment free”工具,例如Sailfish,Kallisto和Salmon,它们将测序reads直接与转录本关联,而无需单独的定量步骤。这些工具在定量高丰度(以及长度更长)的转录本方面表现出很好的性能。但是,它们在定量低丰度或短转录本方面不够准确。(39个工具,120种组合深度评估 (转录组分析工具哪家强))
第2阶段-定量转录本丰度
将reads比对到基因组或转录组后,下一步就是将它们分配给基因或转录本,获得表达矩阵。不同的比较研究表明,定量过程中采用的方法对最终结果的影响最大,甚至比比对工具影响更大。
常用的定量工具包括RSEM,CuffLinks,MMSeq和HTSeq,以及上述的无比对直接定量工具。基于reads计数的工具(例如HTSeq或featureCounts)通常会丢弃许多比对的序列,包括那些具有多个匹配位置或比对到多个表达特征的reads。这可以在随后的分析中消除同源和重叠的转录本。RSEM使用期望最大化模型来分配模糊的reads,而无参考的比对方法(例如Kallisto)则将这些reads用于后续的定量,这可能会导致结果偏差。转录本丰度估计可以转换成等效的read计数,能完成这一转换的部分工具依赖tximport包。量化步骤结束后会得到一个合并的表达矩阵,每个表达特征(基因或转录本)各占一行,每个样品各占一列,中间的值是实际读数 (reads count)或估计的表达丰度。
阶段3-过滤和标准化
通常,基因或转录本的reads count需要进行过滤和标准化,以移除测序深度、表达模式和技术偏差的影响。过滤去除在所有样本中都低丰度表达的基因是很直接的方式,标准化表达矩阵的方法要复杂一些。简单的转换可以校正丰度,降低GC含量和测序深度的影响
阶段4-差异表达分析
获得表达矩阵后,就可以构建统计模型评估哪些转录本发生了显著的表达改变。有几个常用工具可以完成此任务;一些基于基因水平的表达计数,其它的基于转录本水平的表达计数。基因水平的工具通常依赖于比对的reads计数,并使用广义线性模型来进行复杂实验设计的评估。这些工具包括EdgeR,DESeq2和limma + voom等工具,这些工具计算效率高并且彼此之间结果稳定性好。评估差异isoforms表达的工具,例如CuffDiff,MMSEQ和Ballgown,往往需要更多的计算资源,并且结果的变化也更大。但是,在差异表达工具应用之前的操作(即关于比对、定量、过滤和标准化)对最终结果的影响更大。
参考基因组和基因注释文件获取
通常测序生成的reads要与参考基因组或参考转录组进行比对,或Pseudo-alignment。所以首先需要获取参考基因组和参考转录组信息。
Ensembl http://www.ensembl.org/info/data/ftp/index.html 是常用的信息齐全的参考基因组和GTF文件下载网站。
Ensembl提供的参考基因组有2种组装形式和3种重复序列处理方式, 分别是primary, toplevel和unmasked (dna)、soft-masked (dna_sm)和masked (dna_rm)。一般选择dna.primary或dna_sm.primary。
为什么选择Primary
Primary assembly contains all toplevel sequence regions excluding haplotypes and patches. This file is best used for performing sequence similarity searcheswhere patch and haplotype sequences would confuse analysis.
为什么不选择masked
Masked基因组是指所有重复区和低复杂区被N代替的基因组序列,比对时就不会有reads比对到这些区域。一般不推荐用masked的基因组,因为它造成了信息的丢失,由此带来的一个问题是uniquely比对到masked基因组上的reads实际上可能不是unique的。而且masked基因组还会带来比对错误,使得在允许错配的情况下,本来来自重复区的reads比对到基因组的其它位置。 另外检测重复区和低复杂区的软件不可能是完美的,这就造成遮盖住的重复序列和低复杂区并不一定是100%准确和敏感的。
soft-masked基因组是指把所有重复区和低复杂区的序列用小写字母标出的基因组,由于主要的比对软件,比如BWA、bowtie2等都忽略这些soft-mask,直接把小写字母当做大写字母比对,所以使用soft-masked基因组的比对效果和使用unmasked基因组的比对效果是相同的。
基因注释GTF文件在分析转录组数据时会用到,也从这获取,GTF文件的解释见文件格式部分。ENSEMBL的基因注释文件与GeneCode(http://www.gencodegenes.org/)V26版本一致
下载基因功能和结构注释信息
ENSEMBL数据库的BioMart http://www.ensembl.org/biomart/martview工具为下载基因的功能信息、序列信息、结构信息、ID的转换等提供了很大的便利。
注意在BioMart的Attribute选项里如果选择了蛋白相关的选项,得到的结果中只有蛋白编码基因的信息。如果要下载所有基因信息,请不要选择蛋白相关的选项。
GTF/GFF文件格式解读和转换
测序原始数据下载
生物或医学中涉及高通量测序的论文,一般会将原始测序数据上传到公开的数据库,上传方式见测序文章数据上传找哪里;并在文章末尾标明数据存储位置和登录号
NCBI的SRA (Sequence Read Archive) 数据库(http://www.ncbi.nlm.nih.gov/sra/) 是最常用的存储测序数据的数据库。目前SRA数据的组织方式分为下面4个层次:
Studies—研究课题;
Experiments—实验设计;
Runs—测序结果集;
Samples—样品信息。
使用NCBI提供的SRA-toolkit中的工具fastq-dump直接下载SRR文件,并转换为FASTQ格式,--split-3参数表示如果是双端测序就自动拆分,如果是单端不受影响。--gzip转换fastq为压缩文件,节省空间。