github-actions[bot]
commented
4 years ago 表达量矩阵分组很复杂也可以使用limma的3大策略 by 生信技能树
学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！
我把3年前的收费视频课程：3年前的GEO数据挖掘课程你可以听3小时或者3天甚至3个月，免费到B站:
这个课程超级棒，B站免费学习咯：https://m.bilibili.com/video/BV1dy4y1C7jz 配套代码在GitHub哈：https://github.com/jmzeng1314/GSE76275-TNBC TCGA数据库挖掘，代码在：https://github.com/jmzeng1314/TCGA_BRCA GTEx数据库挖掘，代码在：https://github.com/jmzeng1314/gtex_BRCA METABRIC数据库挖掘，代码在：https://github.com/jmzeng1314/METABRIC
然后马上就有了3千多学习量，而且有学员给出来了图文并茂版本万字笔记，让我非常感动！
扫描下面二维码马上就可以学习起来啦，笔记需要至少半个小时来阅读哦！
如何我给大家布置了几十个GSE数据集，让大家练习，熟练掌握我的代码，其中有一些数据集分组会超级复杂。然后优秀学徒根据我的指导提供了一个比较详尽的笔记，分享给大家！
第一次看到这么多分组头都大了。首先要考虑如何分组得到grouplist，其次考虑如何在limma包中分组分析。
听说limma包的官方文档中对这些特殊的情况描述的很细致，于是我找到了这张图，觉得和我目前所面临的情况十分相似
首先下载数据，
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语言读取
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$Subjecttable(g)
于是g变成了这样的
并且很好的各种都排列在一起
获取平台信息
eSet@annotation# Illumina HumanHT-12 V4.0 expression beadchipGPL=eSet@annotation###这里是GPL570##对应的找注释平台和包的网页在http://www.bio-info-trainee.com/1399.htmlif(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=idshead(probes_anno)genes_expr <- filterEM(probes_expr,probes_anno )head(genes_expr)
开始分析
n = genes_expr# 首先是每个组都和第一个组比较，比如时间或者浓度梯度的药物处理library(limma)design=model.matrix(~factor(g))designfit=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)
O我查了下coef的含义，昨天之前我还真没注意这个小东西的作用。之前都是直接处理两个分组，或者从多个分组中选取两个分组分析，昨天处理的数据全是乱七八糟的分组..刚开始直接就做了，3个分组的limma分析也直接用2个分组的套路分析，然后后来某一刻顿悟...发现哦这样不行，于是想起来之前小洁老师讲过一篇GSE474的例子，于是可以完美的借用呀。然后耗费了很多时间，把之前的一一改正。确实有些弯路不走不长记性鸭。
## 然后是每个组都和其它所有样本比较，比如不同的组织gstr(g)table(g)comp1toOther=do.call(cbind,lapply(c('NEC', 'NTC','TC', 'TEC'),function(x) as.numeric( g %in% x)))library(limma)comp1toOtherfit=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老师给的，呃以前就是走流程来分析，现在好歹也多了解了一点点细节，哦统计计算方面的问题不予评价。
## 最后是任意选取两个或者多个分组，组合比较kp=g %in% c('NEC', 'NTC','TC', 'TEC')genes_expr=n[,kp]group_list=g[kp]library(limma)design=model.matrix(~factor(group_list))designfit=lmFit(genes_expr,design)fit=eBayes(fit)###可以调节分组DEG=topTable(fit,coef=2,n=Inf)head(DEG)deg_heatmap(DEG, n, g, topn = 20)
火山图
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)
写到最后
如果你也想开启自己的生物信息学数据处理生涯，加入生信技能树小圈子，还等什么呢，赶快行动起来吧！ 发邮件（jmzeng1314@163.com）给生信技能树创始人jimmy就有惊喜哦！当然了，不能是辣鸡或者骚扰邮件啦，带上自己的简历和想学习交流的诚心吧！
https://mp.weixin.qq.com/s/OIqN_aXx66Y9cGoiVy0wJw