使用salmon和sleuth进行小麦RNA-seq差异表达分析

blast大家一定熟悉吧, 最近有大神爆出一个bug。也即参数“-max_target_seqs”并不是我们理解的那样。下图红色箭头处是NCBI上的设置选项。根据NCBI的官方文件，这个选项是返回前N个最匹配的结果。如果我们设定为1，则表示返回数据库中最匹配的那个。但实际上却不是，该选项意味着返回数据库中第一个匹配良好的结果（simply returns the first good hit found in the database, not the best hit as one would assume ）。据说该bug最早在2015年被发现，但是迄今NCBI还没有任何改动，哪怕是对“-max_target_seqs”参数所表示的实际意思进行修改。

image-20180926112536008

目前的解决办法就是事后处理，即blast运行完之后再进行筛选。

好了，我们言归正传。上次我们在“

评估salmon和kallisto在小麦RNA-seq定量中的异同

”中论证了salmon能够充分区分小麦的同源基因，哪怕只有一个SNP也能够有效区分。同时我们也发现，早期版本的kallisto有重大bug，要尽快升级到最新版本。以前使用kallisto进行定量时，有专门的软件sleuth进行差异表达分析。那么sleuth是何方神圣？2017年sleuth被发表在nature methods。谷歌显示至今大约被引用116次。自然要比我们常用的DESeq2，edgeR等软件的准确率要高。

image-20180926115001362

下面我们就重点谈一谈如何使用这两个软件进行差异基因的表达分析。在实际使用过程中发现里边有很多坑。

第一部分 参考转录组

参考转录组也即中国春参考基因组上注释出来的基因序列。参考转录组当然越全面越准确才好。变好总得需要一个过程。目前1.0版本的数据有如下几个不足：首先，基因的3’和5’端比较短，没有到头；第二是，很多基因具有多个转录本，1.0给的比较少；第三，遗漏了不少转录本。

具体到我们手里的数据，我们可以采用以下做法。也即先找出自己数据中的新转录本或新基因，然后补充进参考转录组进行下一步的差异表达分析。这样有利于发现数据中独特的转录本或基因。当然了，如果只用中国春1.0的参考转录组，对最终的结果影响也不大，例如，功能富集分析的结果等。这些比参考转录组多出来的基因，往往是功能未知基因，没有保守结构域等，可以根据表达等特征筛选比较有趣的进行具体的研究。

下面，我们谈一谈如何获取相对中国春1.0的新转录本和新基因。

1 使用STAR将reads映射至中国春参考基因组。这一步也可以使用其他mapping软件，如hisat2, bowtie2, bbmap等。因为处理的样本较多，我这里使用python写了一个脚本处理。熟悉shell的也可以写shell脚本。

2 筛选bam文件。这里将bam文件里unmapped的reads去掉；因为，我们的样品来自中国春的叶片，所以大于一个mismatch的reads也去掉；去掉mapping长度小于80的reads;只保留proper pair的reads。

顺带交代下样本的背景。这里是对中国春叶片外施ABA后，分别在1h,12h,24h取样。每个时间点3个重复。这是我现编的。

3 组装转录本。输入文件是上面筛选过的bam文件。软件使用scallop。scallop是去年发表在nature biotechnology上的一款有参转录本组装软件。

image-20180927145240879

据作者的评测要比过往的软件好。但这并不意味着使用二代测序reads组装出来的转录本就是无比正确的，这一点要铭记，特别是当要针对具体的某条转录本进行研究时。

On 10 human RNA-seq samples, Scallop produces 34.5% and 36.3% more correct multi-exon transcripts than StringTie and TransComb, and respectively identifies 67.5% and 52.3% more lowly expressed transcripts. Scallop achieves higher sensitivity and precision than previous approaches over a wide range of coverage thresholds.

上述命令运行完之后会产生12个gtf文件。下面要做的就是合并这12个gtf文件中的转录本。这里使用 stringtie merge命令。

4 合并转录本。其中iwgsc_refseqv1.0_HCandLC.gtf来自官方的注释文件。mergelist.txt包含上述12个gtf文件的名字，每行一个，共12行。

5 与官方注释信息比较，找出新增的转录本或基因。

6 获取新增转录本的序列。

7 新转录本和中国春参考转录本合并。

至此新转录本获取完毕。下面选取2个与参考转录本比较下。

上图显示了与参考基因组不同的转录本

新基因，这里也有TGAC的转录本支持

第二部分 基因定量

基因的定量使用salmon这个软件。使用的转录本就是上述我们合并而成的。

#index

salmon index -p 4 -t all_merged.fasta -i all_merged_salmon_index

#定量

for i in ;do salmon quant -i Uall_merged_salmon_index -l A -1 ../clean/$_1_clean.fq.gz -2 ../clean/$_2_clean.fq.gz -p 8 --numBootstraps 100 -o $_quant;done

#参数 --numBootstraps 100 一定要放上

#将生成的以“_quant”结尾的文件夹统一放到salmon_result文件夹下面.

mkdir salmon_result

mv *_quant salmon_result/

至此，定量部分完成

第三部分 差异表达分析

这里使用sleuth。sleuth本来是为kallisto量身定制的，如果salmon想使用sleuth进行差异表达分析，则需要一个R包（wasabi）进行转换。

1 安装wasabi

或者使用 安装

2 结果转换

3 使用sleuth进行差异表达分析

我们这里是按时间点取样，可以按照时序分析，不必进行两两比较。

我们也可以进行两两比较

sleuth还提供了一些画图函数，有兴趣的可以自行尝试。

PCA

samples_density

sleuth结果里只提供了q value 来筛选差异表达基因，没有倍数变化的结果。这里可以使用excel自行进行计算。

计算 Fold change（倍数变化），计算方法如下：

E = mean (group1) B=mean (group2)

FC= (E-B) / min (E, B)

基因 A 和基因 B 的平均值之差与两者中较小的比值。选择 2-3 倍的基因作为差异表达基因即可。可以根据具体数据，斟酌使用。

提供一个小麦里的文献作为参考吧。今年Cristobal Uauy在BMC Plant Biology上以通讯作者发表的文章，题目是“Ubiquitin-related genes are differentially expressed in isogenic lines contrasting for pericarp cell size and grain weight in hexaploid wheat”。

BMC Plant Biology

q value