仍然是兩年前的筆記
1. prepare-reference
如果用RSEM對比對後的bam進行轉錄本定量，則在比對過程中要確保比對用到的索引是由rsem-prepare-reference產生的。
~/software/rsem/rsem-prepare-reference \
--transcript-to-gene-map ~/project/RNA-seq/ref_cds/gene_transcript.txt \ #作用是在後面的定量結果檔案中，新增gene名稱, 轉錄本名稱兩列，該檔案每一行都是gene_id\ttranscript_id的形式，eg: cluster_11236 cluster_11236.1
--bowtie2 \ #RSEM可呼叫bowtie, bowtie2, STAR三種比對工具；這裡選用bowtie2
~/project/RNA-seq/ref_cds/HC_cds_and_8sample_clustercds.fa \
~/project/RNA-seq/ref_cds/cds.byrsem
可以看到，單純用bowtie2建的索引和rsem呼叫bowtie2建的索引是不一樣的。
2. calculate-expression
用法分為兩類，分別是從fa/fq得到表達矩陣，和從sam/bam得到表達矩陣（仍然要求是比對到rsem-prepare-reference生成的索引）。以單端的fq資料為例。
rsem-calculate-expression [options] upstream_read_file(s) reference_name sample_name
rsem-calculate-expression [options] --paired-end upstream_read_file(s) downstream_read_file(s) reference_name sample_name
rsem-calculate-expression [options] --sam/--bam [--paired-end] input reference_name sample_name
cat ~/project/RNA-seq/dir.txt | while read id
do
~/software/rsem/rsem-calculate-expression -p 8 --bowtie2 \
~/project/data/RNA-seq/${id}.fastq.gz \
~/project/RNA-seq/ref_cds/cds.byrsem \
--samtools-sort-mem 2G --fragment-length-mean 50 \ # 單端資料建議使用--fragment-length-mean和--fragment-length-sd
~/project/RNA-seq/map/${id}.rsem
done
完成之後得到這些檔案，其中，rsem.genes.results和rsem.isoforms.results分別表示gene水平和轉錄本水平的定量結果。每一列含義：
less rsem.genes.results|head -n 1
gene_id transcript_id(s) length effective_length expected_count TPM FPKM
less rsem.isoforms.results|head -n 1
transcript_id gene_id length effective_length expected_count TPM FPKM IsoPct
後面用EBseq檢驗差異基因/轉錄本時，會使用到這兩個檔案。
3. Differential Expression Analysis using EBSeq
下面以gene水平差異分析為例。
3.1 generate-data-matrix
這一步提取上一步得到的每個樣本定量結果檔案中的expected_count列，組成資料矩陣。
~/software/rsem/rsem-generate-data-matrix \
SRR1.rsem.genes.results SRR2.rsem.genes.results \
SRR3.rsem.genes.results SRR4.rsem.genes.results \
SRR5.rsem.genes.results SRR6.rsem.genes.results \
SRR7.rsem.genes.results SRR8.rsem.genes.results \
> ~/project/RNA-seq/count/GeneMat.txt
3.2 run-ebseq
呼叫EBseq進行檢驗
~/software/rsem/rsem-run-ebseq \
GeneMat.txt 2,2,2,2 GeneMat.results #2,2,2,2表示4個condition, 每個condition有兩個重複；順序要和3.1中輸入檔案表示的condition的順序一致
#會得到三個檔案
GeneMat.results.condmeans GeneMat.results GeneMat.results.pattern
#GeneMat.results.pattern
"C1" "C2" "C3" "C4"
"Pattern1" 1 1 1 1
"Pattern2" 1 1 1 2
"Pattern3" 1 1 2 1
"Pattern4" 1 1 2 2
"Pattern5" 1 2 1 1
"Pattern6" 1 2 1 2
"Pattern7" 1 2 2 1
"Pattern8" 1 2 2 2
"Pattern9" 1 1 2 3
"Pattern10" 1 2 1 3
"Pattern11" 1 2 2 3
"Pattern12" 1 2 3 1
"Pattern13" 1 2 3 2
"Pattern14" 1 2 3 3
"Pattern15" 1 2 3 4
#以Pattern14為例，1 2 3 3表示某基因表達：C1與C2不同，C3與C4相同
#四種condition如果有基因表達存在差異，就這些情況了
#GeneMat.results
#第一列是各個基因名稱，接著15列是該基因符合該種Parttern的概率
#"MAP"為該基因最可能的模式；"PPDE"：posterior probability of being differentially expressed，越大越好
"Pattern1" "Pattern2" "Pattern3" "Pattern4" "Pattern5" "Pattern6" "Pattern7" "Pattern8" "Pattern9" "Pattern10" "Pattern11" "Pattern12" "Pattern13" "Pattern14" "Pattern15" "MAP" "PPDE"
#GeneMat.results.condmeans
#為每個樣本合併重複之後的定量結果，如下圖，這個結果可以用來控制fold change
3.3 control_fdr
控制FDR(錯誤發現率)來挑選差異基因
~/software/rsem/rsem-control-fdr \
GeneMat.results 0.05 GeneMat.de.txt
將GeneMat.results檔案中，PPDE大於0.95的記錄提取出來
因水平有限，有錯誤的地方，歡迎批評指正！