当下载完GEO数据，并做好数据整理后，得到了行名为SYMBOL，列名为SAMPLE的基因表达矩阵（已去重，log2等），现在就要开始进行差异表达分析。
分组表达矩阵，告诉代码，哪些样本属于那些分组
#懒得下载，自己构建基因表达矩阵
x <- 100#行数
y <- 24#列数
r <- rnorm(x*y)#创建24*100个符合正态分布的数据
dim(r) <- c(x,y)#给数据以行、列维度
class(r)#matrix
dim(r)# 100 24
r[,1:6] <- r[,1:6]+5#因为有些为负数，把他全部变为正数，以下一样
r[,7:12] <- r[,7:12]+5.3
r[,13:18] <- r[,13:18]+6.3
r[,19:24] <- r[,19:24]+7
r[1:4,1:4]
# > r[1:4,1:4]
# sample_1 sample_2 sample_3 sample_4
# gene_1 6.930203 3.445526 6.364362 5.482402
# gene_2 5.676164 4.680503 4.450948 5.608166
# gene_3 3.427881 6.070465 5.540297 4.522752
# gene_4 6.666888 6.002557 5.336065 5.947510
rowname <- c(paste0("gene_",1:nrow(r)))
colname <- c(paste0("sample_",1:ncol(r)))
# > length(rowname)
# [1] 100
# > length(colname)
# [1] 24
colnames(r) <- colname#赋列名
rownames(r) <- rowname#赋行名
r[1:4,1:4]
# > r[1:4,1:4]
# sample_1 sample_2 sample_3 sample_4
# gene_1 6.930203 3.445526 6.364362 5.482402
# gene_2 5.676164 4.680503 4.450948 5.608166
# gene_3 3.427881 6.070465 5.540297 4.522752
# gene_4 6.666888 6.002557 5.336065 5.947510
group_list <- rep(c('normal','tumour_1','tumour_2','tumour_3'),each=rep(6,4))
#group_list <- factor(group_list,levels = c('normal','tumour_1','tumour_2','tumour_3'))
library(limma)
design=model.matrix(~group_list)
design
colnames(design) <- levels(group_list)#给design赋列名
rownames(design) <- colname#给design赋行名
#> design
# normal tumour_1 tumour_2 tumour_3
# sample_4 1 0 0 0
# sample_5 1 0 0 0
# sample_6 1 0 0 0
# sample_10 1 1 0 0
# sample_11 1 1 0 0
# sample_12 1 1 0 0
# sample_16 1 0 1 0
# sample_17 1 0 1 0
# sample_18 1 0 1 0
# sample_22 1 0 0 1
# sample_23 1 0 0 1
# sample_24 1 0 0 1
# attr(,"assign")
# [1] 0 1 1 1
# attr(,"contrasts")
# attr(,"contrasts")$`group_list`
# [1] "contr.treatment"
#design是一个分组对比矩阵，他告诉代码哪个样本是属于哪个分组
#design的列顺序就是group_list这个因子向量的level顺序，这个很重要
fit1=lmFit(r,design)
fit2=eBayes(fit1)
allDiff=topTable(fit2,adjust='fdr',coef=2,number=Inf)#根据自己的需求调参数得到差异表达基因
allDiff[1,]
# > allDiff[1,]
# logFC AveExpr t P.Value adj.P.Val B
# gene_71 2.116048 6.101674 3.697074 0.0002239722 0.02239722 0.3369323
#logFC即logFoldChange，即两组之间的log表达量取平均后进行相除，看两组之间某个基因表达的对比情况
#一般是实验组/对照组，正值代表实验组该基因表达量高于对照组，数值代表高多少倍，负值同理
#但是我们怎么知道是实验组/对照组，还是对照组/实验组呢？
#这跟design的设置有关，design=model.matrix(~group_list)会默认第一列为“截距”，而后面的分组都会跟这个“截距”进行对比，这个“截距”会被默认为分母
#而design的列名的顺序是取决于group_list的level顺序
#这边写得比较乱，mark下
mean(r["gene_71",7:12])-mean(r["gene_71",1:6])#验证一下logFC算法是否正确
# > mean(r["gene_71",7:12])-mean(r["gene_71",1:6])
# [1] 2.116048
mean(r["gene_71",1:24])
# > mean(r["gene_71",1:24])
# [1] 6.101674
#可以看出AveExpr是这24个样本的平均表达量，而不是normal和tumour_1两组的平均表达量
#topTable函数的coef参数，coef=2是指design的第2列，即tumour_1，即把tumour_1与normal进行对比
#同理coef=3是指design的第3列，即tumour_2，即把tumour_2与normal进行对比
#所以，这里coef=2也可以写成coef="tumour_1"
#allDiff=topTable(fit2,adjust='fdr',coef="tumour_1",number=Inf)
那假如我想要tumour_2和tumour_1进行对比呢？
这里有两种方法：
#方法一：
#修改group_list的level的顺序
group_list <- factor(group_list,levels = c('tumour_1','normal','tumour_2','tumour_3'))
design=model.matrix(~group_list)
colnames(design) <- levels(group_list)
rownames(design) <- colname
design
#> design
# tumour_1 normal tumour_2 tumour_3
#6 1 1 0 0
#12 1 0 0 0
#18 1 0 1 0
#24 1 0 0 1
#现在tumour_1变成了“截距”，所以就可以进行比较了
#方法二：
design <- model.matrix(~0+group_list)
colnames(design) <- levels(group_list)
rownames(design) <- colname
design
#> design
# tumour_1 normal tumour_2 tumour_3
#6 0 1 0 0
#12 1 0 0 0
#18 0 0 1 0
#24 0 0 0 1
#这里并没有指定哪一个是“截距”
contrast.matrix <- makeContrasts(contrasts=c('tumour_2-normal',
'tumour_1-normal',
'tumour_2-tumour_1'),levels = design)
#> contrast.matrix
# Contrasts
#Levels tumour_2-normal tumour_1-normal tumour_2-tumour_1
# normal -1 -1 0
# tumour_1 0 1 -1
# tumour_2 1 0 1
# tumour_3 0 0 0
fit=lmFit(r,design)
fit1 <- contrasts.fit(fit, contrast.matrix)
fit2=eBayes(fit1)
allDiffc=topTable(fit2,adjust='fdr',coef=2,number=Inf,p.value=0.05)
#coef=2，即指contrast.matrix的第2列
#所以这个方法是自由度最大的一种
还有配对比较矩阵，太晚了，明天写
网友评论