使用RSEM進行轉錄組測序的差異表達分析

TOP生物資訊發表於2021-02-21

仍然是兩年前的筆記

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的記錄提取出來

因水平有限，有錯誤的地方，歡迎批評指正！

相關文章