基于Salmon的转录组定量流程

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

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

Name    Length    EffectiveLength    TPM    NumReads
ENST00000609179    1196    1052.656    0.000000    0.000
ENST00000492242    1277    1058.126    92.111195    19.918
ENST00000382291    2088    1963.212    1447.695765    580.820
ENST00000382285    1329    1183.099    211.657526    51.174