基于Salmon的轉(zhuǎn)錄組定量流程
為什么使用Salmon？
Salmon是不基于比對計數(shù)而直接對基因進(jìn)行定量的工具，適用于轉(zhuǎn)錄組、宏基因組等的分析。
其優(yōu)勢是：
定量時考慮到不同樣品中基因長度的改變（比如不同isoform的使用）
速度快、需要的計算資源和存儲資源小
敏感性高，不會丟棄匹配到多個基因同源區(qū)域的reads
可以直接校正GC-bias
自動判斷文庫類型
39個轉(zhuǎn)錄組分析工具，120種組合評估表明Salmon的定量準(zhǔn)確性和穩(wěn)定性都比較好。
其原理如下圖所示，概括來講是通過構(gòu)建統(tǒng)計模型來推測已經(jīng)注釋的轉(zhuǎn)錄本呈現(xiàn)出什么表達(dá)模式時我們才會測序產(chǎn)生當(dāng)前的FASTQ數(shù)據(jù)：
怎么使用Salmon？
Salmon定量依賴于cDNA序列和原始的FASTQ序列，新版本也可以提供基因組序列以處理某些能同時比對到已經(jīng)注釋的基因區(qū)和基因間區(qū)的reads，獲得更準(zhǔn)確地定量結(jié)果。
第一步，構(gòu)建索引
從ENSEMBL下載基因組和基因注釋文件，具體參考NGS基礎(chǔ) - 參考基因組和基因注釋文件。
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.fagunzip -c GRCh38.gtf.gz >GRCh38.gtf# 獲取cDNA序列
gffread GRCh38.gtf -g GRCh38.fa -w GRCh38.transcript.fa.tmp
# gffread生成的fasta文件同時包含基因名字和轉(zhuǎn)錄本名字
grep '>' GRCh38.transcript.fa.tmp | head
# 去掉空格后面的字符串，保證cDNA文件中fasta序列的名字簡潔，不然后續(xù)會出錯
cut -f 1 -d ' ' GRCh38.transcript.fa.tmp >GRCh38.transcript.fa
構(gòu)建索引
# 獲取所有基因組序列的名字存儲于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
# 構(gòu)建索引 （更慢，結(jié)果會更準(zhǔn)）
salmon index -t GRCh38_trans_genome.fa -d GRCh38.decoys.txt -i GRCh38.salmon_sa_index定量單樣品FASTQ數(shù)據(jù)
cd ../
fastq-dump -v --split-3 --gzip SRR1039521
rename "SRR1039521" "trt_N061011" SRR1039521*
# -p: 表示若待創(chuàng)建的文件夾已存在則跳過；若不存在，則創(chuàng)建；也可用于創(chuàng)建多層文件夾
# man mkdir 可查看詳細(xì)幫助
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含量，獲得更準(zhǔn)確的轉(zhuǎn)錄本定量結(jié)果
# 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定量后輸出結(jié)果存儲于trt_N061011/trt_N061011.salmon.count目錄中
# 輸出結(jié)果存儲在 trt_N061011/trt_N061011.salmon.count目錄中
# quant.sf 為轉(zhuǎn)錄本表達(dá)定量結(jié)果，第4列為TPM結(jié)果，第5列為reads count
# quant.genes.sf 為基因表達(dá)定量結(jié)果
head -n 30 trt_N061011/trt_N061011.salmon.count/quant.sf | tail定量結(jié)果為
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多樣品定量
未完待續(xù)
往期精品(點擊圖片直達(dá)文字對應(yīng)教程)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
后臺回復(fù)“生信寶典福利第一波”或點擊閱讀原文獲取教程合集
（請備注姓名-學(xué)校/企業(yè)-職務(wù)等）
評論
圖片
表情