Salmon

基于Salmon的转录组定量流程

Posted by CHY on December 12, 2020

Salmon是不基于比对计数而直接对基因进行定量的工具,适用于转录组、宏基因组等的分析。
优势在于:

  1. 定量时考虑到不同样品中基因长度的改变(比如不同isoform的使用)
  2. 速度快、需要的计算资源和存储资源小
  3. 敏感性高,不会丢弃匹配到多个基因同源区域的reads
  4. 可以直接校正GC-bias
  5. 自动判断文库类型 salmon ```

    第一步:构建索引

    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 ```