library(limma)
library(Glimma)
library(edgeR)
library(Mus.musculus)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
# 数据下载
url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"
utils::download.file(url, destfile="GSE63310_RAW.tar", mode="wb")
utils::untar("GSE63310_RAW.tar", exdir = ".")
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
"GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
"GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")
for(i in paste(files, ".gz", sep=""))
R.utils::gunzip(i, overwrite=TRUE)
# 读取文件
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt",
"GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt",
"GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
"GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt",
"GSM1545545_JMS9-P8c.txt")
read.delim(files[1], nrow=5)
#readDGE函数即可一次性读取所有文件,生成DGEList对象中包含一个计数矩阵，
#它的27179行分别对应唯一的Entrez基因标识（ID），九列分别对应此实验中的每个样品。
x <- readDGE(files, columns=c(1,3))
class(x)
# meta data
samplenames <- substring(colnames(x), 12, nchar(colnames(x)))
samplenames
colnames(x) <- samplenames
group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP",
"Basal", "ML", "LP"))
x$samples$group <- group
lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2)))
x$samples$lane <- lane
x$samples
> x$samples
files group lib.size norm.factors lane
10_6_5_11 GSM1545535_10_6_5_11.txt LP 32863052 1 L004
9_6_5_11 GSM1545536_9_6_5_11.txt ML 35335491 1 L004
purep53 GSM1545538_purep53.txt Basal 57160817 1 L004
JMS8-2 GSM1545539_JMS8-2.txt Basal 51368625 1 L006
JMS8-3 GSM1545540_JMS8-3.txt ML 75795034 1 L006
JMS8-4 GSM1545541_JMS8-4.txt LP 60517657 1 L006
JMS8-5 GSM1545542_JMS8-5.txt Basal 55086324 1 L006
JMS9-P7c GSM1545544_JMS9-P7c.txt ML 21311068 1 L008
JMS9-P8c GSM1545545_JMS9-P8c.txt LP 19958838 1 L008
# 基因注释
geneid <- rownames(x)
genes <- select(Mus.musculus, keys=geneid, columns=c("SYMBOL", "TXCHROM"),
keytype="ENTREZID")
head(genes)
> head(genes)
ENTREZID SYMBOL TXCHROM
1 497097 Xkr4 chr1
2 100503874 Gm19938 <NA>
3 100038431 Gm10568 <NA>
4 19888 Rp1 chr1
5 20671 Sox17 chr1
6 27395 Mrpl15 chr1
# 去重复，
# 保留每个基因ID第一次出现的信息
genes <- genes[!duplicated(genes$ENTREZID),]
#使用edgeR中的cpm函数将原始计数转换为CPM和log-CPM值。
#如果可以提供基因长度信息，可以使用edgeR中的rpkm函数计算RPKM值
#对于一个基因，CPM值为1相当于在测序深度最低的样品中（JMS9-P8c,
#序列数量约2千万）有20个计数，
#或者在测序深度最高的样品中（JMS8-3，序列数量约7.6千万）有76个计数。
cpm <- cpm(x)
lcpm <- cpm(x, log=TRUE, prior.count=2)
L <- mean(x$samples$lib.size) * 1e-6
M <- median(x$samples$lib.size) * 1e-6
c(L, M)
[1] 45.48855 51.36862
# 删除低表达的基因
table(rowSums(x$counts==0)==9)
FALSE TRUE
22026 5153
# edgeR包中的filterByExpr函数提供了自动过滤基因的方法，可保留尽可能多的有足够表达计数的基因。
keep.exprs <- filterByExpr(x, group=group)
x <- x[keep.exprs,, keep.lib.sizes=FALSE]
dim(x)
[1] 16624 9
#此函数默认选取最小的组内的样本数量为最小样本数，保留至少在这个数量的样本中有10个或更多序列片段计数的基因。
#对基因表达量进行过滤时使用CPM值而不是表达计数，以避免对总序列数大的样本的偏向性。
#在这个数据集中，总序列数的中位数是5.1千万，且10/51约等于0.2，所以filterByExpr函数保留在至少三个样本中CPM值大于等于0.2的基因。
#此处，一个具有生物学意义的基因需要在至少三个样本中表达，因为三种细胞类型组内各有三个重复。
#所使用的阈值取决于测序深度和实验设计。如果样本总表达计数量增大，那么可以选择更低的CPM阈值，
#因为更大的总表达计数量提供了更好的分辨率来探究更多表达水平更低的基因。
#使用这个标准，基因的数量减少到了16624个，约为开始时数量的60%。
#过滤后的log-CPM值显示出每个样本的分布基本相同（下图B部分）。
#需要注意的是，从整个DGEList对象中取子集时同时删除了被过滤的基因的计数和其相关的基因信息。
#过滤后的DGEList对象为留下的基因保留了相对应的基因信息和计数。
lcpm.cutoff <- log2(10/M + 2/L)
library(RColorBrewer)
nsamples <- ncol(x)
col <- brewer.pal(nsamples, "Paired")
par(mfrow=c(1,2))
plot(density(lcpm[,1]), col=col[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="")
title(main="A. Raw data", xlab="Log-cpm")
abline(v=lcpm.cutoff, lty=3)
for (i in 2:nsamples){
den <- density(lcpm[,i])
lines(den$x, den$y, col=col[i], lwd=2)
}
legend("topright", samplenames, text.col=col, bty="n")
lcpm <- cpm(x, log=TRUE)
plot(density(lcpm[,1]), col=col[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="")
title(main="B. Filtered data", xlab="Log-cpm")
abline(v=lcpm.cutoff, lty=3)
for (i in 2:nsamples){
den <- density(lcpm[,i])
lines(den$x, den$y, col=col[i], lwd=2)
}
legend("topright", samplenames, text.col=col, bty="n")
image.png
- Figure 1: 每个样本过滤前的原始数据（A）和过滤后（B）的数据的log-CPM值密度。竖直虚线标出了过滤步骤中所用阈值（相当于CPM值为约0.2）。
#使用edgeR中的calcNormFactors函数，用TMM(Robinson and Oshlack 2010)方法进行归一化。
#此处计算得到的归一化系数被用作文库大小的缩放系数。
#当我们使用DGEList对象时，这些归一化系数被自动存储在x$samples$norm.factors。
#对此数据集而言，TMM归一化的作用比较温和，这体现在所有的缩放因子都相对接近1。
x <- calcNormFactors(x, method = "TMM")
x$samples$norm.factors
#为了更好地可视化表现出归一化的影响，我们复制了数据并进行了调整，
#使得第一个样品的计数减少到了其原始值的5%，而第二个样品增大到了5倍。
x2 <- x
x2$samples$norm.factors <- 1
x2$counts[,1] <- ceiling(x2$counts[,1]*0.05)
x2$counts[,2] <- x2$counts[,2]*5
par(mfrow=c(1,2))
lcpm <- cpm(x2, log=TRUE)
boxplot(lcpm, las=2, col=col, main="")
title(main="A. Example: Unnormalised data",ylab="Log-cpm")
x2 <- calcNormFactors(x2)
x2$samples$norm.factors
image.png
- Figure 2: 样例数据：log-CPM值的箱线图展示了未经归一化的数据（A）及归一化后的数据（B）中每个样本的表达分布。数据集经过调整，样本1和2中的表达计数分别被缩放到其原始值的5%和500%。
lcpm <- cpm(x, log=TRUE)
par(mfrow=c(1,2))
col.group <- group
levels(col.group) <- brewer.pal(nlevels(col.group), "Set1")
col.group <- as.character(col.group)
col.lane <- lane
levels(col.lane) <- brewer.pal(nlevels(col.lane), "Set2")
col.lane <- as.character(col.lane)
plotMDS(lcpm, labels=group, col=col.group)
title(main="A. Sample groups")
plotMDS(lcpm, labels=lane, col=col.lane, dim=c(3,4))
title(main="B. Sequencing lanes")
image.png
Figure 3: log-CPM值在维度1和2的MDS图，以样品分组上色并标记（A）和维度3和4的MDS图，以测序道上色并标记（B）。图中的距离对应于最主要的倍数变化（fold change），默认情况下也就是前500个在每对样品之间差异最大的基因的平均（均方根）log2倍数变化。
design <- model.matrix(~0+group+lane)
colnames(design) <- gsub("group", "", colnames(design))
design
contr.matrix <- makeContrasts(
BasalvsLP = Basal-LP,
BasalvsML = Basal - ML,
LPvsML = LP - ML,
levels = colnames(design))
contr.matrix
par(mfrow=c(1,2))
v <- voom(x, design, plot=TRUE)
v
vfit <- lmFit(v, design)
vfit <- contrasts.fit(vfit, contrasts=contr.matrix)
efit <- eBayes(vfit)
plotSA(efit, main="Final model: Mean-variance trend")
summary(decideTests(efit))
BasalvsLP BasalvsML LPvsML
Down 4648 4927 3135
NotSig 7113 7026 10972
Up 4863 4671 2517
image.png
- Figure 4: 图中绘制了每个基因的均值（x轴）和方差（y轴），显示了在该数据上使用voom前它们之间的相关性（左），以及当运用voom的精确权重后这种趋势是如何消除的（右）。左侧的图是使用voom函数绘制的，它为进行log-CPM转换后的数据拟合线性模型从而提取残差方差。然后，对方差取平方根（或对标准差取平方根），并相对每个基因的平均表达作图。均值通过平均计数加上2再进行log2转换计算得到。右侧的图使用plotSA绘制了log2残差标准差与log-CPM均值的关系。平均log2残差标准差由水平蓝线标出。在这两幅图中，每个黑点表示一个基因，红线为对这些点的拟合。
summary(decideTests(efit))
#一些研究中不仅仅需要使用校正p值阈值，更为严格定义的显著性可能需要差异倍数的对数（log-FCs）也高于某个最小值。
# treat方法(McCarthy and Smyth 2009)可以按照对最小log-FC值的要求，使用经过经验贝叶斯调整的t统计值计算p值。
#当我们的检验要求基因的log-FC显著大于1（等同于在原本的尺度上不同细胞类型之间差两倍）时，差异表达基因的数量得到了下降，
# basal与LP相比只有3684个DE基因，basal与ML相比只有3834个DE基因，LP与ML相比只有414个DE基因。
tfit <- treat(vfit, lfc=1)
dt <- decideTests(tfit)
summary(dt)
BasalvsLP BasalvsML LPvsML
Down 1632 1777 224
NotSig 12976 12790 16210
Up 2016 2057 190
de.common <- which(dt[,1]!=0 & dt[,2]!=0)
length(de.common)
head(tfit$genes$SYMBOL[de.common], n=20)
vennDiagram(dt[,1:2], circle.col=c("turquoise", "salmon"))
image.png
网友评论