表达量矩阵分组很复杂也可以使用limma的3大策略

2020-09-29 20:35:33 浏览数 (1)

论一个很复杂的分组到底该怎么办

第一次看到这么多分组头都大了。首先要考虑如何分组得到grouplist，其次考虑如何在limma包中分组分析。

听说limma包的官方文档中对这些特殊的情况描述的很细致，于是我找到了这张图，觉得和我目前所面临的情况十分相似

首先下载数据，

代码语言：javascript复制
rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
library(AnnoProbe)
library(GEOquery)
gse=geoChina('GSE51401')
eSet=gse[[1]]

probes_expr <- exprs(eSet);dim(probes_expr)
head(probes_expr[,1:4])
boxplot(probes_expr,las=2)
#probes_expr=log2(probes_expr 1)
phenoDat <- pData(eSet)
head(phenoDat[,1:4])

由于我不会交叉着分组...所以直接把网页上的分组信息复制粘贴存为了TXT格式的GSE51401文件，然后使用R语言读取

代码语言：javascript复制
a = read.table(file ='GSE51401')
# 分组
index1=grep('TEC',a$V2)
b=a[index1,]
b$Subject = 'TEC'

index2=grep('NEC',a$V2)
d=a[index2,]
d$Subject = 'NEC'

e = a[c(36,38,41,48,49,54,57,62),]
e$Subject = 'TC'

index4=grep('NTC',a$V2)
f=a[index4,]
f$Subject = 'NTC'

h = rbind(b,d,e,f)

g = h$Subject
table(g)

于是g变成了这样的

并且很好的各种都排列在一起

获取平台信息

代码语言：javascript复制
eSet@annotation
# 	Illumina HumanHT-12 V4.0 expression beadchip
GPL=eSet@annotation
###这里是GPL570
##对应的找注释平台和包的网页在http://www.bio-info-trainee.com/1399.html
if(F){
  if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
  library(hgu133plus2.db)
  ls("package:hgu133plus2.db")
  ids <- toTable(hgu133plus2SYMBOL)
  head(ids)
}else if(F){
  getGEO(GPL)
  ids = data.table::fread(paste0(GPL,".soft"),header = T,skip = "ID",data.table = F)
  ids = ids[c("ID","Gene Symbol"),]
  colnames(ids) = c("probe_id")
}else if(T){
  ids = idmap(GPL,type = "bioc")
}

probes_anno=ids
head(probes_anno) 
genes_expr <- filterEM(probes_expr,probes_anno )
head(genes_expr)

开始分析

代码语言：javascript复制
n = genes_expr
# 首先是每个组都和第一个组比较，比如时间或者浓度梯度的药物处理
library(limma)
design=model.matrix(~factor(g))
design
fit=lmFit(n,design)
fit=eBayes(fit)
DEG=topTable(fit,coef=2,n=Inf)
head(DEG)
deg_heatmap(DEG, n, g, topn = 20)
DEG=topTable(fit,coef=3,n=Inf)
head(DEG) 
deg_heatmap(DEG, n, g, topn = 20)
DEG=topTable(fit,coef=4,n=Inf)
head(DEG) 
deg_heatmap(DEG, n, g, topn = 20)

我查了下coef的含义，昨天之前我还真没注意这个小东西的作用。之前都是直接处理两个分组，或者从多个分组中选取两个分组分析，昨天处理的数据全是乱七八糟的分组..刚开始直接就做了，3个分组的limma分析也直接用2个分组的套路分析，然后后来某一刻顿悟...发现哦这样不行，于是想起来之前小洁老师讲过一篇GSE474的例子，于是可以完美的借用呀。

代码语言：javascript复制
## 然后是每个组都和其它所有样本比较，比如不同的组织
g
str(g)
table(g)
comp1toOther=do.call(cbind,lapply(c('NEC', 'NTC','TC', 'TEC'),function(x) as.numeric( g %in% x)))
library(limma)
comp1toOther 
fit=lmFit(n,comp1toOther)
fit=eBayes(fit)
DEG=topTable(fit,coef=1,n=Inf)
head(DEG)
deg_heatmap(DEG, n, g, topn = 20)
DEG=topTable(fit,coef=3,n=Inf)
head(DEG) 
deg_heatmap(DEG, n, g, topn = 20)

代码都是Jimmy老师给的，呃以前就是走流程来分析，现在好歹也多了解了一点点细节，哦统计计算方面的问题不予评价。

代码语言：javascript复制
## 最后是任意选取两个或者多个分组，组合比较
kp=g %in% c('NEC', 'NTC','TC', 'TEC')
genes_expr=n[,kp]
group_list=g[kp]

library(limma)
design=model.matrix(~factor(group_list))
design
fit=lmFit(genes_expr,design)
fit=eBayes(fit)
###可以调节分组
DEG=topTable(fit,coef=2,n=Inf)
head(DEG)
deg_heatmap(DEG, n, g, topn = 20)

火山图

代码语言：javascript复制
need_deg=data.frame(symbols=rownames(DEG), logFC=DEG$logFC, p=DEG$P.Value)
deg_volcano(need_deg,1)
deg_volcano(need_deg,2)
deg_volcano(need_deg, style = 1, p_thred = 0.05, logFC_thred = 1)

txt

0 人点赞