快速计算基因表达软件：Salmon

2021-01-31 本文已影响0人 BINBINCC

我们常见的转录组表达分析一般都是将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