• 如果您觉得本站非常有看点,那么赶紧使用Ctrl+D 收藏吧

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

开发技术 开发技术 3天前 4次浏览

仍然是两年前的笔记

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_idttranscript_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

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

可以看到,单纯用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进行转录组测序的差异表达分析

完成之后得到这些文件,其中,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

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

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

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

3.3 control_fdr

控制FDR(错误发现率)来挑选差异基因

~/software/rsem/rsem-control-fdr 
GeneMat.results 0.05 GeneMat.de.txt

将GeneMat.results文件中,PPDE大于0.95的记录提取出来

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


喜欢 (0)