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