RNA-seq从入门到自闭(Kallisto和Salmon)
这是RNA-seq上游分析的最后一站,seq数据定量。这一篇文章会介绍基于k-mer定量两软件:kallisto和salmon。其中关于kallisto的部分我会附上TBtools插件的用法。
抱歉又更新晚了,之前一直想尝试selected alignment method来定量RNA-seq数据。电脑不给力,试了好几次都失败了,只好放弃……
1. 什么是k-mer
简单来说,k-mer是一段长度为k的序列,而后面的mer即为monomeric unit(单体单元),也就是每个碱基。因k-mer包含k个碱基,若一段核酸序列长度为L,以一个碱基为步长滑动,那么根据这个核酸序列就可以得到L-k+1个k-mer;由于每个位点的碱基可以为(A、T、C、G)中的任意一个,因此k-mer理论上说有个不同的序列。原本一条长片段,就变成了很多短的片段,因此计算机处理的碱基数量也会增加很多倍。而且,每次取k-mer是同一条reads正反取两次,这就是对这条reads的反向互补序列再取一次k-mer。下面的图就形象化了这一过程,长度为15的序列,选取k-mer为5,那么就会得到11(15-5+1=11)个5-mer。
如果你有兴趣,可以跳转进一步了解。
https://cloud.tencent.com/developer/article/1613847
2 利用kallisto定量
安装kallisto还是很轻松
conda install kallisto
定量需要两步,第一步是对你的数据建立index。之后就能用建立好的index做RNA-seq数据定量。
建立目录的命令很简单:
kallisto index -i kallisto_IWGSC transcripts.fasta
由于可变剪切等原因,同一个mRNA可能有不止一个isoform,如果你只在乎某个基因转录了多少,不在乎有多少同一个mRNA有多少个isoform的话,那么可以用TBtools提取每个mRNA的最独特的序列。一般来讲提取的是所有isoform中最长的那个序列。听起来很合理,但有些时候会出问题,比如那个最长序列本身不太对的时候。
定量之后只需要一个for语句循环就能完成RNA-seq的定量
#kallisto定量分析
for i in `seq 56 79`;
do
mkdir quants/SRR51316${i};
kallisto quant -i '/mnt/d/onecloud/OneDrive/wheat geno/wheat_index/kallisto_IWGSC' --single -o quants/SRR51316${i} -l 200 -s 30 SRR51316${i}.fastq.gz;
done
salmon定量分析
salmon的逻辑跟kallisto是一样的,都是先建立index再定量。
跟kallisto不同得地方是,salmon支持更多种的index模式
小麦因为基因组太大了,试了好几次建库都不完整,索性只说最简单的,既基于转录组数据建立index。如果你目标基因组比较小(如水稻,拟南芥),推荐你们尝试一下selective alignment。
地址为:
https://combine-lab.github.io/alevin-tutorial/2019/selective-alignment/
废话不多说,salmon的建立index得命令差不太多。就是把index和输出参数换了个位置
具体命令:
$ salmon index -t transcripts.fasta -i IWGSC_1.1
其中
-t 是转录组数据文件
-i 是输出地址
之后再进行定量就好
for i in `seq 56 79`;
do
mkdir quants/SRR51316${i};
salmon quant -i IWGSC_1.1 -l A \
-r /SRR51316${i}.fastq_clean.fq \
-p 8 --validateMappings --gcBias -o quants/SRR51316${i};
done
如果是双端测序数据
salmon quant -i index -l A\
-1 read1.fastq.gz -2 read2.fastq.gz \
-p 8 --validateMappings -o read_quants
其中
-i 对应的是index地址
-1和-2对应双端测序两个fastq文件
-p 是设定核心数量
-o 是输出地址
TBtools的众筹插件
这部分可能是最没必要讲的了。
都是最简单的东西了,需要注意的是,TBtools插件每次定量前都默认会重新建立一次index,所以……数据多的时候耗时会比较久……
不过这个插件最后会自动统计并整理好gene counts和TPM文件,方便后续DEseq2的操作。
最后alignment-free到底准不准
首先老版本好像有问题,新版本修复了。具体参考马省伟大佬的文章。
http://blog.sciencenet.cn/blog-1094241-1133526.html
其次,好像是不是alignment-free最终准确率都差不多。但是基于k-mer明显快很多,而且普通pc就能跑……具体参考发表在NC上的文章
https://www.nature.com/articles/s41467-017-00050-4
抱歉鸽了这么久……也不知道为啥会鸽……