使用salmon进行alignment-free的定量 编辑

Salmon可以快速从fastq快速得到基因表达

Patro R , Duggal G , Love M I , et al. Salmon provides fast and bias-aware quantification of transcript expression[J]. Nature Methods, 2017, 14(4):417-419.

salmon的下载安装 编辑

conda install -y salmon

salmon的使用 编辑

$ salmon
Usage:  salmon -h|--help or 
        salmon -v|--version or 
        salmon -c|--cite or 
        salmon [--no-version-check] <COMMAND> [-h | options]

Commands:
     index 构建salmon index
     quant 对一个样本进行定量
     alevin 单细胞分析
     swim  软件开发者命令
     quantmerge 将多个定量文件合并成单个文件

首先构建cDNA序列的索引,salmon index构建索引的参数:

-v [ --version ]       		 打印版本号
-h [ --help ]           	 打印帮助信息
-t [ --transcripts ] arg   转录本fasta文件
-k [ --kmerLen ] arg (=31) quasi索引的kmers大小 
-i [ --index ] arg         salmon index
-p [ --threads ] arg (=2)  计算bias特征时使用的线程数
--perfectHash              [quasi index only] 
  												 使用perfect hash(而非dense hash)构建index
  												 对内存需求小(尤其定量时), 
  												 但需要更长的时间来构建index
-d [ --decoys ] arg        将这些序列是为一些已知转录本的同源序列
--type arg (=quasi)        构建index的类型; 该版本下的参数只有 "quasi"

salmon quant包括两种定量模式 --- 使用原始FastQ格式的read和比对后的BAM/SAM格式read,如果设置-a [ --alignments ]参数,salmon quant执行基于比对的定量模式,不设置则执行原始read的定量模式。使用salmon quant --help-reads命令查看quasi-mapping-based模式的帮助文档,使用salmon quant --help-alignment命令查看alignment-based模式定量的帮助文档。 salmon Quasi-mapping模式的参数:

## mapping 输入选项:
  -l [ --libType ] arg                  描述文库类型的字符串,A代表自动搜索
  -i [ --index ] arg                    salmon index
  -r [ --unmatedReads ] arg             非mate文库,如单末端测序
  -1 [ --mates1 ] arg                   双末端1 mates
  -2 [ --mates2 ] arg                   双末端2 mates
  
## Quasi-mapping 模式选项,以下参数只支持Quasi-mapping模式:
  --discardOrphansQuasi                 [只支持Quasi-mapping模式] :quasi-mapping模式下,只有双末端的两个read都比对上的read对用于表达量评估。如果没有检测到双末端序列,则对单个read进行计数,使用 --writeOrphanLinksThis 可以将双末端只有一个read比对上的转录本计数结果输出到文件中。
  --validateMappings                    [Quasi-mapping模式] :使用比对来验证Quasi-mapping,保证结果合理。
  --consensusSlack arg (=0)             [Quasi-mapping模式] : The amount of slack allowed in the quasi-mapping  consensus mechanism.  一般情况下,一个转录本所有的hits都被认为是mapping,如果设置成分数X,区间在[0,1)之间,表示一个transcript有(100 * X)%的hits没有mapping,使用 --validateMappings 选项时,该选项的默认值是0.2,0表示所有hits都是mapping。
  --minScoreFraction arg                [Quasi-mapping模式 (w / mapping validation)] : 达到 "valid" 的最小比对分数,区间为(0,1]。Salmon默认0.65,Alevin模式下默认0.87。
  --maxMMPExtension arg (=7)            [Quasi-mapping模式 (w / mapping 
                                        validation) ] : 设置最大的MMP值,防止出现过大的MMSP(maximum mappable safe prefix),该值越小越敏感,但会降低mapping的速度。
  --ma arg (=2)                         [Quasi-mapping mode (w / mapping validation) ] : read比对到参考序列的match值。
  --mp arg (=-4)                        [Quasi-mapping mode (w / mapping validation) only] : read比对到参考序列的mis-match值。
  --go arg (=4)                         [Quasi-mapping mode (w / mapping validation) only] : 比对中允许gap数。
  --ge arg (=2)                         [Quasi-mapping mode (w / mapping 
                                        validation) only] : 比对中gap延伸数。
  --bandwidth arg (=15)                 [Quasi-mapping mode (w / mapping 
                                        validation) only] : 传递给ksw2的bandwidth的值,bandwidth值越小,比对验证的越快,也有可能会错过有效的比对
  --allowDovetail                       [Quasi-mapping mode] : 允许dovetailing mapping。
  --recoverOrphans                      [Quasi-mapping mode] : 尝试使用edlib修复mate文库的只比对上一个的orphan read,增加了计算时间,提高了敏感性。
  --mimicBT2                            [Quasi-mapping mode (w / mapping  validation)] : 类似于Bowtie2的 --no-discordant  --no-mixed 选项,不允许dovetailing reads并去除orphans read。和RSEM+Bowtie2的模式下的更严格筛选参数不同,如无gap比对参数,详见 --mimiStrictBT2 选项。
  --mimicStrictBT2                      [Quasi-mapping mode (w / mapping  validation)] : 设置用于RSEM+Bowtie2的严格过滤参数, --minScoreFraction会被提高到0.8,比对中不允许dovetailing read和gap,并去除未配对read。
  --hardFilter                          [Quasi-mapping模式 (w / mapping 
 validation)] : 不考虑次优的比对结果,默认是根据比对打分进行soft-filtering,如果想要得到更精确的丰度评估结果,而非一个区间值,需要设置此参数。
  --allowOrphansFMD                     [FMD-mapping模式] : 执行轻量级比对时,将未配对的read视为有效的hits。该选项能提高软件的敏感性(比对上更多的read、检测到更多的转录本),可能会降低准确性,未配对的比对可能是不正确的。
  -z [ --writeMappings ] [=arg(=-)]     将quasi-mapping的结果写入到SAM格式中,默认情况下,结果会直接输出到屏幕上,可以自定义输出文件名。
  -c [ --consistentHits ]               quasi-mapping下,强制将hits聚集(如共线性和) gathered during quasi-mapping to be "consistent" (i.e. co-linear and approximately the right distance apart).
  --fasterMapping                       [Developer]: 关闭quasi-mapping的额外检查。会加快mapping的速度,但可能会为某些read返回过多的映射(一些次优的mapping)。
  -x [ --quasiCoverage ] arg (=0)       [实验]: 被MMP(长度大于31)覆盖到的read分值,read “mapped”的阈值。帮助避免"spurious"比对。默认值0表述不设置覆盖阈值,精确match的覆盖度,较大MMP是一个严格条件,使用时设置为较小的值。

salmon alignment模式的参数:

## alignment输入选项:
--discardOrphans                      [Alignment-based模式] : 只对成对的比对进行表达量评估,默认:如果没有成对的mapping存在,则统计不成对的比对。
  -l [ --libType ] arg                  描述文库类型的字符串,A代表自动搜索
  -a [ --alignments ] arg               指定输入为BAM格式文件
  -t [ --targets ] arg                  用于定量的FASTA格式的转录本序列文件

## alignment模式选项,只支持alignment模式:
  --noErrorModel                        关闭alignment误差模型,该模型会统计比对中不同mismatch/indel类型的频率。启用该选项可加速alignment模式下salmon的运行,可能会降低定量的准确性。
  --numErrorBins arg (=6)               应用error模型时,将每个read分成的区间数。例如,10表示error模型将read分成10个区段进行高效处理,3表示将read分成read起始三分之一,中间三分之一,最后三分之一进行处理。
  -s [ --sampleOut ]                    根据统计的转录本丰度对比对结果进行抽样,并写入"postSample.bam"。该选项会过滤掉一些有歧义的片段。
  -u [ --sampleUnaligned ]              对比对的read进行额外的抽样,将未比对的read写入"postSample.bam"  --gencode                             指定输入格式为GENCODE,转录本名称通过第一个 '|' 字符提取,这些名称用于输出,以及在GTF文件的基因中,搜索这些转录本。
  --mappingCacheMemoryLimit arg (=2000000)
                                        如果文件中包含的read数低于此值,则会将文件中所有数据读入内存中用于计算。如果不想消耗过多内存,该值不要设置过大。对于只有几百万的read文件,可以将该值设置为文件read数,可以加快软件的运行速度。

salmon quant的选项,两种模式下的通用选项

# salmon quant选项:
## 基础选项:
  -v [ --version ]                      打印版本信息
  -h [ --help ]                         打印帮助信息
  -o [ --output ] arg                   输出结果文件夹
  --seqBias                             执行序列偏好性校正
  --gcBias                              执行fragment GC偏好校正(单末端read的测序版)
  -p [ --threads ] arg (=8)            	使用的线程数
  --incompatPrior arg (=0)              特定文库类型(--libType)下,比对结果与实际fragment来源不同的先验概率,设置为0表示比对结果与实际文库相同,不存在比对错误的情况,设置为1表示比对结果与文库类型完全不对应
  -g [ --geneMap ] arg                  基因和转录本的对应关系。使用该选项,salmon在结果中提供quant.sf和quant.genes.sf文件,quant.genes.sf文件中包含了gene水平的丰度评估结果。输入文件可以是GTF或制表符分隔的文件(每行包含一个transcript和gene的对应关系,以制表符分隔),文件的扩展名用于判断文件的处理方式。以'.gtf', '.gff' '.gff3'按照GTF格式处理;其他扩展名的文件都被当作制表符分隔文本文件处理。GTF/GFF格式中, "transcript_id" 包含transcript名,"gene_id"包含对应的gene名。
  --meta                                如果使用Salmon分析metagenomic数据集,可以考虑设置此参数,禁用对metagenomic数据意义不大的分析。

## 高级选项:
  --alternativeInitMode                 [实验]: 使用其他可选策略(不同于简单插入法)初始化EM/VBEM程序,用于online和归一化丰度的评估。
  --auxDir arg (=aux_info)              输出结果的子文件夹中输出辅助信息如 bootstraps, bias parameters等。
  --skipQuant                           跳过transcript定量(包括Gibbs抽样或bootstrap)  --dumpEq                              quasi-mapping下统计equivalence class的数目。
  -d [ --dumpEqWeights ]                在输出的文件中,包含更丰富的equivalence class 权重信息。
  --minAssignedFrags arg (=10)          用于转录组定量的最小fragment数。
  --reduceGCMemory                      降低用于计算片段GC含量的内存数,但会降低运行速度,结果不变。
  --biasSpeedSamp arg (=5)              评估序列特异和GC fragment bias时,fragment长度PMF下采样的值。值越大,有效长度校正速度越快,但会降低bias模型的可靠性。
  --fldMax arg (=1000)                  用于计算经验分布的最大fragment长度。
  --fldMean arg (=250)                  用于计算fragment长度分布的平均值
  --fldSD arg (=25)                     fragment长度分布的标准偏差
  -f [ --forgettingFactor ] arg (=0.65000000000000002)
                                        online学习策略下遗忘因子。值越小,学习速度越快,但不稳定,会产生更高的方差。较大的值会减慢学习的速度,但结果更稳定,该值的区间为(0.5, 1.0]  --initUniform                         offline推断前,初始化使用的归一化参数。
  -w [ --maxReadOcc ] arg (=200)        比对上的区域数超过此值的read不用于计算。
  --noLengthCorrection                  [experimental] : 在进行转录本丰度评估时,禁用长度校正,该选项用于不需要考虑目标片段长度的测序技术,如QuantSeq。
  --noEffectiveLengthCorrection         在计算片段是否来自转录本时,关闭有效长度校正,设置该参数后,不在考虑fragment长度分布。
  --noFragLengthDist                    [experimental] : 在判断某个fragment是否来源于特定位置时,不考虑fragment的长度分布。fragment长度观测值对fragment的鲜艳概率没有影响。
  --noBiasLengthThreshold               [experimental] : 不设置bias校正有效长度的最小阈值,会提高bias校正的正确性,但会降低鲁棒性。默认会执行有阈值的校正。
  --numBiasSamples arg (=2000000)       构建序列特异bias模型时使用的fragment数。
  --numAuxModelSamples arg (=5000000)   用于auxiliary模型参数训练的样本数(如片段长度分布,bias等),第一次训练后,这些模型参数收敛并固定。
  --numPreAuxModelSamples arg (=5000)   在不使用上个选项提供的auxiliary模型参数的情况下,评估分配的可能性和对转录本丰度计算的贡献。该选项的目的是:在确定auxiliary模型参数被充分训练之前,避免使用该模型。
  --useEM                               使用常规EM程序优化批处理过程。
  --useVBOpt                            使用Variational Bayesian EM[默认]
  --rangeFactorizationBins arg (=0)     基于不同转录本fragment的条件概率,使用新的equivalence class,将定量过程中的可能性分解为区间。默认值0对应标准富集,值越大,分解的区间会更多更细。如果启用此参数,该参数值为4。
  --numGibbsSamples arg (=0)            Gibbs抽样执行的次数。
  --noGammaDraw                         Gibbs抽样时,该参数会关闭的Gamma分布的统计transcript分数的操作。取样不考虑shot-noise,只考虑assignment ambiguity。
  --numBootstraps arg (=0)              bootstrap样本数,与Gibbs抽样互斥。 
  --bootstrapReproject                  从每个样本的bootstrap计数中学习参数分布,但会将这些学习的参数重新用原始equivalence class计数。
  --thinningFactor arg (=16)            Gibbs链中每个样本丢弃的步骤数,该数字越大,后续样本相关的计划越小,但会减慢取样速度。
  -q [ --quiet ]                        不打印提示信息,除非软件运行出错
  --perTranscriptPrior                  transcript水平的prior (使用默认或--vbPrior,每个transcript都会有prior值的read)
  --sigDigits arg (=3)                  输出EffectiveLength和NumReads列时要写入的有效位数
  --vbPrior arg (=1.0000000000000001e-05)
  																			VBEM程序中使用的prior值,用于per-nucleotide prior,如果设置--perTranscriptPrior选项,该值则用于transcript水平 prior值
  --writeOrphanLinks                    输出双末端中只有一个read支持的transcript。
  --writeUnmappedNames                  在auxiliary文件中输出unmapped_names.txt,其中是未map的read的统计信息。

使用salmon alevin模式进行单细胞定量分析 编辑

Alevin - Salmon 1.5.2 documentation

alevin选项:

mapping输入参数:
  -l [ --libType ] arg                  指定文库类型
  -i [ --index ] arg                    salmon index
  -r [ --unmatedReads ] arg             非mate文库read(单末端read)
  -1 [ --mates1 ] arg                   包含 #1 mates 的文件
  -2 [ --mates2 ] arg                   包括 #2 mates 的文件
alevin-specific Options:
  -v [ --version ]                      打印版本信息
  -h [ --help ]                         打印帮助信息
  -o [ --output ] arg                   指定定量输出结果储存的文件夹
  -p [ --threads ] arg (=2)             使用的线程数
  --tgMap arg                           tsv格式的转录本和基因的对应关系
  --hash arg                            Alevin使用Big freaking Hash (bfh.txt)文件作为二级输入。必须与 --chromium选项一起使用
  --dropseq                             处理DropSeq单细胞文库
  --chromiumV3                         	处理10x chromium v3单细胞文库
  --chromium                            处理10x chromium v2单细胞文库
  --gemcode                             处理0x gemcode v1单细胞文库
  --celseq                              处理CEL-Seq单细胞文库
  --celseq2                             处理CEL-Seq2单细胞文库
  --whitelist arg                       包含white-list barcodes的文件
  --noQuant                             不运行下游的barcode-salmon模型
  --numCellBootstraps arg (=0)          计算cell x基因矩阵定量的平均值和方差
  --forceCells arg (=0)                 指定细胞的数目
  --expectCells arg (=0)                定义细胞的期望数目
  --mrna arg                            线粒体RNA基因文件的路径,每行一个
  --rrna arg                            核糖体RNA文件的路径,每个一行
  --keepCBFraction arg (=0)             CB保留的分值,范围在(0,1], 
  																			 1 定量所有CB
  --dumpfq                              多比对使用coin toss,为下游分析,
  																			保留barcode修正后的fastq 文件     
	--dumpBfh                             保留所有barcode和UMI序列的big hash
  --dumpUmiGraph                        保留细胞水平的UmiGraph
  --dumpFeatures                        保留whitelist和下游分析
  --dumpMtx                             以mtx格式保存cell v转录本计数矩阵
  --lowRegionMinNumBarcodes arg (=200)  CB最小值作为学习的最低执行区域 
                                        (Default: 200)
  --maxNumBarcodes arg (=100000)        处理细胞barcode的最大值
  																			(Default: 100000)

使用salmon quantmerge模式对多个定量文件进行合并 编辑

salmon quantmerge 选项:

基础选项:
  -v [ --version ]            打印版本号
  -h [ --help ]               打印帮助信息
  --quants arg                列出定量文件夹
  --names arg                 可选的样本名
  -c [ --column ] arg (=TPM)  需要合并的列名,
  														可选参数:{len, elen, tpm, numreads}
  --genes                     使用基因定量代替转录本定量
  --missing arg (=NA)         NA代替缺失值The value of missing values.
  -o [ --output ] arg         输出定量文件

salmon运行示例 编辑

salmon quant -i hg38_salmon -l A -1 ${i}_1_val_1.fq.gz -2 ${i}_2_val_2.fq.gz -p 5 -o ${i}_quant

salmon结果解析 编辑

quant.sf结果文件:

Name    Length  EffectiveLength TPM     NumReads
ENST00000631435.1       12      13.000  0.000000        0.000
ENST00000434970.2       9       10.000  0.000000        0.000
ENST00000448914.1       13      14.000  0.000000        0.000
ENST00000415118.1       8       9.000   0.000000        0.000
ENST00000604446.1       23      24.000  0.000000        0.000

各列含义解析:

  • Name:target transcript 名称,由输入的 transcript database (FASTA file)所提供。
  • Length:target transcript 长度,即有多少个核苷酸。
  • EffectiveLength:target transcript 计算的有效长度。此项考虑了所有建模的因素,这将影响从这个转录本中取样片段的概率,包括片段长度分布和序列特异性和gc片段偏好
  • TPM:估计转录本的表达量
  • NumReads:估计比对到每个转录本的reads数。

Salmon输出其他文件:

  • 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:在辅助文件夹中,该文件记录的是观察到的片段长度分布的近似值
    • observed_bias_3p.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分隔符的文本文件,含有两列。记录的是每个转录本对应的唯一mapping的read数和模糊mapping的read总数