Salmon是不基于比对计数而直接对基因进行定量的工具,适用于转录组、宏基因组等的分析。
优势在于:
- 定量时考虑到不同样品中基因长度的改变(比如不同isoform的使用)
- 速度快、需要的计算资源和存储资源小
- 敏感性高,不会丢弃匹配到多个基因同源区域的reads
- 可以直接校正GC-bias
- 自动判断文库类型
```
第一步:构建索引
ENSEMBL下载基因组和基因注释文件
mkdir -p genome cd genome
GRCh38.fa 人基因组序列,从Ensembl下载
GRCh38.gtf 人基因注释序列,从Ensembl下载
wget ftp://ftp.ensembl.org/pub/release-100/gtf/homo_sapiens/Homo_sapiens.GRCh38.100.gtf.gz -O GRCh38.fa.gz wget ftp://ftp.ensembl.org/pub/release-100/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz -O GRCh38.gtf.gz gunzip -c GRCh38.fa.gz >GRCh38.fa gunzip -c GRCh38.gtf.gz >GRCh38.gtf
获取cDNA序列
gffread GRCh38.gtf -g GRCh38.fa -w GRCh38.transcript.fa.tmp
gffread生成的fasta文件同时包含基因名字和转录本名字
grep ‘>’ GRCh38.transcript.fa.tmp | head
去掉空格后面的字符串,保证cDNA文件中fasta序列的名字简洁,不然后续会出错
cut -f 1 -d ‘ ‘ GRCh38.transcript.fa.tmp >GRCh38.transcript.fa
获取所有基因组序列的名字存储于decoy中
grep ‘^>’ GRCh38.fa | cut -d ‘ ‘ -f 1 | sed ‘s/^>//g’ >GRCh38.decoys.txt
合并cDNA和基因组序列一起
注意:cDNA在前,基因组在后
cat GRCh38.transcript.fa GRCh38.fa >GRCh38_trans_genome.fa
构建索引 (更慢,结果会更准)
salmon index -t GRCh38_trans_genome.fa -d GRCh38.decoys.txt -i GRCh38.salmon_sa_index
定量单样品FASTQ数据
cd ../ fastq-dump -v –split-3 –gzip SRR1039521 rename “SRR1039521” “trt_N061011” SRR1039521*
-p: 表示若待创建的文件夹已存在则跳过;若不存在,则创建;也可用于创建多层文件夹
man mkdir 可查看详细帮助
mkdir -p trt_N061011
-l: 自动判断文库类型,尤其适用于链特异性文库
The library type -l should be specified on the command line
before the read files (i.e. the parameters to -1 and -2, or -r).
This is because the contents of the library type flag is used to determine how the reads should be interpreted.
–gcBias: 校正测序片段GC含量,获得更准确的转录本定量结果
One can simply run Salmon with –gcBias in any case,
as it does not impair quantification for samples without GC bias,
it just takes a few more minutes per sample.
For samples with moderate to high GC bias, correction for this bias at the
fragment level has been shown to reduce isoform quantification errors
salmon quant –gcBias -l A -1 trt_N061011_1.fq.gz -2 trt_N061011_2.fq.gz -i genome/GRCh38.salmon_sa_index -g genome/GRCh38.gtf -o trt_N061011/trt_N061011.salmon.count -p 10
输出结果存储在 trt_N061011/trt_N061011.salmon.count目录中
quant.sf 为转录本表达定量结果,第4列为TPM结果,第5列为reads count
quant.genes.sf 为基因表达定量结果
head -n 30 trt_N061011/trt_N061011.salmon.count/quant.sf | tail
循环定量多个样品的表达量
for samp in tail -n +2 sampleFile | cut -f 1
; do salmon quant –gcBias -l A -1 ${samp}_1.fq.gz -2 ${samp}_2.fq.gz -i genome/GRCh38.salmon_sa_index -o ${samp}/${samp}.salmon.count -p 4 >${samp}.salmon.log 2>&1; done &
整理Salmon定量文件用于DESeq2差异基因鉴定
列出salmon的输出文件
find . -name quant.sf
这个压缩包下载解压到本地
zip quant.sf.zip find . -name quant.sf
生成一个两列文件方便R导入
xargs接收上一步的输出,按批次提供给下游程序作为输入
-i: 用{}表示传递的值
cut -f 1 sampleFile | xargs -i echo -e “{}\t{}/{}.salmon.count/quant.sf” >salmon.output head salmon.output
获得基因和转录本的对应关系,获取基因的表达量
如果没有GTF文件,可以用其他文件,只需获取转录本和基因名字对应关系就可以
如果不知道对应关系,也可以把每个转录本当做一个基因进行分析
Trinity拼装时会生成这个文件
注意修改$14, $10为对应的信息列,
tx2gene为一个两列文件,第一列是转录本没名字,第二列是基因名字。
sed ‘s/”/\t/g’ genome/GRCh38.gtf | awk ‘BEGIN{OFS=FS=”\t”}{if(FNR==1) print “TXname\tGene”; if($3==”transcript”) print $14, $10}’ >GRCh38.tx2gene head GRCh38.tx2gene ```