关于转录本定量

刘小泽写于19.2.9
转录本定量一般有两种方式，一个是基于比对定量（如RSEM、eXpress）；另一种是不基于比对的方法（转录本划分成kmer, 用kmer出现的次数来衡量转录本丰度，如kallisto、salmon，优点就是快）

注意一：如果自己有多组数据（例如同一物种的不同组织），最好先全部拼接成一个转录本，然后再分开进行定量

注意二：Trinity软件虽然可以调用RSEM、eXpress等，但是并不是内置，还是需要我们自己先手动安装好，放入环境变量；或者直接conda安装

Trinity为了定量，推出了一个脚本align_and_estimate_abundance.pl

帮助文档很详细，需要注意的就是： 尽量使用--samples_file这个参数，因为它会将不同的处理及重复的结果文件自动放入不同的目录；--prep_reference这个参数是很有帮助的，它先对拼接的转录本构建索引，然后再对不同样本进行并行计算表达量，可以视作一个提速的手段；--trinity_mode或者--gene_trans_map除了得到isoform的表达矩阵，还可以得到gene层面的表达矩阵【前面提到过，拼接的转录本中一个gene对应多个isoform】

怎么定量？

基于比对的模式生成的结果中会有FPKM（fragments per kilobase transcript length per million fragments mapped）、TPM（transcripts per million transcripts）两种标准化方法；不基于比对的结果只有TPM。

所谓的标准化就是：去除测序深度和基因长度的的影响

这两种哪种更合理呢？关于FPKM和TPM之争，许多文章都有介绍：Wagner et al. 2012. Theory Biosci ， Li and Dewey, BMC Bioinf. 2011. 目前一般都推荐TPM

首先，回归主要问题：定量。我们来看看什么是表达量？

• 如果要明确知道一个细胞中、或者是一定摩尔量的 RNA中，有多少条某种转录本，就是绝对定量（如果能够对样品进行细胞计数或者添加 spike-in, 转录组也可以实现绝对定量）
• 但我们平时做的转录组测序，要求并没那么高，并不知道用于建库测序的 RNA 来自多少个细胞，总共有多少转录本。因此只能相对定量，也就是描述某基因的转录本占样本中所有转录本的百分比。

好了，前面说的定量就是统计的read count数（下图左一），但是这还并不能说明真实情况，因为没有考虑长度的分布（比如下图的3和4：在read count总数中确实4比3要多许多，但是注意观察，3是两段短的转录本，4是一个长一短。有没有注意到，其实3和4的read分布是一个水平的，它们都很均匀，和1相比1就显得比较稀疏）。为了更好了了解表达量（Expression）和片段数（read count）的差别，我们再来看看下面图中的2和4，总量的话（左一）2和4是一样的，但是2明显每个转录本分布的reads更多，因此换算到表达量时，2的表达量要比4高不少（右一）

raw-count and expression

上面知道了要同时考虑raw count和转录本长度的影响，因此就需要公式来进行标准化，公式很多，那么哪一种更能反映真实问题呢？

对比

对公式不清楚的话，其实简单看发表年份就能知道，TPM比FPKM更先进一些

下面具体从公式角度了解，看看两种方法的计算公式：

TPMvsFPKM

Xi表示比对到某基因上的Fragment数目，单位是Millions Fragments；Li表示这个基因外显子的长度，单位是Kilobase。因此Xi/Li 的商就是外显子长度为1Kb的基因Fragment数目，即转录本丰度，因此TPM可以校正转录本长度分布的影响，用于可以更好地进行样本间比较。

【注意，二者分母不同，TPM的分母中考虑了长度的影响，可以理解成总的转录本丰度；而FPKM中，没有考虑长度影响，可以理解成总的Fragment数目】

如果对每个样品总 表达量求和，会发现 TPM中各个样本的总和是相等的，而FPKM则不相等。因此，能用TPM就不用FPKM

注意最后加和

【目前DESeq2和edgeR不用RPKM/FPKM或TPM做均一化，而是直接用原始的read counts做均一化处理】

参考：https://haroldpimentel.wordpress.com/2014/05/08/what-the-fpkm-a-review-rna-seq-expression-units/

https://www.youtube.com/watch?v=TTUrtCY2k-w

定量结果都有啥？

基于比对的定量

比对过程会产生bowtie.bam，这个文件会直接导入RSEM或eXpress，如果定量时选择--coordsort_bam参数，会对bam文件进行sort，生成一个bowtie.csorted.bam文件，为了方便IGV可视化

• RSEM（RNA-Seq by Expectation Maximization）：结果生成两个文件，均包含定量信息：

  结果产生两个文件
  RSEM.isoforms.results : 每个拼接转录本的最大期望count数[见下表]
  RSEM.genes.results : 一个基因对应多个转录本(此处基因非严格意义上的基因)，统计了每个"基因“的最大期望count数
  

  transcript_id	gene_id	length	effective_length	expected_count	TPM	FPKM	IsoPct
  TRINITY_DN100_c0_g1_i1	TRINITY_DN100_c0_g1	443	181.37	4.84	1670.06	12311.85	100.00
  TRINITY_DN101_c0_g1_i1	TRINITY_DN101_c0_g1	251	19.37	1.00	3231.22	23820.87	100.00
  TRINITY_DN103_c0_g1_i1	TRINITY_DN103_c0_g1	1219	957.37	0.00	0.00	0.00	100.00
  TRINITY_DN103_c0_g2_i1	TRINITY_DN103_c0_g2	414	152.41	0.00	0.00	0.00	100.00

  关于RSEM：https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-12-323

• eXpress：

  结果产生两个文件(10 tab-delimited columns, sorted by the bundle_id)
  
  results.xprs  : transcript abundance estimates (generated by eXpress)
  results.xprs.genes : gene abundance estimates (provided here by summing up transcript values per gene)
  

  包含了： estimated counts, FPKM, and TPM 等信息

  关于eXpress：https://pachterlab.github.io/eXpress/overview.html#top

  帮助文档：https://pachterlab.github.io/eXpress/manual.html

不基于比对的定量

• Kallisto：

  结果产生两个文件
  abundance.tsv  : transcript abundance estimates (generated by kallisto)
  abundance.tsv.genes : gene abundance estimates (provided here by summing up transcript values per gene)
  

  输出结果短小精悍

  target_id	length	eff_length	est_counts	tpm
  TRINITY_DN10_c0_g1_i1	334	100.489	13	4186.62
  TRINITY_DN11_c0_g1_i1	319	87.9968	0	0
  TRINITY_DN12_c0_g1_i1	244	38.2208	2	1693.43
  TRINITY_DN17_c0_g1_i1	229	30.2382	5	5351.21

  https://pachterlab.github.io/kallisto/manual

• Salmon：

  结果产生两个文件
  quant.sf : transcript abundance estimates (generated by salmon)
  quant.sf.genes : gene-level abundance estimates (generated here by summing transcript values)
  

  结果也是很简单

  Name	Length	EffectiveLength	TPM	NumReads
  TRINITY_DN10_c0_g1_i1	334	135.084	0	0
  TRINITY_DN11_c0_g1_i1	319	120.926	6158.51	13
  TRINITY_DN12_c0_g1_i1	244	57.9931	0	0
  TRINITY_DN17_c0_g1_i1	229	47.8216	2395.84	2

组合定量结果

因为之前是每个sample生成一个结果文件，现在需要组合在一起，使用abundance_estimates_to_matrix.pl

需要指定的参数：
          
--est_method <string>  RSEM|eXpress|kallisto|salmon  (needs to know what format to expect)需要和之前定量方法一致才能识别并组合

--gene_trans_map <string>  the gene-to-transcript mapping file. (if you don't want gene estimates, indicate 'none'）. 

例如：

ls */RSEM.genes.results >genes.quant_files.txt
$TRINITY_HOME/util/abundance_estimates_to_matrix.pl 
    --est_method RSEM \
    --cross_sample_norm TMM \
    --out_prefix rsem-gene \
    --name_sample_by_basedir \ 
    --quant_files genes.quant_files.txt

得到结果

rsem-gene.isoform.counts.matrix  : the estimated RNA-Seq fragment counts (raw counts)

rsem-gene.isoform.TPM.not_cross_norm  : a matrix of TPM expression values (not cross-sample normalized)

rsem-gene.isoform.TMM.EXPR.matrix : a matrix of TMM-normalized expression values

其中counts.matrix文件是用来下游差异表达分析的 ， TMM.EXPR.matrix文件是基因表达矩阵

关于TMM的重要性： Robinson & Oshlack, Genome Biology 2010 and [Dillies et al., Brief Bioinf, 2012

---

欢迎关注我们的公众号～_～　　
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到Bioplanet520@outlook.com

Welcome to our bioinfoplanet!