RNA-seq:Salmon 快速定量
1RNA-seq:Salmon 快速定量
Salmon 原译鲑鱼,又称三文鱼,是深海鱼类的一种,也是一种非常有名的溯河洄游鱼类,通常认为鲑鱼就是靠味觉从海洋重新游回到它们原来的出生地,这往往要逆流而上经过很多的障碍。 ——百度百科
学习一件新的东西很难,若是没有兴趣,从 0 到 1 的可能性微乎其微。 —— Jun 言 Jun 语
2017 年 3 月 6 日在 Nature Methods 上发表了一篇转录本快速定量的 salmon 软件,影响因子 30.823 分,到现在引用次数达 2846 次!
2软件介绍
特点:
1、快速、准确定量转录本丰度。 2、自动矫正 GC 含量的偏差。 3、低消耗、低内存使用。 4、结合了一种新的双阶段并行推理算法和特征丰富的偏差模型与一个超快速的读取映射程序。 5、直接对 fastq 数据对转录本定量,不需要经比对过程。
软件一直在更新,从 0.5.1 版本更新到现在的 1.5.0 版本,在 2020 年 9 月 7 号发表在 Genome Biology 上推出了一种新的算法 selective alignment (SA),是对之前 1.0.0 所有版本的重大改进,进一步提高了定量的准确性。
3软件安装
1、下载编译好的最新二进制版:https://github.com/COMBINE-lab/salmon/releases
2、conda 安装最新版:conda install -c bioconda salmon
# 添加镜像源
$ conda config --add channels conda-forge
$ conda config --add channels bioconda
# 安装salmon
$ conda install -c bioconda salmon
比对模式
A、mapping-based mode :基于对转录组序列建立索引再进行定量两个步骤,有 quasi-mapping 和 Selective Alignment 两种类型
B、alignment-based mode :不需要建立索引,只需要提供转录本序列和比对好的 bam 或 sam 文件即可定量
第一种比对模式又有三种不同的类型:
1、cDNA-only index :又叫 salmon_index,属于 quasi-mapping 类型,这种方法将产生最小的索引,并需要最少的资源来构建,但是最容易出现虚假的比对结果。 2、SA mashmap index :又叫 salmon_partial_sa_index (regions of genome that have high sequence similarity to the transcriptome) ,属于 Selective Alignment 类型,含有基因组与转录组上高度相似性的序列的索引,为 partial decoy
索引,提高比对准确性。3、SAF genome index :又叫 salmon_sa_index (the full genome is used as decoy) ,利用整个基因组的信息来进一步的提高比对的准确性,减少假比对结果,属于 Selective Alignment 类型,索引文件也是三个里最大的,准确效果是最好的,为 full decoy
索引。作者推荐 Selective Alignment 的类型 :
decoy sequence : 基因组上未被注释的区域与转录组上的某些注释的区域含有较高的序列相似性,这些序列就为诱饵序列,salmon 比对产生的假比对就是这些序列造成的。
查看帮助文档
index 建立索引、quant 定量、alevin 单细胞分析和 quantmerge 合并定量结果等用法。
$ salmon
salmon v1.5.0
Usage: salmon -h|--help or
salmon -v|--version or
salmon -c|--cite or
salmon [--no-version-check] <COMMAND> [-h | options]
Commands:
index : create a salmon index
quant : quantify a sample
alevin : single cell analysis
swim : perform super-secret operation
quantmerge : merge multiple quantifications into a single file
作者还提供了已经构建好的索引下载:Pre-built versions of both the partial decoy
and full decoy
(i.e. using the whole genome) salmon indices for some common organisms are available via refgenie here:http://refgenomes.databio.org/
4软件使用实战
首先是工作目录文件内容:
/mnt/d/rnaseq/salmon$ ls -l |awk '{print $9}'
1.fastq-data/
GRCm39.primary_assembly.genome.fa
decoy_file/
decoys.txt
decoys.txt.bak
gencode.vM27.annotation.gtf
gencode.vM27.transcripts.fa.gz
generateDecoyTranscriptome.sh
gentrome.fa.gz
quant.sh
quants/
quants_fd/
quants_merge.sh
salmon-1.0.0_linux_x86_64.tar.gz
salmon-latest_linux_x86_64/
salmon_index/
salmon_index_full_decoy/
salmon_index_partial_decoy/
A、mapping-based mode ——cDNA-only index
1、建索引
现在小鼠的基因组和注释文件都更新到了 GRCm39,那我们就用这个吧。
# gencode下载转录本序列
/mnt/d/rnaseq/salmon$ wget http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M27/gencode.vM27.transcripts.fa.gz
# 建立索引,-p线程数量,-k为k-mer的数量,--gencode为转录本格式为gencode
# -t转录本序列,-i为输出所有文件夹名称,自动创建
/mnt/d/rnaseq/salmon$ salmon index -p 10 -k 31 --gencode \
-t gencode.vM27.transcripts.fa.gz \
-i salmon_index
# 只用了几分钟就建好了,查看产生了哪些文件
/mnt/d/rnaseq/salmon$ ls -lh salmon_index/
total 629M
-rwxrwxrwx 1 root root 550K Jun 14 16:28 complete_ref_lens.bin
-rwxrwxrwx 1 root root 36M Jun 14 16:30 ctable.bin
-rwxrwxrwx 1 root root 2.3M Jun 14 16:30 ctg_offsets.bin
-rwxrwxrwx 1 root root 66K Jun 14 16:28 duplicate_clusters.tsv
-rwxrwxrwx 1 root root 1003 Jun 14 16:30 info.json
-rwxrwxrwx 1 root root 76M Jun 14 16:30 mphf.bin
-rwxrwxrwx 1 root root 403M Jun 14 16:30 pos.bin
-rwxrwxrwx 1 root root 496 Jun 14 16:30 pre_indexing.log
-rwxrwxrwx 1 root root 18M Jun 14 16:30 rank.bin
-rwxrwxrwx 1 root root 1.1M Jun 14 16:30 refAccumLengths.bin
-rwxrwxrwx 1 root root 17K Jun 14 16:30 ref_indexing.log
-rwxrwxrwx 1 root root 550K Jun 14 16:30 reflengths.bin
-rwxrwxrwx 1 root root 58M Jun 14 16:30 refseq.bin
-rwxrwxrwx 1 root root 35M Jun 14 16:30 seq.bin
-rwxrwxrwx 1 root root 126 Jun 14 16:30 versionInfo.json
-k 参数为 k-mer 的数量,为 31 时适合 75bp 或者更长的 reads,若小于 75bp 可以尝试小一点的 k-mer 来提高灵敏度。
2、定量
数据还是我们上一期的文章里的 rnaseq 数据:RNA-seq : Hisat2+Stringtie+DESeq2
# 查看测序数据
/mnt/d/rnaseq/salmon$ ls 1.fastq-data/
test1_R1_rep1.fq.gz test1_R2_rep1.fq.gz test2_R1_rep1.fq.gz test2_R2_rep1.fq.gz
test1_R1_rep2.fq.gz test1_R2_rep2.fq.gz test2_R1_rep2.fq.gz test2_R2_rep2.fq.gz
# 建立输出文件夹
/mnt/d/rnaseq/salmon$ mkdir quants
# 单样本定量
# -i为索引文件夹,-l为测序文库类型,A为自动检测,-1、-2跟双端测序文件,单端测序为-r参数跟上测序文件
# 注意:-l参数必须在-1、-2或-r参数前面,因为文库类型标志用于确定应该如何解析reads的内容
# -p为线程数量,-o为输出文件夹
# --validateMappings打开selective alignment当比对到转录组上时,增加比对准确性和精确性
/mnt/d/rnaseq/salmon$ salmon quant -i salmon_index -l A \
-1 1.fastq-data/test1_R1_rep1.fq.gz \
-2 1.fastq-data/test1_R2_rep1.fq.gz \
-p 10 --validateMappings -o quants/test1_rep1_quant
# 查看产生了哪些文件,salmon_quant.log为日志文件,包含了比对率,quant.sf为定量结果
/mnt/d/rnaseq/salmon$ tree quants/
quants/
└── test1_rep1_quant
├── aux_info
│ ├── ambig_info.tsv
│ ├── expected_bias.gz
│ ├── fld.gz
│ ├── meta_info.json
│ ├── observed_bias.gz
│ └── observed_bias_3p.gz
├── cmd_info.json
├── libParams
│ └── flenDist.txt
├── lib_format_counts.json
├── logs
│ └── salmon_quant.log
└── quant.sf
# 查看定量结果文件,第1列为转录本id,2、3列为转录本长度和估算的有效长度
# 4、5列为tpm值和估计的来自改转录本的reads数量
/mnt/d/rnaseq/salmon$ less -S quants/test1_rep1_quant/quant.sf
Name Length EffectiveLength TPM NumReads
ENSMUST00000193812.2 1070 770.293 0.000000 0.000
ENSMUST00000082908.3 110 19.816 0.000000 0.000
ENSMUST00000162897.2 4153 3853.293 0.000000 0.000
比对率为 87.3207% :
批量定量:
# 创建一个脚本文件然后输入批量定量的脚本
/mnt/d/rnaseq/salmon$ vi quant.sh
#!/bin/bash
# batch salmon quant
for i in test1 test2
do
for j in rep1 rep2
do
echo "Processing sample ${i}_${j}"
salmon quant -i salmon_index -l A \
-1 1.fastq-data/${i}_R1_${j}.fq.gz \
-2 1.fastq-data/${i}_R2_${j}.fq.gz \
-p 10 --validateMappings -o quants/${i}_${j}_quant
done
done
# 保存后挂后台运行
/mnt/d/rnaseq/salmon$ nohup ./quant.sh &
大概十分钟左右就结束了。
如果定量时只想对
基因定量
,就需要把转录本合并为基因,这个参数可以设置:-g [ --geneMap ]
,定量后会产生 quant.sf 和 quant.genes.sf 文件,后者就是基因的表达定量结果,参数跟上 gtf、gff、gff3 后缀结尾的注释文件。
3、合并输出定量结果
我们需要使用 quantmerge 功能,首先查看一下帮助文档
/mnt/d/rnaseq/salmon$ salmon quantmerge -h
Version Server Response: Not Found
quantmerge
==========
Merge multiple quantification results into
a single file.
salmon quantmerge options:
basic options:
-v [ --version ] print version string
-h [ --help ] produce help message
--quants arg List of quantification directories.
--names arg Optional list of names to give to the samples.
-c [ --column ] arg (=TPM) The name of the column that will be merged
together into the output files. The options are
{len, elen, tpm, numreads}
--genes Use gene quantification instead of transcript.
--missing arg (=NA) The value of missing values.
-o [ --output ] arg Output quantification file.
1、-v [ --version ]
# 输出版本信息2、-h [ --help ]
# 输出帮助信息3、--quants arg
# 定量结果文件夹的 list,花括号包围4、--names arg
# 合并结果的样本名称的 list,花括号包围5、-c [ --column ] arg (=TPM)
# 指定哪一列合并,len、elen、tpm 或者 numreads,-c=tpm 格式6、--genes
# 使用基因定量代替转录本定量,首先要有进行基因的定量7、--missing arg (=NA)
# 缺失值8、-o [ --output ] arg
# 输出文件
开始合并:
/mnt/d/rnaseq/salmon$ cd quants/
# 合并输出tpm
/mnt/d/rnaseq/salmon/quants$ salmon quantmerge \
--quants {test1_rep1_quant,test1_rep2_quant,test2_rep1_quant,test2_rep2_quant} \
--names {test1_rep1,test1_rep2,test2_rep1,test2_rep2} \
--column=tpm \
-o ./merge_tpm.txt
Version Server Response: Not Found
[2021-06-15 09:35:42.246] [mergeLog] [info] samples: [ test1_rep1_quant, test1_rep2_quant, test2_rep1_quant, test2_rep2_quant ]
[2021-06-15 09:35:42.246] [mergeLog] [info] sample names : [ test1_rep1, test1_rep2, test2_rep1, test2_rep2 ]
[2021-06-15 09:35:42.246] [mergeLog] [info] output column : TPM
[2021-06-15 09:35:42.246] [mergeLog] [info] output file : ./merge_tpm.txt
[2021-06-15 09:35:42.252] [mergeLog] [info] Parsing test1_rep1_quant/quant.sf
[2021-06-15 09:35:43.161] [mergeLog] [info] Parsing test1_rep2_quant/quant.sf
[2021-06-15 09:35:44.008] [mergeLog] [info] Parsing test2_rep1_quant/quant.sf
[2021-06-15 09:35:44.805] [mergeLog] [info] Parsing test2_rep2_quant/quant.sf
# 查看产生文件
/mnt/d/rnaseq/salmon/quants$ head -4 merge_tpm.txt
Name test1_rep1 test1_rep2 test2_rep1 test2_rep2
ENSMUST00000082422.1 0 0 0 0
ENSMUST00000082421.1 6685.15 6750.1 5357.29 6426.71
ENSMUST00000082420.1 20.5699 0 0 0
# 同样的提取 NumReads
/mnt/d/rnaseq/salmon/quants$ salmon quantmerge \
--quants {test1_rep1_quant,test1_rep2_quant,test2_rep1_quant,test2_rep2_quant} \
--names {test1_rep1,test1_rep2,test2_rep1,test2_rep2} \
--column=numreads \
-o ./merge_reads.txt
A、mapping-based mode——SA mashmap index
这个模式 salmon 将基因组的外显子做了 mask 处理,就是用 N 代替,然后利用 mashmap 比对到转录组上获取最低相似度为 80%的基因间隔区的同源序列,得到的称为 decoy sequence,然后把这些序列和转录组序列合并在一起建立索引。
1、准备文件
需要可执行的 mashmap 和 bedtools 软件 salmon 提供了一个 shell 脚本来帮我们做这些事:https://raw.githubusercontent.com/COMBINE-lab/SalmonTools/master/scripts/generateDecoyTranscriptome.sh
# 下载脚本文件
/mnt/d/rnaseq/salmon$ wget --no-check-certificate https://raw.githubusercontent.com/COMBINE-lab/SalmonTools/master/scripts/generateDecoyTranscriptome.sh
# 安装mashmap和bedtools
/mnt/d/rnaseq/salmon$ conda install -y -c bioconda mashmap bedtools
# 查看脚本文件用法
/mnt/d/rnaseq/salmon$ bash generateDecoyTranscriptome.sh
****************
*** getDecoy ***
****************
Error: missing required argument(s)
Usage: generateDecoyTranscriptome.sh [-j <N> =1 default] [-b <bedtools binary path> =bedtools default] [-m <mashmap binary path> =mashmap default] -a <gtf file> -g <genome fasta> -t <txome fasta> -o <output path>
***************
*** ABORTED ***
***************
# 参数解释如下
> 1、-j # 线程数,设置不要太大,mashmap比较耗运行内存,小鼠的需要10多个G!大了容易崩
> 2、-b # bedtools路径
> 3、-m # mashmap路径
> 4、-a # gtf文件
> 5、-g # 基因组文件
> 6、-t # 环路组序列文件
> 7、-o` #输出文件夹,自动创建
# 下载gtf注释文件
/mnt/d/rnaseq/salmon$ wget http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M27/gencode.vM27.annotation.gtf.gz
# 解压
/mnt/d/rnaseq/salmon$ gunzip gencode.vM27.annotation.gtf.gz
# 下载基因组文件
/mnt/d/rnaseq/salmon$ wget http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M27/GRCm39.primary_assembly.genome.fa.gz
# 解压
/mnt/d/rnaseq/salmon$ gunzip GRCm39.primary_assembly.genome.fa.gz
小插曲--conda 安装遇到报错及解决:
/mnt/d/rnaseq/salmon$ conda install -y -c bioconda mashmap bedtools
Collecting package metadata (current_repodata.json): failed
CondaHTTPError: HTTP 000 CONNECTION FAILED for url <http://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/linux-64/current_repodata.json>
Elapsed: -
An HTTP error occurred when trying to retrieve this URL.
HTTP errors are often intermittent, and a simple retry will get you on your way.
'http://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/linux-64'
# 修改condarc文件,在镜像地址后面添加win-64就可以了
/mnt/d/rnaseq/salmon$ cd && vi .condarc
channels:
- http://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/win-64
- http://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/win-64
- http://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/win-64
- http://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/win-64
show_channel_urls: true
ssl_verify: true
mashmap 报错:
mashmap 需要依赖这 4 个包,如果报错可以尝试把这几个软件安装一下,太难了哈哈。
/mnt/d/rnaseq/salmon$ conda config --add channels conda-forge
/mnt/d/rnaseq/salmon$ conda install gsl
/mnt/d/rnaseq/salmon$ conda install libcxx libcxxabi
/mnt/d/rnaseq/salmon$ conda install zlib
2、获取 decoy 信息
# 查看bedtools可执行文件路径
/mnt/d/rnaseq/salmon$ which bedtools
/root/miniconda3/bin/bedtools
# 查看mashmap可执行文件路径
/mnt/d/rnaseq/salmon$ which mashmap
/root/miniconda3/bin/mashmap
# 获取decoy文件
/mnt/d/rnaseq/salmon$ bash generateDecoyTranscriptome.sh \
-j 1 \
-b /root/miniconda3/bin/bedtools \
-m /root/miniconda3/bin/mashmap \
-a ./gencode.vM27.annotation.gtf \
-g ./GRCm39.primary_assembly.genome.fa \
-t ./gencode.vM27.transcripts.fa.gz \
-o decoy_file
# 后面直接崩了,试了几次还是崩!内存还是太小,去大服务器跑完后
/mnt/d/rnaseq/salmon$ tree decoy_file
decoy_file
├── decoys.txt
└── gentrome.fa
# 完整日志文件
****************
*** getDecoy ***
****************
-j <Concurrency level> = 10
-b <bedtools binary> = /home/zhoulab/anaconda3/envs/salmon/bin/bedtools
-m <mashmap binary> = /home/zhoulab/anaconda3/envs/salmon/bin/mashmap
-a <Annotation GTF file> = /home/zhoulab/salmon/gencode.vM27.annotation.gtf
-g <Genome fasta> = /home/zhoulab/salmon/GRCm39.primary_assembly.genome.fa
-t <Transcriptome fasta> = /home/zhoulab/salmon/gencode.vM27.transcripts.fa.gz
-o <Output files Path> = decoy_file
[1/10] Extracting exonic features from the gtf
[2/10] Masking the genome fasta
[3/10] Aligning transcriptome to genome
>>>>>>>>>>>>>>>>>>
Reference = [reference.masked.genome.fa]
Query = [/home/zhoulab/salmon/gencode.vM27.transcripts.fa.gz]
Kmer size = 16
Window size = 5
Segment length = 500 (read split allowed)
Alphabet = DNA
Percentage identity threshold = 80%
Mapping output file = mashmap.out
Filter mode = 1 (1 = map, 2 = one-to-one, 3 = none)
Execution threads = 10
>>>>>>>>>>>>>>>>>>
INFO, skch::Sketch::build, minimizers picked from reference = 844163332
INFO, skch::Sketch::index, unique minimizers = 276702648
INFO, skch::Sketch::computeFreqHist, Frequency histogram of minimizers = (1, 141693231) ... (2608258, 1)
INFO, skch::Sketch::computeFreqHist, With threshold 0.001%, ignore minimizers occurring >= 7364 times during lookup.
INFO, skch::main, Time spent computing the reference index: 346.634 sec
INFO, skch::Map::mapQuery, [count of mapped reads, reads qualified for mapping, total input reads] = [111792, 111958, 142375]
INFO, skch::main, Time spent mapping the query : 3487.89 sec
INFO, skch::main, mapping results saved in : mashmap.out
[4/10] Extracting intervals from mashmap alignments
[5/10] Merging the intervals
[6/10] Extracting sequences from the genome
index file reference.masked.genome.fa.fai not found, generating...
[7/10] Concatenating to get decoy sequences
[8/10] Making gentrome
[9/10] Extracting decoy sequence ids
[10/10] Removing temporary files
**********************************************
*** DONE Processing ...
*** You can use files `$outfolder/gentrome.fa`
*** and $outfolder/decoys.txt` with
*** `salmon index`
**********************************************
我用 一个线程跑都崩了
,我电脑配置是 12 线程 12G 运行内存,主要是 mashmap 需要很多的内存,也和物种基因组大小有关
,后来直接拿到大服务器跑(200G运行内存
),消耗高达 64G!跑了大概半个多小时吧。结束后在 decoy_file 会输出 decoy.txt 文件
3、建立索引
/mnt/d/rnaseq/salmon$ salmon index -p 10 -k 31 --gencode \
-t gencode.vM27.transcripts.fa.gz \
-i salmon_index_partial_decoy \
--decoys decoy_file/decoys.txt
# 报错了
[2021-06-15 20:24:24.893] [puff::index::jointLog] [critical] The decoy file contained the names of 49 decoy sequences, but 0 were matched by sequences in the reference file provided. To prevent unintentional errors downstream, please ensure that the decoy file exactly matches with the fasta file that is being indexed.
[2021-06-15 20:24:24.944] [puff::index::jointLog] [error] The fixFasta phase failed with exit code 1
不知道啥原因,现在用的是最新版 1.5.0 的,是不是版本问题?于是下载了 1.0.0 的二进制版试一试:
/mnt/d/rnaseq/salmon$ wget https://github.com/COMBINE-lab/salmon/releases/download/v1.0.0/salmon-1.0.0_linux_x86_64.tar.gz
# 解压
/mnt/d/rnaseq/salmon$ tar -zxvf salmon-1.0.0_linux_x86_64.tar.gz
# 可执行文件在salmon-latest_linux_x86_64/bin/目录下
/mnt/d/rnaseq/salmon$ salmon-latest_linux_x86_64/bin/salmon
salmon v1.0.0
Usage: salmon -h|--help or
salmon -v|--version or
salmon -c|--cite or
salmon [--no-version-check] <COMMAND> [-h | options]
# 再尝试一下竟然可以了,是不是1.0版本以后都推荐使用第三种模式了?时间也要半个多小时
/mnt/d/rnaseq/salmon$ salmon-latest_linux_x86_64/bin/salmon \
index -p 10 -k 31 --gencode \
-t gencode.vM27.transcripts.fa.gz \
-i salmon_index_partial_decoy \
--decoys decoy_file/decoys.txt
...
[2021-06-15 20:57:08.770] [puff::index::jointLog] [info] finished populating pos vector
[2021-06-15 20:57:08.770] [puff::index::jointLog] [info] writing index components
[2021-06-15 20:57:10.139] [puff::index::jointLog] [info] finished writing dense pufferfish index
[2021-06-15 20:57:10.191] [jLog] [info] done building index
# 查看输出文件
/mnt/d/rnaseq/salmon$ ls -lh salmon_index_partial_decoy/
total 647M
-rwxrwxrwx 1 root root 550K Jun 15 20:59 complete_ref_lens.bin
-rwxrwxrwx 1 root root 36M Jun 15 21:45 ctable.bin
-rwxrwxrwx 1 root root 2.3M Jun 15 21:45 ctg_offsets.bin
-rwxrwxrwx 1 root root 66K Jun 15 20:59 duplicate_clusters.tsv
-rwxrwxrwx 1 root root 19M Jun 15 21:45 eqtable.bin
-rwxrwxrwx 1 root root 988 Jun 15 21:45 info.json
-rwxrwxrwx 1 root root 76M Jun 15 21:45 mphf.bin
-rwxrwxrwx 1 root root 403M Jun 15 21:45 pos.bin
-rwxrwxrwx 1 root root 115 Jun 15 21:45 pre_indexing.log
-rwxrwxrwx 1 root root 18M Jun 15 21:45 rank.bin
-rwxrwxrwx 1 root root 1.1M Jun 15 21:45 refAccumLengths.bin
-rwxrwxrwx 1 root root 17K Jun 15 21:45 ref_indexing.log
-rwxrwxrwx 1 root root 550K Jun 15 21:45 reflengths.bin
-rwxrwxrwx 1 root root 58M Jun 15 21:45 refseq.bin
-rwxrwxrwx 1 root root 35M Jun 15 21:45 seq.bin
-rwxrwxrwx 1 root root 96 Jun 15 21:45 versionInfo.json
我发现-t 参数跟 decoy_file 文件夹下的 gentrome.fa 也可以建索引,参数为-t decoy_file/gentrome.fa
,这样就有点像下面要介绍的模式了。
4、定量
定量基本和上面是一样的,只需要把 salmon quant 的 -i 参数 salmon_index 替换成 salmon_index_partial_decoy 就可以了。
|| A、mapping-based mode——SAF genome index ||
这个模式将构建最大的索引,也是最准确的方法,作者推荐哦。
1、准备文件
# 获取基因组序列名称
/mnt/d/rnaseq/salmon$ grep "^>" <(gunzip -c GRCm39.primary_assembly.genome.fa.gz) | cut -d " " -f 1 > decoys.txt
# 把>符号去掉
/mnt/d/rnaseq/salmon$ sed -i.bak -e 's/>//g' decoys.txt
# 合并基因组和转录组
/mnt/d/rnaseq/salmon$ cat gencode.vM27.transcripts.fa.gz GRCm39.primary_assembly.genome.fa.gz > gentrome.fa.gz
2、构建索引
可试读81%