从RNA-seq结果到差异表达

2011-09-12 ~ ADMIN

翻译自：From RNA-seq reads to differential expression, Oshlack et al. Genome Biology 2010, 11:220

高通量测序技术，也就是下一代测序技术已经成为现代生物学研究的一个较为常规的实验手段了。这一技术的发展极大地推动了基因组学，表观基因组学以及翻译组学的研究。RNA-seq通过测定稳定状态下的RNA样品的序列来对RNA样品进行研究，从而避免了许多之前研究手段的不足，比如象基因芯片或者PCR就需要背景知识。而且RNA-seq还可以触及以前无法研究的领域，比如复杂结构的转录体。RNA-seq可以应用于以下几个方面的研究，

1. SNPs；

2. novel transcripts；

3. alternative splicing；

4. RNA editing。

无论如何，使用RNA-seq最多的还是比较两组样品基因水平表达差异，比如野生型与突变型，用药组与对照组，不同组织之间，癌细胞与正常细胞，等等。我们把这种基因水平差异表达，简称为DE (differential expression，注，不是ED啊〜〜〜)。

常用的RNA-seq操作平台有Illumina GA/ HiSeq, SOLiD 还有Roche 454。它们都是提取RNA后，纯化，打碎，逆转录成cDNA，然后测序。测序的结果被称为short reads，短序。通常一个短序的长度为25-300bp之间。如果测序只测一端可能会带来比对时的困难，于是这些操作平台提供了两端都测的办法，这样的结果成对出现，中间有一定的间隔，但是因为测序长度一下子提高了一倍，所以比对会精准很多。人们把这种测序结果称为’paired-end’ reads，成对短序。一般来讲，测序结果会直接转换成一行一行的由字母组成的短序列，可能是fasta,fastq等等不同格式。

然而，这一技术产生的海量数据分析却给生物学家带来了难题。一个测序的结果文件少则几Gb，多则几十Gb，单独对比拼接，就会用去几个小时，而后再得出差异表达的结果，其耗时耗力，并非实验生物学家可以应付得了的。于是生物信息学的研究人员努力做出一些软件，以降低结果分析的难度。但是，即使这样，还是必须对分析过程有个较为细致地了解，才能正确地使用这些软件，从而得到比较接近事实的结果。

一般的来讲，RNA-seq后DE的工作流程是这样的（图1），首先，将短序映射到基因组相应的位置上去，其次，对映射的结果进行基因水平，外显子水平，以及转录水平的拼接，而后对结果进行数据统计，标准化之后生成表达水平报告文件，最后由生物学者依据系统生物学相关知识，来对数据结果进行分析。

RNA-seq分析工作流程

不同步骤涉汲的软件和方法：

分析步骤	方法	软件
mapping	General aligner	GMAP/GSNAP
		BFAST
		BOWTIE
		CloudBurst
		GNUmap
		MAQ/BWA
		PerM
		RzaerS
		Mrfast/mrsfast
		SOAP/SOAP2
		SHRiMP
	De novo annotator	QPALMA/GenomeMapper/PALMapper
		SpliceMap
		SOAPals
		G-Mo.R-Se
		TopHat
		SplitSeek
	De novo transcript assembler	Qases
		MIRA
Summarization	Isoform-based	Cufflinks
		ALEXA-seq
	Gene-based	Count exons only
		Exon junction libraries
Normalization	library size
	RPKM: reads per kilobase of exon model per million mapped reads	ERANGE
	TMM: trimmed mean of M-values	edgeR
	Upper quartile	Myrna
Differential expression	Poisson GLM (generalized linear model)	DEGseq
		Myrna
	Negative binomial	edgeR
		DESeq
		baySeq
Systems biology	Gene Ontology analysis	GOseq

映射至基因组（Mapping）

第一步的工作是比对(alignment)。对于RNA-seq的比对，从来都不是一件容易的事情。其难点如下：

没有很好的比对模板。现在的比对模板都是基因组模板，而不是真正的转录组模板，也就是说，这对本来就不是很长的短序来说，它很有可能是界于两个exon之间。我们在比对junction的时候，一般还是假设它如果没能在基因组模板中找到合适的位置的时候，才考虑它是否是界于junction上。这种人为的假设可能并不准确。
SNPs，碱基插入，删除，错配，或者质量不高的测序结果，从模板至比对序列本身，都存在着比基因比对更为复杂的问题。
短序可能会有多个100％的匹配位点。
有些基因组可能需要庞大的内存空间。

为了解决最后一个问题，人们使用了很多办法，但基本上都会基于事先建立的引索库。即所谓“启发式”比对(heuristic match)。首先使用一定长度的（通常是11个碱基）的序列做为索引用的关键字，在匹配这一索引字之后，就很大程度地缩小了其需要匹配的模板范围。但是这一办法的问题在于不容易解决问题2中的空格，错配问题。所以在很多软件使用时，会要求人工确认高保真区，以及最高允许2〜3个错配。

现在比较快的“启发式”比对主要有两种算法，一种是哈希表(hash table)，一种是BW压缩转换(Burrows Wheeler transform, BWT)。前者速度快，但是对内存要求比后者要高。

对于问题3，一般而言，大部分软件使用的办法是只保留一个匹配位点，其中，有些是只保留第一个匹配位点，有些是按照概率分布选取保留的位点。当然，前面已经提到过，可以使用paired-end read来尽量避免问题3的出现。

对于问题1，可以使用外显子库来确定junction reads。有两种办法，一种是依靠已知的外显子库来构建，另一种办法就是依据已经匹配好的短序来构建外显子库(de novo assembly of transcriptome)。后者的不足是运算量大，对测序覆盖范围要求高，最好是使用paired-end reads。

还有人发现，对于ploy(A)的处理会减少不能映身的短序数。比如，Pickrell et al.就发现，对于46bp的Illumina reads，87％的短序可以映射至模板，7％可以映射至junction library。如果对那些不能映射的短序，将在头或者尾含有的超过连续4个的A或者T去除，就可以得到约0.005%的映射。

综合评价（Summarizing mapped reads）

这一步，主要是基本于不同水平（外显子水平，转录水平，或者基因水平）进行统计。最简单的办法就是统计落在每个外显上的短序数。但是有研究表明，很多（可能超过15％）的短序会落在外显子两侧，这会影响统计的结果。另一种办法就是统会落在内显子区域的短序数。

无论如何，即使是基因水平的综合评价，也还是有其它的一些问题。比如overlapping的基因的统计。比如junction的统计。

标准化（Normalization）

标准化对于样品内及样品间的比较而言是非常重要的。标准化被分为两类，样品内及样品间（between- and within-library）。

样品内标准化使得在同一样品内不得基因之间的表达差异变得有意义。最常用到的一个办法就是使用落在同一基因内的短序数除以单位基因长度。比较常用的单位是RPKM (reads per kilobase of exon model per million mapped reads)。但是这一方法也受到样品制备和测序方法的干扰。

而对于样品间标准化，最简单而直接的办法使用短序总数来平衡表达量。然而短序总数受测序深度的干扰，而且单个基因的短序数与实际的表达量并不一定会呈线性比较关系。人们又使用四分位(quantile normlization)标准化的办法。但是有研究说这一办法并没有实际的价值。还有提出使用对数分布法则(power law distributions)来进行样品间标准化。但没有研究对这一处理方式进行验证。

差异表达（Differential expression）

差异表达分析的最终目的是将那些差异表达的基因（外显子等等）从海量数据中提取出来。最终的结果显示一般来说是表格化的，这一表格按照一定的规则排序，让人们能够尽可能简单地拿到想要的结果。

由于RNA-seq结果的离散性，人们一般都会使用统计模型来拟合实验得到的结果。一般而言，RNA-seq的结果是比较附合伯松分布(poisson distribution)的。这一结果得到了单通道Illumina GA测序结果的实验验证。但是，伯松分布分析结果常常在多组重复的样品间带来较高的假阳性，因为它低估了生物取样的样品间误差。所以RNA-seq如何设置重复是一个很重要的问题。为了平衡重复样品所带来的误差，人们使用了serial analysis of gene expression (SAGE) data。

现有的软件一般都是针对较为简单的实验设计的。而对于复杂的实验设计，比如说成对样品，时间依赖样品等等，还没有专门的，较好的解决方案。大多数都使用edgeR的线性模型来进行分析。

后期系统生物学分析

简单地讲，前景是广阔的，但目前为止手段还是比较有限的，基本上就是GO分析。