Overview:
RNA测序结果首先在Linux/Unix中通过命令行得到每个基因的读数(read counts),接着在R中利用R包进行统计分析,找到差异表达的基因,从而发掘生物意义。这部分内容所要介绍的就是如何通过R来找到差异表达的基因。
开始前的设置
在开始数据处理与分析前,有条理地建立文件夹与目录非常重要,下面将简述开始前的设置流程:
1.建立R的工程(Project),建立新的Directory。
2.依次在Directory中建立三个子集:data,meta和results,分别用来储存原数据、中间处理的数据以及最后的结果。
3.载入library,对于没有安装的R包,用install.package("")指令安装。
4.载入数据
## Load in data
data<-read.table("data/Mov10_full_counts.txt",header=T,row.names=1)
meta<-read.table("meta/Mov10_full_meta.txt",header=T,row.names=1)
5.读取数据:知道数据的类型,并简单查看数据。
RNA-seq读数分布
The goal of differential expression analysis is to determine, for each gene, whether the differences in expression (counts) between groups is significant given the amount of variation observed within groups (replicates).
将基因表达读数(counts)绘制成柱状图,就会发现,很大一部分基因的读数接近于0,这就是RNA测序数据的特点:与大量基因的读数很低,由于缺乏表达上限而导致的长右尾。
ggplot(data) +
geom_histogram(aes(x = Mov10_oe_1), stat = "bin", bins = 200) +
xlab("Raw expression counts") +
ylab("Number of genes")
若将x轴放大一些,可以更加明显地看出低读数基因数目之多。
ggplot(data) +
geom_histogram(aes(x = Mov10_oe_1), stat = "bin", bins = 200) +
xlim(-5, 500) +
xlab("Raw expression counts") +
ylab("Number of genes")
RNA测序数据的分布拟合
单个样本的RNA测序数据通常认为是符合泊松分布的,满足分布的一个前提是均值和方差相等,但是生物学样本总存在差异,在重复组之间存在变异性,此时负二项式模型较为合适,特点就是均值小于方差。
泊松分布容易低估变异性,导致假阳性差异表达基因的出现,这也是为什么需要有生物学重复的原因,用生物学重复能够提高平均值、减少方差,避免假阳性的出现。
本学习笔记主要参考了Github上的DGE_workshop教程