Salmon應(yīng)用

查看幫助文檔

#查看可用的命令
###Salmon v0.9.1
salmon -h
#查看幫助文檔之Salmon's quasi-mapping-based mode
salmon --no-version-check quant --help-reads
#查看幫助文檔之Salmon's alignment-based mode
salmon --no-version-check quant --help-alignment

Quasi-mapping-based mode (including lightweight alignment)

目前，Salmon有兩種不同的方法---支持將reads比對(duì)到transcriptomes，分別是(SMEM-based) lightweight-alignment 和 quasi-mapping。SMEM-based mapping是Salmon自帶的輕量級(jí)的比對(duì)方法；而 quasi-mapping這種比對(duì)方法相對(duì)較新較快。這兩種方法的使用均是通過quant命令來應(yīng)用的，但是兩者比對(duì)方法所應(yīng)用到的索引不同。

quasi-mapping 應(yīng)用到的數(shù)據(jù)結(jié)構(gòu)是：a combination of data structures—a hash table, suffix array (SA) and efficient rank data structure

quasi-mappings的詳情見：RapMap: a rapid, sensitive and accurate tool for mapping RNA-seq reads to transcriptomes

#構(gòu)建quasi-mapping-based index
##注意，此處的-k值的設(shè)定。
####當(dāng)reads長度長于75bp時(shí)，k值設(shè)為31是個(gè)理想的選擇；
####當(dāng)reads長度短于75bp時(shí)，應(yīng)該嘗試使用降低k值。
####另外，參數(shù)k能夠被傳遞給命令quant。因此，在構(gòu)建索引階段提供的參數(shù)k能夠覆蓋在估計(jì)表達(dá)量階段的參數(shù)k的設(shè)定。
######即如果用戶使用quasi-mapping-based index，則估計(jì)表達(dá)量階段的參數(shù)k的設(shè)定對(duì)結(jié)果沒有任何影響。
####因?yàn)閕ndex命令中，--type參數(shù)的默認(rèn)值是quasi，所以在構(gòu)建quasi-mapping-based index時(shí)，--type quasi可以省略不寫。
salmon index -t transcripts.fa -i transcripts_index --type quasi -k 31

#構(gòu)建SMEM-based mapping index 
#### 雖然此處未給定參數(shù)k,但是程序的默認(rèn)值是19。
#### 同上，在構(gòu)建索引階段提供的參數(shù)k能夠覆蓋在估計(jì)表達(dá)量階段的參數(shù)k的設(shè)定。
salmon index -t transcripts.fa -i transcripts_index --type fmd

#對(duì)雙端測(cè)序數(shù)據(jù)reads表達(dá)量的估計(jì)
salmon quant -i transcripts_index -l <LIBTYPE> -1 reads1.fq -2 reads2.fq -o transcripts_quant

#對(duì)單端測(cè)序數(shù)據(jù)reads表達(dá)量的估計(jì)
salmon quant -i transcripts_index -l <LIBTYPE> -r reads.fq -o transcripts_quant

### 命令quant均適用于這兩個(gè)index(quasi-mapping or SMEM-based)，此外，Salmon能夠自動(dòng)檢測(cè)到使用的是哪種index，從而采用與之匹配的比對(duì)方法。
### 注意：參數(shù)-l必須指定在參數(shù)-1，-2和-r的前面。
### 生成一個(gè)目錄，內(nèi)含文件quant.sf

Alignment-based mode

#使用aligner比對(duì)生成的結(jié)果文件aln.bam，對(duì)其進(jìn)行表達(dá)量分析
### 生成一個(gè)目錄，內(nèi)含文件quant.sf
salmon quant -t transcripts.fa -l <LIBTYPE> -a aln.bam -o salmon_quant

salmon分析RNAseq數(shù)據(jù)實(shí)戰(zhàn)

#創(chuàng)建分析文件夾
mkdir salmon_tutorial
cd salmon_tutorial

#下載參考轉(zhuǎn)錄組數(shù)據(jù)
curl ftp://ftp.ensemblgenomes.org/pub/plants/release-28/fasta/arabidopsis_thaliana/cdna/Arabidopsis_thaliana.TAIR10.28.cdna.all.fa.gz -o athal.fa.gz

#創(chuàng)建索引
salmon index -t athal.fa.gz -i athal_index

#獲取測(cè)序數(shù)據(jù)
mkdir data
cd data
for i in `seq 25 40`; 
do 
  mkdir DRR0161${i}; 
  cd DRR0161${i}; 
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/DRR016/DRR0161${i}/DRR0161${i}_1.fastq.gz; 
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/DRR016/DRR0161${i}/DRR0161${i}_2.fastq.gz; 
  cd ..; 
done
cd .. 

#估計(jì)樣本的表達(dá)量
for fn in data/DRR0161{25..40};
do
samp=`basename ${fn}`
echo "Processing sample ${samp}"
salmon quant -i athal_index -l A \
         -1 ${fn}/${samp}_1.fastq.gz \
         -2 ${fn}/${samp}_2.fastq.gz \
         -p 8 -o quants/${samp}_quant
done 

### 其他常用的參數(shù)
--seqBias   #糾偏序列的bias
--gcBias    #糾偏GC含量的bias[beta for single-end reads]
-g/--geneMap    #提供轉(zhuǎn)錄本與基因之間對(duì)應(yīng)關(guān)系的文件。如果提供了該文件，Salmon將會(huì)輸出文件：quant.sf 和 quant.genes.sf。
                #其中quant.genes.sf 文件中含有對(duì)基因表達(dá)水平的估計(jì)。
                #提供的文件可以是GTF文件，也可以是由TAB分隔符分隔文件，文件中的每一行含轉(zhuǎn)錄本名稱和對(duì)應(yīng)的基因名。
                #以'.gtf', '.gff' 或 '.gff3'結(jié)尾的文件名均被當(dāng)做GTF格式做處理解析。
                #以其他字符結(jié)尾的被當(dāng)做由TAB分隔符分隔的簡單文件格式做處理解析。
-z/--writeMappings      #若提供了該參數(shù)，則quasi-mapping的結(jié)果也會(huì)以SAM格式輸出。
                        #默認(rèn)情況下，輸出結(jié)果是標(biāo)準(zhǔn)輸出，建議提供文件名存儲(chǔ)至文件中，以便后續(xù)分析的需要。
### 其他常用高級(jí)的參數(shù)                       
--useVBOpt #使用Variational Bayesian EM算法，而非傳統(tǒng)的EM算法。目的是批處理優(yōu)化
--numBootstraps num     #Number of bootstrap samples to generate。
                        #注意，該參數(shù)的使用與Gibbs sampling是相斥的。
                        #如果后續(xù)的分析差異表達(dá)搭配sleuth使用，則必須應(yīng)用此參數(shù)。


#下游的差異表達(dá)分析
### 差異表達(dá)分析的包：DESeq2, edgeR, limma, or sleuth
### 使用tximport包導(dǎo)入salmon轉(zhuǎn)錄本表達(dá)量分析的結(jié)果文件
### 再根據(jù) DESeq2 vignette（https://bioconductor.org/packages/DESeq2）提供的范例，將salmon轉(zhuǎn)錄本表達(dá)量分析的結(jié)果文件讀入DESeq2

Salmon輸出文件

Salmon輸出文件之Quantification File

Quantification File： quant.sf
quant.sf文件有5列，分別是Name，Length ，EffectiveLength，TPM和NumReads。分別表示的含義如下所述：

• Name — target transcript 名稱， 由輸入的 transcript database (FASTA file)所提供。
• Length — target transcript 長度，即有多少個(gè)核苷酸
• EffectiveLength — target transcript 計(jì)算的有效長度。此項(xiàng)考慮了所有被建模的因素，這將影響從這個(gè)轉(zhuǎn)錄本中取樣片段的概率，包括片段長度分布和序列特異性和gc片段偏差（如果這些因素在建模時(shí)均被考慮的話）。 （It takes into account all factors being modeled that will effect the probability of sampling fragments from this transcript, including the fragment length distribution and sequence-specific and gc-fragment bias (if they are being modeled)）。
• TPM — 估計(jì)轉(zhuǎn)錄本的表達(dá)量。
• NumReads — 估計(jì)比對(duì)到每個(gè)轉(zhuǎn)錄本的reads數(shù)。

補(bǔ)充：FPKM，RPKM，RPM以及TPM的關(guān)系之見解
RPKM (Reads Per Kilobase Million)
FPKM (Fragments Per Kilobase Million)
TPM(Transcripts Per Kilobase Million)
RPM (Reads per million)
CPM (Counts per million)

TPM：
Transcripts Per Kilobase of exonmodel per Million mapped reads (每千個(gè)堿基的轉(zhuǎn)錄每百萬映射讀取的Transcripts)，優(yōu)化的RPKM計(jì)算方法，可以用于同一物種不同組織的比較。
TPM概括了基因的長度、表達(dá)量和基因數(shù)目。TPM可以用于同一物種不同組織間的比較，因?yàn)閟um值總是唯一的。

Salmon輸出文件之其他

cmd_info.json： JSON格式文件，記錄salmon程序運(yùn)行的命令和參數(shù)
lib_format_counts.json： Observed library format counts。當(dāng)運(yùn)行salmon是 mapping-based mode時(shí)，則會(huì)生成改文件。 JSON格式文件，記錄有關(guān)文庫格式和reads比對(duì)的情況。
eq_classes.txt： Equivalence class file。當(dāng)Salmon運(yùn)行時(shí)，應(yīng)用參數(shù)--dumpEq，則會(huì)生成此文件。
aux_info： 輔助文件夾，內(nèi)含多個(gè)文件
fld.gz：在輔助文件夾中，該文件記錄的是觀察到的片段長度分布的近似值
obs5_seq.gz, obs3_seq.gz, exp5_seq.gz, exp5_seq.gz： Sequence-specific bias files
expected_gc.gz, observed_gc.gz： 當(dāng)Salmon運(yùn)行時(shí)，應(yīng)用fragment-GC bias correction，在輔助文件夾中則會(huì)生成這兩個(gè)文件。記錄Fragment-GC bias。
meta_info.json： JSON格式文件，記錄salmon程序運(yùn)行的統(tǒng)計(jì)信息
ambig_info.tsv： tab分隔符的文本文件，含有兩列。記錄的是每個(gè)轉(zhuǎn)錄本對(duì)應(yīng)的 the number of uniquely-mapping reads 和 the total number of ambiguously-mapping reads

參考：Salmon---a tool for wicked-fast transcript quantification from RNA-seq data