mmu-macrophages-GSE69607-GPL1261 / 开普饭

五月份的学徒专注于GEO数据库里面的表达量芯片数据处理，主要的难点是表达量矩阵获取和探针的基因名字转换，合理的分组后就是标准的差异分析，富集分析。主要是参考我八年前的笔记：

下面是James的投稿

一、数据集介绍

GEO链接：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE69607
芯片平台：GPL1261（https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL1261） [Mouse430_2] Affymetrix Mouse Genome 430 2.0 Array

样品列表：

GSM1704683 WT1M0, biological replicate 1GSM1704684 WT2M0, biological replicate 2GSM1704685 WT3M0, biological replicate 3GSM1704686 WT1M1, biological replicate 1GSM1704687 WT2M1, biological replicate 2GSM1704688 WT3M1, biological replicate 3

文章链接是：Novel Markers to Delineate Murine M1 and M2 Macrophages. PLoS One 2015;10(12):e0145342. PMID: 26699615(https://pubmed.ncbi.nlm.nih.gov/26699615/)

二、核心步骤

1. R包加载

rm(list = ls())library(AnnoProbe)library(GEOquery) library(ggplot2) library(ggstatsplot) library(reshape2)library(patchwork)

2. 获取并且检查表达量矩阵

主要是得是否需要log

# 获取表达量矩阵gse_number <- 'GSE69607'gset <- geoChina(gse_number)a=gset[[1]] dat=exprs(a) dim(dat)

# 检查，判断需不需要取logdat[1:4,1:4] dat=log2(dat)library(limma)dat=normalizeBetweenArrays(dat)

# 画图，使用ggplot需宽数据变长数据class(dat)data <- as.data.frame(dat)data <- melt(data)head(data)title <- paste (gse_number, "/", a@annotation, sep ="")p1 <- ggplot(data,aes(x=variable,y=value))+ geom_boxplot()+ theme_ggstatsplot()+ theme(panel.grid = element_blank(), axis.text=element_text(size=10,face = 'bold'), axis.text.x=element_text(angle=90), plot.title = element_text(hjust = 0.5,size =15))+ xlab('')+ ylab('')+ ggtitle(title)p1#可以看到，处理前后我们的表达量矩阵的表达量范围箱线图如下所示：

3. 根据生物学背景及研究目的人为分组

pd=pData(a) #通过查看说明书知道取对象a里的临床信息用pData## 挑选一些感兴趣的临床表型。head(pd)[,1:4]library(stringr)

##～～～分组信息按需修改～～～# Classically (M1) and alternatively activated (M2) macrophages play distinct roles # in various physiological and disease processes. pd_m1 <- pd[!str_detect(pd$title,'M2'),]

dat_m1 <- dat[,colnames(dat) %in% rownames(pd_m1)]

group_list=ifelse(grepl('M1',pd_m1$title),'M1','M0')

dat <- dat_m1

#以最简单的2分组做展示table(group_list)## group_list## M0 M1 ## 3 3

4. probe_id 和symbol的转换至表达矩阵

获取芯片注释信息

ids=idmap('GPL1261','soft')head(ids)#探针基因ID对应以及去冗余ids=ids[ids$symbol != '',]dat=dat[rownames(dat) %in% ids$ID,]ids=ids[match(rownames(dat),ids$ID),]head(ids) colnames(ids)=c('probe_id','symbol') ids$probe_id=as.character(ids$probe_id)rownames(dat)=ids$probe_iddat[1:4,1:4]

ids=ids[ids$probe_id %in% rownames(dat),]dat[1:4,1:4] dat=dat[ids$probe_id,]

整理芯片注释信息，使探针与基因symbol无重复且一一对应

#接下来，使探针与基因symbol一一对应

ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行

ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids

ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s

dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat

rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名

fivenum(dat['Actb',])fivenum(dat['Gapdh',])

可以看到此芯片的探针与基因ID或者symbol的对应关系如下所示：

dat[1:4,1:4] #保留每个基因ID第一次出现的信息## GSM1704683 GSM1704684 GSM1704685 GSM1704686## Zzz3 8.485861 8.773211 8.675836 8.495483## Zzef1 9.236859 9.187063 9.174501 8.608217## Zyx 10.451070 10.335393 10.254424 9.114873## Zyg11b 8.352753 8.315857 8.421117 8.708668

保存为R数据文件：step1-output.Rdata

save(gse_number,dat,group_list, file = 'step1-output.Rdata')

三、标准步骤之质控

得到标准的3张图，包括主成分分析，高变基因的表达量热图，样品相关性热图

rm(list = ls()) ## 魔幻操作，一键清空~options(stringsAsFactors = F)library(ggplot2)library(ggstatsplot)library(patchwork)library(ggplotify)load(file = 'step1-output.Rdata')table(group_list)#每次都要检测数据dat[1:4,1:4]## 下面是画PCA的必须操作，需要看说明书。exp <- datexp=t(exp)#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换exp=as.data.frame(exp)#将matrix转换为data.frame library("FactoMineR")#画主成分分析图需要加载这两个包library("factoextra")

dat.pca <- PCA(exp , graph = FALSE)#现在exp最后一列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的

# 画图，主成分分析图p2this_title <- paste0(gse_number,'_PCA')p2 <- fviz_pca_ind(dat.pca, geom.ind = "point", # show points only (nbut not "text") col.ind = group_list, # color by groups palette = "Dark2", addEllipses = TRUE, # Concentration ellipses legend.title = "Groups")+ ggtitle(this_title)+ theme_ggstatsplot()+ theme(plot.title = element_text(size=15,hjust = 0.5))p2

# 下面是1000_sd热图library(pheatmap)

cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个

n=t(scale(t(dat[cg,]))) n[n>2]=2 n[n< -2]= -2n[1:4,1:4]ac=data.frame(Group=group_list)rownames(ac)=colnames(n)# 画图，高变基因的表达量热图p3p3 <- pheatmap::pheatmap(n, show_colnames =F, show_rownames = F, main = gse_number, annotation_col=ac, breaks = seq(-3,3,length.out = 100))#因为已经手动设置了表达量最大值，所以，可以不用设置break

# 画图，样品相关性热图p4colD=data.frame(Group=group_list)exprSet=datrownames(colD)=colnames(exprSet)#问题-exprSet设置成转置后的expp4 <- pheatmap::pheatmap(cor(exprSet),#热图对样本-列操作 annotation_col = colD, show_rownames = F, main = gse_number )

四、标准步骤之limma差异分析

library(limma)design=model.matrix(~factor( group_list ))fit=lmFit(dat,design)fit=eBayes(fit)options(digits = 4) #设置全局的数字有效位数为4deg = topTable(fit,coef=2,adjust='BH', n=Inf)

有了差异分析就可以进行标准的可视化，包括火山图和上下调的差异基因热图

nrDEG=deghead(nrDEG)attach(nrDEG)df=nrDEGdf$v= -log10(P.Value) #df新增加一列'v',值为-log10(P.Value)cut_logFC <- with(deg,mean(abs(deg$logFC)) + 2*sd(abs(deg$logFC)) )cut_logFClogFC_t = 1#设定上下调基因df$g=ifelse(df$P.Value>0.05,'stable', ifelse( df$logFC >logFC_t,'up', ifelse( df$logFC < -logFC_t,'down','stable') ))#统计上下调基因数量table(df$g)#给绘制火山图用的数据新增一列symboldf$name=rownames(df)head(df)#设置可循环使用的plot标题this_tile <- paste0('Cutoff for logFC is ',round(logFC_t,3), '\nThe number of up gene is ',nrow(df[df$g == 'up',]) , '\nThe number of down gene is ',nrow(df[df$g == 'down',]))#画图，火山图p5p5 <- ggplot(data = df, aes(x = logFC, y = -log10(P.Value))) + geom_point(alpha=0.6, size=1.5, aes(color=g)) + ylab("-log10(Pvalue)")+ scale_color_manual(values=c("#34bfb5", "#828586","#ff6633"))+ geom_vline(xintercept= 0,lty=4,col="grey",lwd=0.8) + xlim(-3, 3)+ ylim(0,6)+ theme_classic()+ ggtitle(this_tile )+ theme(plot.title = element_text(size=8,hjust = 0.5), legend.title = element_blank(), )#热图library(pheatmap)x=deg$logFC names(x)=rownames(deg) cg=c(names(head(sort(x),100)),

names(tail(sort(x),100)))#对x进行从小到大排列，取前100及后100，并取其对应的探针名，作为向量赋值给cg

n=t(scale(t(dat[cg,])))n[n>2]=2n[n< -2]= -2n[1:4,1:4]ac=data.frame(group=group_list)rownames(ac)=colnames(n) # 画图，上下调的差异基因热图p6p6 <- pheatmap(n,show_colnames =F, show_rownames = F, cluster_cols = T, main = gse_number, annotation_col=ac)

五、标准步骤之生物学功能数据库注释

我们这里不根据任何武断的阈值来区分统计学显著的上下调基因，而是直接根据基因的变化情况排序进行gsea分析，而且仅仅是展示kegg这个生物学功能数据库的注释情况！

gsea分析需要基因的ENTREZID，需要根据物种进行转换

# 加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）

deg$symbol=rownames(deg)library(org.Mm.eg.db)library(clusterProfiler)df <- bitr(deg$symbol, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Mm.eg.db)#人#bitr()用于SYMBOL转ENTREZID#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDbhead(df)DEG=degDEG=merge(DEG,df,by.y='SYMBOL',by.x='symbol')save(DEG,file = 'anno_DEG.Rdata')

进行gsea富集

library(dplyr)library(ggplot2) geneList=DEG$logFCnames(geneList)=DEG$ENTREZIDgeneList=sort(geneList,decreasing = T)head(geneList)library(clusterProfiler)kk_gse <- gseKEGG(geneList = geneList, organism = 'mmu',#按需替换 #nPerm = 1000, minGSSize = 10, pvalueCutoff = 0.9, verbose = FALSE)tmp=kk_gse@resultdim(tmp)kk=DOSE::setReadable(kk_gse, OrgDb='org.Mm.eg.db',keyType='ENTREZID')#按需替换#DOSE::setReadable():mapping geneID to gene Symboltmp=kk@resultdim(tmp)pro='comp1'write.csv(kk@result,paste0(pro,'_kegg.gsea.csv'))save(kk,file = 'gsea_kk.Rdata')

上面的kk这个变量就存储了kegg这个生物学功能数据库的gsea分析结果，我们进行简单可视化，代码如下：

# 展现前6个上调通路和6个下调通路down_k <- kk_gse[tail(order(kk_gse$enrichmentScore,decreasing = F)),];down_k$group=-1up_k <- kk_gse[head(order(kk_gse$enrichmentScore,decreasing = F)),];up_k$group=1

dat=rbind(up_k,down_k)colnames(dat)dat$pvalue = -log10(dat$pvalue)dat$pvalue=dat$pvalue*dat$group dat=dat[order(dat$pvalue,decreasing = F),]# gsea分析结果p7p7<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) + geom_bar(stat="identity") + scale_fill_gradient(low="#34bfb5",high="#ff6633",guide = FALSE) + scale_x_discrete(name ="Pathway names") + scale_y_continuous(name ="log10P-value") + coord_flip() + theme_ggstatsplot()+ theme(plot.title = element_text(size = 15,hjust = 0.5), axis.text = element_text(size = 12,face = 'bold'), panel.grid = element_blank())+ ggtitle("Pathway Enrichment") p7

具体看上面条形图里面的每个通路的gsea分布情况p8

library(enrichplot)p8 <- gseaplot2(kk, geneSetID = rownames(down_k))+ gseaplot2(kk, geneSetID = rownames(up_k))p8

mmu-macrophages-GSE69607-GPL1261

一、数据集介绍

二、核心步骤

1. R包加载

2. 获取并且检查表达量矩阵

3. 根据生物学背景及研究目的人为分组

4. probe_id 和symbol的转换至表达矩阵

三、标准步骤之质控

四、标准步骤之limma差异分析

五、标准步骤之生物学功能数据库注释

如果你也有类似的数据分析需求，却苦于不会写代码，可以考虑找我们的工程师帮忙哦！

转录组产品线

公共数据库产品线

mmu-macrophages-GSE69607-GPL1261

一、数据集介绍

二、核心步骤

1. R包加载

2. 获取并且检查表达量矩阵

3. 根据生物学背景及研究目的人为分组

4. probe_id 和symbol的转换至表达矩阵

三、标准步骤之质控

四、标准步骤之limma差异分析

五、标准步骤之生物学功能数据库注释

如果你也有类似的数据分析需求，却苦于不会写代码，可以考虑找我们的工程师帮忙哦！

转录组产品线

公共数据库产品线

相关推荐