博文
使用salmon和kallisto进行RNA-seq定量
||
salmon软件于2017年发表在Nature Methods，论文题目是“Salmon provides fast and bias-aware quantification of transcript expression” 。不到两年的时间被引250多次。
根据文中的描述，salmon要优于kallisto和express软件。
今天我们结合小麦最新的IWGSCv1.1版本的基因集，介绍下其定量流程，与kallisto进行一个定量结果的比较，以及看一看基因拷贝数变化时TPM如何变化。
首选要对转录本进行index
salmon index --keepDuplicates -t IWGSC_v1.1_HC_LC_20170706_transcripts.fasta -i IWGSC_v1.1_HC_LC
下载测试数据,该测试数据来自中国春的叶片。
axel -a -n 12 ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR947/SRR947005/SRR947005_1.fastq.gz && axel -a -n 12 ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR947/SRR947005/SRR947005_2.fastq.gz
进行定量
salmon quant -i IWGSC_v1.1_HC_LC -l A -1 SRR947005_1.fastq.gz -2 SRR947005_2.fastq.gz -p 8 -o SRR947005_quant
-i 指要使用的转录组index 这里是我们第一步生成的IWGSCv1.1版本的转录组index；
-l 这里需要指定测序模式，这里-A是指自动判定。
-1 -2 分别要输入左端测序数据和右端测序数据
-p 则是要指定运行程序的cpu线程数
上述定量命令使用测序的原始reads进行的，这里还有另外一种基于比对的方式，即输入文件是BAM或SAM文件。
输出的结果文件如下所示，其中quant.sf里有每个基因的表达量。
salmon定量输出文件kallisto定量过程
首先也是要对转录组进行index
kallisto index -i IWGSC_v1.1_HC_LC_kallisto IWGSC_v1.1_HC_LC_20170706_transcripts.fasta
其次是定量过程。
kallisto quant -i IWGSC_v1.1_HC_LC_kallisto -t 10 -o SRR947005_kallisto SRR947005_1.fastq.gz SRR947005_2.fastq.gz
同样的结果文件如下图所示，其中abundance.tsv里有每个基因的表达量。
kallisto定量结果输出文件salmon文章比较的内容，这里不再比较。下面我们做一个比较初步的分析。
首先看TPM大于>0, >0.5, >1, >10的基因数目。注意两者转录本的总数是一致的。
| TPM | >0 | >0.5 | >1 | >10 |
|---|---|---|---|---|
| Salmon | 109998 | 58690 | 43160 | 7172 |
| kallisto | 116087 | 60805 | 44721 | 7443 |
| kallisto-Salmon | 6089 | 2115 | 1561 | 271 |
从基因数量上来看，kallisto有更多的基因有TPM值，而且相同转录本获得的TPM值，kallisto获得的往往要大于salmon获得的。但这也不是绝对的，如下图。以转录本TraesCS3B02G080100.1为例，这里使用salmon获得的TPM值是338.82，而使用kallisto获得的是168.69。那么针对这一个基因，哪个更准确一点呢？
image-20180902220909122而查阅信息可知，主要是两者的 effecting length 不同导致，而count数量很相似。
| Name | Length | EffectiveLength | TPM(salmon) | NumReads | target_id | length | eff_length | est_counts | tpm(kallisto) |
|---|---|---|---|---|---|---|---|---|---|
| TraesCS3B02G080100.1 | 171 | 35.605 | 338.816942 | 531.769 | TraesCS3B02G080100.1 | 171 | 78.4896 | 532.115 | 168.686 |
This is the computed effective length of the target transcript. It takes into account all factors being modeled that will effect the probability of sampling fragments from this transcript, including the fragment length distribution and sequence-specific and gc-fragment bias (if they are being modeled).
而kallisto里的定义为：
“Effective length” is a scaling of transcript length by the fragment length distribution
具体的可参见“https://haroldpimentel.wordpress.com/2014/05/08/what-the-fpkm-a-review-rna-seq-expression-units/”。
查过一些资料，越看越糊涂。不过，如果使用同一软件处理所有的数据应该是没有问题的。
因为IWGSC官方也有这个数据集的count和TPM值，我去查可一下，竟然和我的不一样。
| COUNT | TPM | |
TraesCS3B02G080100.1 | 294.42 | 96.0559 |
而我这边给的count是532.115，tpm是168.686。事情变得复杂起来。我细想了原因，可能有两方面的原因。其一，我们使用的软件版本不同，其二我们使用的数据SRR947005，他们可能经过一些预处理。而我直接下下来使用，没有经过预处理。
关于软件版本，他们用的是0.42.3，而我用的是kallisto 0.44.0 。所以，接下来我们就要排除软件版本的差异。
To calculate the read counts and TPMs we used Kallisto version 0.42.3. The index was built with the default values and K=31. Kallisto was run with the default parameters. All the fastq files from a single sample are run together.
The reference used was RefSeqv1.1, including the High and Low quality genes (IWGSC_v1.1_ALL_20170706_transcripts.fasta)
结果有点让我接受不了。即如果我使用0.42.3版本，得到的count和TPM就和IWGSC上下载的一样了。但显然与我0.44.0的版本不一致。所以中间出现了什么变故？
| kallisto | target_id | length | eff_length | est_counts | tpm(kallisto) |
|---|---|---|---|---|---|
| 0.44.0 | TraesCS3B02G080100.1 | 171 | 78.4896 | 532.115 | 168.686 |
| 0.42.3 | TraesCS3B02G080100.1 | 171 | 79.8953 | 294.42 | 96.0563 |
有什么变故我们不得而知，但显然让我有点接受不了，我的结果里有多少这样任性的基因？做出来的差异基因有多少又是可信的？下面是作者的回答。软件作者建议使用最新版本。
image-20180903155107203作者都这样说了，加上最新版本的count也与salmon一致，所以，这让本打算使用850 RNA-seq 样本做些分析的我有点迟疑。
第二部分，探究基因拷贝数对定量结果的影响
尽管最近science上转录组那篇文章也有相似的讨论，不过我还是想眼见为实。这个测试也非常简单，那就是人为增加某个基因的拷贝数，观察对这个基因本身的影响以及另外另个同源基因的影响。这里选取TraesCS3D02G074400.1 ，TraesCS3B02G080100.1 和TraesCS1D02G227700.1 。TraesCS3D02G074400.1是孤家寡人。而TraesCS3B02G080100.1则有多个同源基因，TraesCS3B02G080700.1， TraesCS3D02G067400.1, TraesCS3A02G067400.1, TraesCS3B02G080600.1, TraesCS3D02G067100.1, TraesCS7B02G443100.1, TraesCS3D02G067200.1, TraesCS3D02G219800.1, TraesCS3A02G233600.1 ；而TraesCS1D02G227700.1对应的有TraesCS1B02G240400.1，TraesCS1A02G226300.1,同时该基因还有可变剪切，TraesCS1D02G227700.2，TraesCS1B02G240400.2，TraesCS1A02G226300.2。
所以我们比较的重点就放在这几个基因上，看这几个基因的变化。
kallisto index -i IWGSC_v1.1_HC_LC_kallisto_add IWGSC_v1.1_HC_LC_add.fasta
kallisto quant -i IWGSC_v1.1_HC_LC_kallisto_add -t 10 -o SRR947005_kallisto_add SRR947005_1.fastq.gz SRR947005_2.fastq.gz
首先对这三个基因来说，其表达量减半。因为我们在序列上一字未改的复制了这3个序列。所以从结果看，复制的基因分享了一半的reads。
| target_id | length | eff_length | est_counts | tpm(kallisto) |
|---|---|---|---|---|
| add|TraesCS3B02G080100.1 | 171 | 79.2217 | 266.054 | 83.6109 |
| TraesCS3B02G080100.1 | 171 | 79.2217 | 266.054 | 83.6109 |
| add|TraesCS3D02G074400.1 | 752 | 561.072 | 4175.4 | 185.275 |
| TraesCS3D02G074400.1 | 752 | 561.072 | 4175.4 | 185.275 |
| add|TraesCS1D02G227700.1 | 1867 | 1675.73 | 284.687 | 4.22962 |
| TraesCS1D02G227700.1 | 1867 | 1675.73 | 284.687 | 4.22962 |
同样的我们再比较多拷贝的情况。多拷贝的我们看这个基因TraesCS3B02G080100。复制之后，TraesCS3B02G080100的表达减半，其他拷贝几乎不变。
那么是否对可变剪切的转录本有影响呢？结果和上面一样，只是TraesCS1D02G227700.1自身TPM值减半。
既然拷贝数增加，只改变其自身的TPM值。那么同样的道理，如果减少一个拷贝数，其他拷贝数的TPM也不会发生显著变化。当然这里是有一个前提，那就是不同拷贝之间序列不是100%一致的。
刚才的复制我们一字未改，现在我们改变其中一个碱基，我们再观察结果如何。统一在这3个转录本序列的第80位改变一个碱基。理论上，改变之后这3个转录本的TPM值为0。
Kallisto的结果如下
| target_id | length | eff_length | est_counts | tpm(kallisto) |
|---|---|---|---|---|
| 1snp|TraesCS3D02G074400.1 | 752 | 561.072 | 0 | 0 |
| 1snp|TraesCS3B02G080100.1 | 171 | 79.2217 | 0 | 0 |
| 1snp|TraesCS1D02G227700.1 | 1867 | 1675.73 | 0.698188 | 0.010373 |
| TraesCS1D02G227700.1 | 1867 | 1675.73 | 568.676 | 8.44886 |
| TraesCS3B02G080100.1 | 171 | 79.2217 | 532.107 | 167.222 |
| TraesCS3D02G074400.1 | 752 | 561.072 | 8350.8 | 370.551 |
如上表所示，实际情况也差不多如此，只有1snp|TraesCS1D02G227700.1有很小的TPM。
在改一个SNP的情况下，salmon表现似乎更好，修改的全部是0，而原来的TPM值也不变。
| Name | Length | EffectiveLength | TPM | NumReads |
|---|---|---|---|---|
| 1snp|TraesCS3D02G074400.1 | 752 | 550.209 | 0 | 0 |
| 1snp|TraesCS3B02G080100.1 | 171 | 35.447 | 0 | 0 |
| 1snp|TraesCS1D02G227700.1 | 1867 | 1665.119 | 0 | 0 |
| TraesCS1D02G227700.1 | 1867 | 1665.119 | 7.729668 | 567.493 |
| TraesCS3B02G080100.1 | 171 | 35.447 | 340.231283 | 531.757 |
| TraesCS3D02G074400.1 | 752 | 550.209 | 343.102088 | 8323.484 |
所以总结来说，在同一个项目中，使用的软件版本要一致。至于哪一个salmon和kallisto哪个更好，我个人可能更偏向salmon，但kallisto毕竟是science使用的，所以做小麦的话，kallisto应该是公认的。kallisto算出来的差异基因数量要多于salmon。
但同时也要注意IWGSC注释的转录本的正确率。特别是对于低可信度的基因，基因结构不正确，做出来的TPM可能也不对。而且IWGSC注释的转录本其3’和5’的非翻译区有不少缺失。转录组里还可以组装出新的转录本或者新的可变剪切。
目前IWGSC给出的转录本是26万左右，而其中高可信的基因只有10万。我们进行差异表达分析时使用的是26个转录本，但只有那10万个高可信基因的功能注释较为可信。所以后续的功能富集分析，就要慎重。差异表达可以包括低可信的基因，但功能富集分析时可以不包括在里面，特别是不能用来充当总数，science上的那篇文章也正是这样做的。另外，功能富集分析时使用的背景基因应该是在某条件下表达的基因集，而不是总基因集。比如，总共有10个基因，其中只有5个基因表达，而这5个表达的基因里有2个差异表达，此时使用的总数是5而不是10。当然，如果差异基因里包括低可信的基因，那表达的总基因里也可以包括低可信度的基因。不过，我还是建议使用表达的高可信度的基因作为整体。
昨天我在推送里表达了在基因富集分析中背景基因如何选择的观点？我的观点是需要使用表达的基因，确切的说是使用表达的并且有功能注释的基因作为背景。比如，做GO分析，背景基因都得有GO注释信息，那些没有GO注释信息的显然不宜拿来当作背景。同样的拿来做GO分析的差异表达基因也得都GO注释信息，那些差异表达但没有GO注释信息的基因就不必包括在里边。为啥？道理很简单，我要做的是GO信息的富集分析，那些没有GO注释信息的基因即使放在里边也没毛用。
这是我的观点，也许不对。如果有大神跳出来吊打，我还是欢迎的。群里也有小伙伴发表自己的看法，看过之后很受启发。接下来，且听我细细道来。
这里所说的表达的基因，是指只要在分析所用样本中的一个里有表达（如TPM>0.5就认为表达），比如在叶片和根里，只要基因在这两个组织中的任何一个里表达，我们就认为基因是表达的。即取的是不同样本间的并集。这里只是拿叶片和根说事，如果在文章里有人这样做，你可以将这篇文章扔到垃圾桶了。
同时我在google上也查了一些资料。在biostars（https://www.biostars.org/p/17628/）也有很多人在讨论。
基本的结论是，使用表达的基因作为背景是合适的，没有人明确反对这一做法不妥。
根据上面的回答，找到这篇文献“Multiple sources of bias confound functional enrichment analysis of global -omics data”。该文2015年发表在genome biology上。
文中的结论是：
The generation of an estimated background ‘universe’ in RNA-seq data could be achieved by removing zero-count genes, but the nature of this ‘universe’ will still depend on many factors.
This list should contain only factors (RNA or protein) that are both robustly ‘probed or sequenced’ (to avoid technology and detection bias) and ‘called’ as expressed (to avoid biological bias) in the experiment. After all, we do not want to carry out functional enrichment analysis of a specific tissue simply to be informed that we are studying that tissue!
In conclusion, functional enrichment analysis must not be considered proof of biological plausibility or validity in the analysis of high-throughput -omics data. We strongly advocate for efforts to generate appropriate background expression ‘universes’. We also urge that background gene lists are provided for any functional enrichment analysis, and that a higher statistical threshold is used as a default, given the scale of the pre-existing biases, to avoid marginal (e.g., 1 × 10^−3^ or 1 × ^10−4^) enrichments being relied on to drive the interpretation of an experiment.
做差异基因功能富集分析的一定要看看本文。里边讨论了各种有可能导致最终结果不准确的情况。
一些在线进行GO富集分析的网站，如果不能自己背景序列，那就是耍流氓。
However, somewhat surprisingly, the Gene Ontology consortium website analysis tool does not allow for this option, making any analysis thoroughly unreliable.
实际上任何实验手段都不是完美的。我们现在做RNA-seq比较方便，RNA-seq并不是万能呢，蛋白水平就要比转录水平要好。还有组织毕竟是混杂，那单细胞结果又如何呢？代谢组呢？
通过差异表达基因功能富集结果我们可以说参与了某些通路，或者在哪些通路上富集。但最好不要说与哪些通路无关，但可以说这些通路上的基因没有富集。
https://wap.sciencenet.cn/blog-1094241-1133526.html
上一篇：大麦铝耐受性基因HvAACT1研究
下一篇：使用salmon和sleuth进行小麦RNA-seq差异表达分析
