使用RSEM进行转录组测序的差异表达分析

仍然是两年前的笔记

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的记录提取出来

因水平有限，有错误的地方，欢迎批评指正！