Salmon 进行转录本定量

Salmon 是转录本定量软件，使用转录本定量主要优点一是更准确，比如不同样本同一基因使用不同 isoform 此时基因长度并不相等，直接基因定量无法顾及；二是方便进行可变剪切分析。Salmon 提供2种运行模式，一是直接读取 reads 文件(姑且称为 fastq 模式)；二是读取比对好文件 sam/bam (姑且称为比对模式)。

在文章《StringTie + DESeq2 进行 RNA-seq 差异基因分析》和《featureCounts 计算基因 reads 数目》我介绍了一个基本的 RNA-seq 差异基因分析流程，Salmon 能用来替代这里的 StringTie 和 featureCounts, DESeq2 也专门有导入 Salmon 数据的函数。这篇文章我用人种 RNA-seq 数据演示一遍 Salmon 流程，下篇文章将用 Salmon 数据进行可变剪切分析。

Reads 乱序
Salmon 要求 reads 是无序的，如果你的文件将 reads 排序过，用 Salmon 前要先重新乱序。

线程数
fastq 模式差不多使用线程越多运行越快，不过实际操作也没必要给整个几十个线程去跑；比对模式线程到达一定数目后，总体速度就不怎么加快了，设置 8-12 线程是相对合适的。

文库类型(library type, LIBTYPE)
RNA-seq 有多种多样文库类型，使用 Salmon 应设置好这个参数，因为 Salmon 自动检测也许会出现错误。下面图示是 Salmon 定义的文库类型分类。

Salmon 文库类型

所以 LIBTYPE 设定包含了三个部分，一是 reads 的相对方向，分别用 I/O/M 表示；二是文库有无链特异性，分别用 S/U 表示；三是，假如文库有链特异性，指明链方向，分别用 F/R 表示。如果第二部分是 U 就不需要设置这里。下面是这些字母代表的全称。

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

要注意 Salmon 用 F/R 表示链方向，其他软件可能喜欢用 F2R1/F1R2形式。

流程示例

第一步是建立索引，建议先用 generateDecoyTranscriptome.sh 脚本产生一个 decoys.txt 文件，下面是我命令示例。
注：generateDecoyTranscriptome.sh 脚本在 GitHub - COMBINE-lab/SalmonTools: Useful tools for working with Salmon output 下载，或者 Github 搜索 SalmonTools.

/Example/Project/SalmonTools/scripts/generateDecoyTranscriptome.sh -a gencode.v33.annotation.gtf -g bigZips/hg38.fasta -t gencode.v33.transcripts.fa -o hsa_decoy

这里 -a 和 -t 参数用到文件可以在 Gencode 网站下载(人种)，这个命令我跑了2天时间，所以别以为出错了强制中断，没报错就等待。结束后 -o 目录会输出 decoys.txt 文件。
完成后建立 index 命令：

salmon index -t gencode.v33.transcripts.fa -i hsa_transcripts_index --decoys hsa_decoy/decoys.txt -k 31 --gencode

设置 -k 31 适合 75 bp 及更长 reads，如果你的 reads 更短需要设置小点。或者你觉得比对率有些偏低，可以将 -k 设小些，提高灵敏度。建立索引完成后 -i 参数指定目录会有许多文件。

$ls hsa_transcripts_index
complete_ref_lens.bin  ctg_offsets.bin         eqtable.bin  mphf.bin  pre_indexing.log  refAccumLengths.bin  reflengths.bin  seq.bin
ctable.bin             duplicate_clusters.tsv  info.json    pos.bin   rank.bin          ref_indexing.log     refseq.bin      versionInfo.json

fastq 模式
注意 -l 参数放 -1/-2 前面。

salmon quant -l MU -i /share/database/openData/GRCh38_hg38/hsa_transcripts_index --validateMappings -1 5637_50nm_1_Clean_Data1.fq.gz -2 5637_50nm_1_Clean_Data2.fq.gz -o Test1

运行后有 quant.sf 文件就是转录本定量信息。

Name    Length  EffectiveLength TPM     NumReads
ENST00000456328.2       1657    1372.105        0.034568        1.000
ENST00000450305.2       632     348.480 0.000000        0.000
ENST00000488147.1       1351    1066.105        3.604749        81.025
ENST00000619216.1       68      69.000  0.000000        0.000
ENST00000473358.1       712     427.508 0.000000        0.000
ENST00000469289.1       535     255.958 0.000000        0.000
ENST00000607096.1       138     5.557   0.000000        0.000
ENST00000417324.1       1187    902.105 0.000000        0.000

定量后 Salmon 工作也完成了，后面就是 DESeq2, limma 等软件进行下游分析。

部分参数解释

• --validateMappings 启用选择性比对到转录本。能够提高比对的灵敏度和特异性从而提高了定量的准确度。
• --mimicBT2 和 --mimicStrictBT2 让比对模仿一些 Bowtie2 行为。
• --recoverOrphans 这个参数应与 --validateMappings 一起用。如果只有一端 reads 能比对或者两端 reads 比对不在同一个转录本，Salmon 会在比对上位置的上下游寻找另一端 reads 的比对。
• --hardFilter 这个参数仅在 --validateMappings 模式使用，关闭软过滤(soft filter)。
• --allowDovetail 默认丢弃“燕尾比对”(Dovetailing mapping)，开启这参数后不丢弃。“燕尾比对”见下面 Bowtie 的示意。
• -p/--threads 线程数，默认会使用所有线程，想要用少些线程就设置。
• --incompatPrior 与文库类型(library type)不兼容的比对先验概率。Salmon 默认设置了一个非常小但非0数值，如果一个插入片段(fragment)的唯一比对是这种非兼容比对，Salmon 会把这个比对计数。如果希望不计入，将这个参数设置数值 0.0 即可。
• --minScoreFraction 最低比对得分限制，低于此比对被丢弃。
• --rangeFactorizationBins 提高对“困难”转录本的定量准度。配合比对模式的 --useErrorModel 参数或者 fastq 模式的 --validateMappings 参数使用。
• --seqBias 从输入数据学习然后矫正测序偏差。
• --gcBias 进行 GC 偏差矫正。建议根据质控的 GC 偏差情况设置。
• -g/--geneMap 提供基因注释文件 Salmon 就额外生成 quant.genes.sf 文件给出基因的定量。
  Bowtie2 mates

[参考]
Overview – Salmon: Fast, accurate and bias-aware transcript quantification from RNA-seq data
Bowtie 2: fast and sensitive read alignment