我们常见的转录组表达分析一般都是将reads比对至参考基因组或者转录组上,然后在基因或者转录本水平上定量表达丰度。
但最近在做小RNA分析时却遇到了没有参考基因组注释文件(gtf/gff文件)的情况,而注释文件的缺失则意味传统的转录组定量分析是无法进行的。那在缺少注释文件的情况下,该如何进行定量分析呢?在各种搜索后发现了一款无需mapping便可进行定量的软件——Salmon。
一、基本情况
Salmon软件于2017年发表在Nature Methods,其题目为《Salmon provides fast and bias-aware quantification of transcript expression》
Salmon 提供2种运行模式,一是quasi-mapping直接读取 reads 文件;二是读取比对文件 sam/bam 进行mapping。
1、quasi-mapping-based mode的运行有两阶段:构建索引和用户想要定量的reads文件。
2、alignment-based mode的运行则不需要构建索引,而是仅需提供一个转录本的 FASTA文件和用户想要定量的 SAM/BAM 文件。
二、软件使用:
1、quasi-mapping-based mode
构建索引:
salmon index -t transcripts.fa -i transcripts_index -k 31
参数说明:
-t:转录本的fasta文件
-i:输出目录
-k:K-mers,默认值为31
#如果你的reads大于75bp,那么k设置为31是较好的选择,如果reads低于75可略微减少K值
名词解释:
简单来说,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。
定量分析:
#双端测序数据reads表达量的估计
salmon quant -i transcripts_index -l <LIBTYPE> -1 reads1.fq -2 reads2.fq -o transcripts_quant
#单端测序数据reads表达量的估计
salmon quant -i transcripts_index -l <LIBTYPE> -r reads.fq -o transcripts_quant
参数说明:
-1/2:双端数据
-r:单端数据
-l:--libType,测序文库类型,一般不知道什么文库的话用参数 A 让软件自动检测
#I = inward
#O = outward
#M = matching
#S = stranded
#U = unstranded
#F = read 1 (or single-end read) comes from the forward strand
#R = read 1 (or single-end read) comes from the reverse strand
#A = automatically determine
2、alignment-based mode
该模式下无需创建索引
salmon quant -t transcripts.fa -l <LIBTYPE> -a aln.bam -o salmon_quant
3、输出文件
主要输出文件为quant.sf
,该文件共有5列,分别是Name,Length ,EffectiveLength,TPM和NumReads。
- Name — target transcript 名称, 由输入的 transcript database (FASTA file)所提供。
- Length — target transcript 长度,即有多少个核苷酸
- EffectiveLength — target transcript 计算的有效长度。此项考虑了所有被建模的因素,这将影响从这个转录本中取样片段的概率,包括片段长度分布和序列特异性和gc片段偏差(如果这些因素在建模时均被考虑的话)。 (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))。
- TPM — 估计转录本的表达量。
- NumReads — 估计比对到每个转录本的reads数。
其他输出文件:
cmd_info.json
: JSON格式文件,记录salmon程序运行的命令和参数
lib_format_counts.json
: Observed library format counts。当运行salmon是 mapping-based mode时,则会生成改文件。 JSON格式文件,记录有关文库格式和reads比对的情况。
eq_classes.txt
: Equivalence class file。当Salmon运行时,应用参数--dumpEq,则会生成此文件。
aux_info
: 辅助文件夹,内含多个文件
fld.gz
:在辅助文件夹中,该文件记录的是观察到的片段长度分布的近似值
obs5_seq.gz, obs3_seq.gz, exp5_seq.gz, exp5_seq.gz
: Sequence-specific bias files
expected_gc.gz, observed_gc.gz
: 当Salmon运行时,应用fragment-GC bias correction,在辅助文件夹中则会生成这两个文件。记录Fragment-GC bias。
meta_info.json
: JSON格式文件,记录salmon程序运行的统计信息
ambig_info.tsv
: tab分隔符的文本文件,含有两列。记录的是每个转录本对应的 the number of uniquely-mapping reads 和 the total number of ambiguously-mapping reads
三、补充
TPM:
Transcripts Per Kilobase of exonmodel per Million mapped reads (每千个碱基的转录每百万映射读取的Transcripts),优化的RPKM计算方法,可以用于同一物种不同组织的比较。
TPM概括了基因的长度、表达量和基因数目。TPM可以用于同一物种不同组织间的比较,因为sum值总是唯一的。
计算公式:PMi=(Ni/Li)*1000000/sum(Ni/Li+……..+ Nm/Lm)
其中:Ni:mapping到基因i上的read数; Li:基因i的外显子长度的总和
http://blog.sciencenet.cn/blog-1113671-1038659.html
参考:
https://www.bioinfo-scrounger.com/archives/411/
Salmon 进行转录本定量https://www.jianshu.com/p/f62fd85113d3
tximport 将 Salmon 定量结果导入 DESeq2https://www.jianshu.com/p/e0acb957b351
salmon分析RNA-seq实战https://www.jianshu.com/p/5ffbe89d3b6b