一、 基本介绍

Salmon是一款新的、极快的转录组计数软件。它与Kallisto和Sailfish类似，可以不通过mapping而获得基因的counts值。Salmon的结果可由edgeR/DESeq2等进行counts值的下游分析。

Salmon是一种准确快速定量转录本丰度的方法；它是第一个可校正转录组片段GC含量范围内偏差的定量工具，大大提高了丰度估计的准确性以及后续差异表达分析的可靠性；其计算速度快，需要的计算资源小，硬盘空间占用少，不会生成巨型SAM格式比对文件。该软件2015年发布在BioRxiv上面，两年只有10来次的引用，4年也只有37个引用。而文章2017年被Nature Methods接收发表后，短短两年引用达891次，相同时间内引用高达上百倍。

我们常见的转录组表达分析一般都是将reads比对至参考基因组或者转录组上，然后在基因或者转录本水平上定量表达丰度。如果是在基因水平上的定量分析，那么一般可以选择用HISAT2比对，HTSeq-count定量；如果是在转录本水平，可以先用HISAT2比对，然后Stringtie组装转录本（假如只对已知isoforms定量，那么这步一般是可以省略的），再用RSEM或者eXpress定量表达丰度。其实以上步骤在转录本定量中属于同一类：Alignment-based transcript quantification，还有一类则是：Alignment-free transcript quantification；后者相比前者而言，省略了比对这一步骤（有些软件是只进行轻度比对），从而直接将reads分配到转录本上；后者相比前者spliced alignment能极大的减少计算机资源的使用，代表软件有Sailfish，Salmon，kallisto。

Salmon使用了expressive and realistic模型以及结合experimental attributes and biases使得其推测的表达量更加符合实际的RNA-seq数据。Salmon原理如下图所示，概括来讲是通过构建统计模型来推测已经注释的转录本呈现出什么表达模式时我们才会测序产生当前的FASTQ数据：

现在Salmon支持两种模式将reads mapping到转录本上：第一种是quasi-mapping-based mode，是一种最新的以及快速准确的定量转录本的方法，不需要像传统的那样需要通过比对才能将reads mapping到转录本上；第二种则是alignment-based mode，跟类似RSEM软件一样，提供比对后的bam/sam文件即可对转录本进行定量。第一种方法需要对转录本建index，第二种方法则不需要。

—— Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., & Kingsford, C. (2017). Salmon provides fast and bias-aware quantification of transcript expression. Nature Methods.

二、 使用方法

(1) 构建索引

salmon index  -t cdna.all.fa  -i salmon_index  --type quasi  -k 31

-t：参考转录组。
-i：在salmon_index/文件夹下得到索引。
-k：是指k-mers长度，默认是31，减少其值能略微增加mapping的灵敏度。如果reads大于75bp，那么k设置为31是较好的选择，如果reads低于75可略微减少k值。
-type：索引类型，分为fmd、quasi，建议quasi。新版只能用puff。

(2) 构建decoy-aware transcriptome file

在0.14.0 新版本中，salmon引入了新的索引方式，且旧版本产生的索引不能在0.14.0版本后通用。一种方法是使用generateDecoyTranscriptome.sh脚本。另一种方法是用生物体的整个基因组作为诱饵序列，这可以通过将基因组连接到你想索引的转录组的末尾，并用染色体名称填充decoys.txt文件来实现。第二个方案提供了一组更全面的诱饵，但是，显然需要更多的内存来构建索引。

• 方法一：SalmonTools的generateDecoyTranscriptome.sh（https://github.com/COMBINE-lab/SalmonTools/archive/master.zip）。这是一个预处理脚本，用于为salmon index创建增强混合fasta文件。它需要 -g genome fasta, -t transcriptome fasta, -a the annotation GTF file 来创建一个新的hybrid fasta文件，包含decoy sequences连接到transcriptome (gentrome.fa)。它转储的decoys.txt文件，其中包含decoy sequences的名称/id。

generateDecoyTranscriptome.sh  -a annotation.gtf  -g hg38.fasta  -t cdna.all.fa  -o decoy

-o：输出目录decoy/，在这个目录下获得decoys.txt和gentrome.fa两个文件。

salmon index  -t gentrome.fa  -i salmon_index  --decoys decoys.txt  -k 31

-i：在salmon_index/文件夹下得到新版本的索引。

• 方法二：用生物体的整个基因组作为诱饵序列，构建一组更全面的诱饵。

# 去掉空格后面的字符串，保证cDNA文件中fasta序列的名字简洁，不然后续会出错。
cut  -f  1  -d  ' '  cdna.all.fa  >  cdna.new.fa
# 获取所有序列的名字存储于decoy中。
grep  '^>'  cdna.new.fa  |  cut  -d  ' '  -f  1  |  sed  's/^>//g'  >  cdna.decoys.txt
# 合并cDNA和基因组序列一起。注意cDNA在前，基因组在后。
cat  cdna.new.fa  GRCh38.fa  >  trans_genome.fa
# 构建索引。更慢，结果会更准。
salmon index  -t trans_genome.fa  -d cdna.decoys.txt  -i salmon_index

(3) 定量

# quasi-mapping-based mode
salmon quant  -i salmon_index  -p 6  --libType IU  -1 reads.clean_1.fastq  -2 reads.clean_2.fastq  -o salmon_quant  #双端测序
salmon quant  -i salmon_index  -p 6  --libType U  -r reads.clean.fastq  -o salmon_quant  #单端测序
salmon quant  -i index  -l IU  -1 lib_1_1.fq lib_2_1.fq  -2 lib_1_2.fq lib_2_2.fq  -o out  #将多个样本放在一起作为一个库定量
salmon quant  -i index  -l IU  -1 <(cat lib_1_1.fq lib_2_1.fq)  -2 <(cat lib_1_2.fq lib_2_2.fq)  -o out  #将多个样本放在一起作为一个库定量

-o：输出文件夹。其中最主要的文件是quant.sf文件，其记录了每个转录本的count数以及对应的TPM值。是一个纯文本文件，用制表符分隔。

--libType是文库类型。A自动检测。如果是mRNA双端测序普通建库，那么选择的就是IU。单端测序普通建库为U。-l 参数放 -1/-2 前面。

--validateMappings 启用选择性比对到转录本。能够提高比对的灵敏度和特异性从而提高了定量的准确度。

# alignment-based mode
salmon quant -t transcripts.fa -l <LIBTYPE> -a aln.bam -o salmon_quant

需要输入一个aln.bam文件，然后直接与转录组进行比对，而不是与基因组比对。

(4) 下游tximport包

使用tximport包，可以导入salmon的转录本水平定量，并可选择将它们整合到基因水平，用于基因水平差异表达分析。Salmon生成的伪计数表示为标准化TPM计数，并映射到转录本。需要将这些转换为非正态化的计数估计，以执行DESeq2分析。要使用DESeq2，还需要将我们的丰度估计值从转录水平转换到基因水平。可以使用R Bioconductor软件包tximport完成上述所有操作。tximport可以纠正不同样本基因长度的潜在改变，比如differential isoform usage；能够用于导入Salmon、Sailfish、kallisto程序的输出文件；能够避免丢弃那些比对到多个基因的同源序列，从而提高灵敏度。

# 假设salmon输出文件为Salmon/quants/${sample}/，里面含有quant.sf文件。
sampleList <- c("sample1", "sample2", " sample3"…)
fileList <- file.path("/…/Salmon/quants", sampleList, "quant.sf")
names(fileList) <- sampleList
# 还要准备一个基因名和转录本名字相关的数据框，用于tximport的tx2gene必须第一列是转录本ID第二列是基因ID，对列名无要求。
txd <- makeTxDbFromGFF(file = "/Example/GRCh38_hg38/gencode.v33.annotation.gtf", format = "gtf", dataSource = "gencode.v33.annotation.gtf", organism = "Homo sapiens")
txTogene <- AnnotationDbi::select(txd, keys=keys(txd, "TXNAME"), keytype="TXNAME", columns=c("TXNAME", "GENEID"))
# tximport整合数据
txi <- tximport(fileList, type = "salmon", tx2gene = txTogene, countsFromAbundance="lengthScaledTPM")

尽管我们的quantum.sf文件中有一列对应于每个转录本的估计计数值，但这些值与有效长度相关。我们要做的是使用countsFromAbundance = ”lengthScaledTPM”参数。这将使用TPM列，并计算与原始计数相同尺度的数量，但不再与样本中的转录本长度相关。

files: 转录本水平丰度的文件名字符串向量。

type: 产生丰度所使用的软件名。包括"salmon", "sailfish", "alevin", "kallisto", "rsem", "stringtie", or "none"。该参数用于自动填写下面的参数（geneIdCol等）。"none "表示用户将指定这些列。

tx2gene: 一个两列的数据框，联系转录本id（第1列）和基因id（第2列）。列名可以无所谓，但一定要先是转录本再是基因的。这个参数是基因水平汇总所需要的。

countsFromAbundance: 字符串，可以是 "no"（默认）、"scaledTPM"、"lengthScaledTPM "或 "dtuScaledTPM"。是否使用丰度估计生成估计的计数。
a. 缩放为文库大小（scaledTPM）。
b. 使用样本上的平均转录本长度和文库大小（lengthScaledTPM）进行缩放。
c. 使用基因同种异构体之间的转录长度中位数，再加上文库大小（dtuScaledTPM）进行缩放。

ignoreTxVersion = TRUE：从Ensembl获取的参考转录组文件的标识符中包含版本号，而注释数据库通常不包含版本号，这将与tx2gene文件产生差异。

# 使用DESeqDataSetFromTximport将数据导入DESeq2
dds <- DESeqDataSetFromTximport(txi, colData = sampleGroup, design = ~ Group)

三、 举个栗子

(1) decoy

./SalmonTools-master/scripts/generateDecoyTranscriptome.sh  -a Drosophila_melanogaster.BDGP6.28.102.gtf  -g Drosophila_melanogaster.BDGP6.28.dna.toplevel.fa  -t Drosophila_melanogaster.BDGP6.28.cdna.all.fa  -o decoy
#输入-a注释GTF、-g基因组fasta、-t转录组fasta
#得到decoy/decoys.txt和decoy/gentrome.fa两个文件

(2) index

salmon index  -t ./decoy/gentrome.fa  -i salmon_index  --decoys ./decoy/decoys.txt  -k 31
#在salmon_index/文件夹下得到索引

(3) quant

for i in `ls cleandata/*.clean.fastq.gz`
do

    i=`basename $i`
    i=${i%.clean.fastq.gz}
    echo $i >> salmon_log.txt
    
    salmon quant  -i salmon/salmon_index  -p 4  --libType U  -r cleandata/$i\.clean.fastq.gz  -o salmon_quant/$i/  2>> salmon_log.txt

done

(4) tximport

# 需要整合的文件列表
sampleList <- dir("salmon_quant/")
fileList <- file.path("salmon_quant",sampleList,"quant.sf")
names(fileList) <- sampleList

# 自制一个TxDb
txd <- makeTxDbFromGFF(file="../salmon/Drosophila_melanogaster.BDGP6.28.102.gtf", format = "gtf")

# 转录本到基因的对应关系
txTogene <- AnnotationDbi::select(txd, 
                                  keys=keys(txd, "TXNAME"), 
                                  keytype="TXNAME", 
                                  columns=c("TXNAME", "GENEID"))

# 整合salmon数据
txi <- tximport(fileList, type = "salmon", 
                tx2gene = txTogene, 
                countsFromAbundance="lengthScaledTPM")
names(txi)
head(txi$counts)

# 导入DESeq2
#dds <- DESeqDataSetFromTximport(txi, colData = sampleGroup, design = ~ Group)

# 保存数据
save(txd, txi, file = "txd_txi_data.RData")