GEO数据挖掘-第三期-口腔鳞状细胞癌（OSCC）

文章来源：生信技能树
原文：GEO数据挖掘-第三期-口腔鳞状细胞癌（OSCC）
文章标题

A pplication of weighted gene co-expression network analysis to identify key modules and hub genes in oral squamous cell carcinoma tumorigenesis

关键词

口腔鳞状细胞癌

第一步 安装包

## 整理流程中的软件包bioPackages <- c(   "stringi",  # 处理字符串  "GEOquery", # 下载包  "limma",    # 差异分析  "ggfortify", "ggplot2", "pheatmap", "ggstatsplot", "VennDiagram", # 作图  "clusterProfiler", "org.Hs.eg.db",                                # 注释  "AnnotationDbi", "impute", "GO.db", "preprocessCore", "WGCNA"     # WGCNA)## 设置软件包的下载镜像local({  # CRAN的镜像设置  r <- getOption( "repos" );   r[ "CRAN" ] <- "https://mirrors.tuna.tsinghua.edu.cn/CRAN/";   options( repos = r )  # bioconductor的镜像设置  BioC <- getOption( "BioC_mirror" );   BioC[ "BioC_mirror" ] <- "https://mirrors.ustc.edu.cn/bioc/";   options( BioC_mirror = BioC )})## 安装未安装的软件包#source( "https://bioconductor.org/biocLite.R" ) #第一次使用bioconductor需要运行lapply( bioPackages,   function( bioPackage ){    if( !require( bioPackage, character.only = T ) ){      CRANpackages <- available.packages()      if( bioPackage %in% rownames( CRANpackages) ){        install.packages( bioPackage )      }else{        BiocInstaller::biocLite( bioPackage, suppressUpdates = FALSE, ask = FALSE)      }    }  })

第二步 整理数据集和临床表型

## 数据集和注释平台的下载步骤略过，详见第一期options( stringsAsFactors = F )load( './gset.Rdata' )library( "GEOquery" )## 取表达矩阵和样本信息表{  gset = gset[[1]]  exprSet = exprs( gset )  pdata = pData( gset )  chl = length( colnames( pdata ) )  group_list = as.character( pdata[, chl] )}colnames(pdata)tmp = pdata[,10:12]library(stringr)status= str_split(as.character(tmp$characteristics_ch1),':',simplify = T)[,2]age = str_split(as.character(tmp$characteristics_ch1.1),':',simplify = T)[,2]sex = str_split(as.character(tmp$characteristics_ch1.2),':',simplify = T)[,2]tmp$status = statustmp$age = agetmp$sex = sexhclust_table = tmp[,4:6]group_list = tmp$statusload('./ID2gene.Rdata')## 去除没有注释的数据集{  exprSet = exprSet[ rownames(exprSet) %in% ID2gene[ , 1 ], ]  ID2gene = ID2gene[ match(rownames(exprSet), ID2gene[ , 1 ] ), ]}dim( exprSet )dim( ID2gene )tail( sort( table( ID2gene[ , 2 ] ) ), n = 12L )## 相同基因的表达数据取最大值{  MAX = by( exprSet, ID2gene[ , 2 ],               function(x) rownames(x)[ which.max( rowMeans(x) ) ] )  MAX = as.character(MAX)  exprSet = exprSet[ rownames(exprSet) %in% MAX , ]  rownames( exprSet ) = ID2gene[ match( rownames( exprSet ), ID2gene[ , 1 ] ), 2 ]}exprSet = log(exprSet)dim(exprSet)exprSet[1:5,1:5]save( exprSet, group_list, file = 'exprSet_for_WGCNA.Rdata')

第三步 绘制进化树

library("WGCNA")## WGCNAload('./exprSet_for_WGCNA.Rdata')load('./final_exprSet.Rdata')colnames( exprSet ) = paste( group_list, 1:ncol( exprSet ), sep = '_' )dim( exprSet )exprSet[1:5,1:5]datExpr = t(exprSet[order(apply(exprSet,1,mad), decreasing = T)[1:5000],])dim( datExpr )datExpr[1:4,1:4]{  nGenes = ncol(datExpr)  nSamples = nrow(datExpr)  #系统聚类树  datExpr_tree<-hclust(dist(datExpr), method = "average")  par(mar = c(0,5,2,0))  plot(datExpr_tree, main = "Sample clustering", sub="", xlab="",       cex.axis = 0.9, cex.main = 1.2, cex.lab=1, cex = 0.7)  status_colors <- numbers2colors(as.numeric(factor(hclust_table$status)),                                   colors = c("red","white","pink"),signed = FALSE)  age_colors <- numbers2colors(as.numeric(factor(hclust_table$age)),                                colors = c("red","white","pink"),signed = FALSE)  sex_colors <- numbers2colors(as.numeric(factor(hclust_table$sex)),                                colors = c("red","white"),signed = FALSE)  par(mar = c(1,4,3,1),cex=0.8)  plotDendroAndColors(datExpr_tree, cbind(status_colors, age_colors, sex_colors),                      groupLabels = colnames(hclust_table),                      cex.dendroLabels = 0.8,                      marAll = c(1, 4, 3, 1),                      cex.rowText = 0.01,                      main = "Sample dendrogram and trait heatmap")}

image

第四步 确定软阀值

powers = c(c(1:10), seq(from = 12, to=20, by=2))## [1]  1  2  3  4  5  6  7  8  9 10 12 14 16 18 20sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)# Plot the results:# Scale independencepar(mfrow = c(1,2));cex1 = 0.9;plot(sft$fitIndices$Power,     -sign(sft$fitIndices$slope) * sft$fitIndices$SFT.R.sq,     xlab="Soft Threshold (power)",     ylab="Scale Free Topology Model Fit,signed R^2",     type="n",      main = paste("Scale independence"))text(sft$fitIndices$Power,     -sign(sft$fitIndices$slope) * sft$fitIndices$SFT.R.sq,     labels=powers,     cex=cex1,     col="red")abline(h=0.90,col="RED")# Mean connectivityplot(sft$fitIndices$Power,     sft$fitIndices$mean.k.,     xlab="Soft Threshold (power)",     ylab="Mean Connectivity",     type="n",     main = paste("Mean connectivity"))text(sft$fitIndices$Power,     sft$fitIndices$mean.k.,     labels=powers,     cex=cex1,     col="red")best_beta=sft$powerEstimatebest_beta

image

第五步 构建共表达矩阵

net = blockwiseModules(datExpr, power = 8,                       TOMType = "unsigned", minModuleSize = 30,                       reassignThreshold = 0, mergeCutHeight = 0.25,                       numericLabels = TRUE, pamRespectsDendro = FALSE,                       saveTOMs = TRUE,                       saveTOMFileBase = "femaleMouseTOM",                       verbose = 3)## 这个诡异的问题必须记录一下，第一次运行blockwiseModules很顺利，后来就开始报错了，网上给的解决方案是重启电脑，实在不行## 把WGCNA包重装一次，详情看这里https://www.biostars.org/p/305714/table(net$colors)moduleColors <- labels2colors(net$colors)# Recalculate MEs with color labelsMEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes# 计算根据模块特征向量基因计算模块相异度：MEDiss = 1-cor(MEs0);# Cluster module eigengenesMETree = hclust(as.dist(MEDiss), method = "average");# Plot the resultplot(METree,      main = "Clustering of module eigengenes",     xlab = "",      sub = "")# 在聚类图中画出剪切线MEDissThres = 0.25abline(h=MEDissThres, col = "red")merge_modules = mergeCloseModules(datExpr, moduleColors, cutHeight = MEDissThres, verbose = 3)mergedColors = merge_modules$colors;mergedMEs = merge_modules$newMEs;plotDendroAndColors(net$dendrograms[[1]], cbind(moduleColors, mergedColors),                    c("Dynamic Tree Cut", "Merged dynamic"),                    dendroLabels = FALSE, hang = 0.03,                    addGuide = TRUE, guideHang = 0.05)

image image

第六步 TOM图

nGenes = ncol(datExpr)nSamples = nrow(datExpr)geneTree = net$dendrograms[[1]]dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6)nSelect = 400set.seed(10)select = sample(nGenes, size = nSelect)selectTOM = dissTOM[select, select]selectTree = hclust(as.dist(selectTOM), method = "average")selectColors = moduleColors[select]sizeGrWindow(9,9)plotDiss = selectTOM^7diag(plotDiss) = NATOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")

image

第七步 模块与形状的关系

design=model.matrix(~0+ datTraits$subtype)colnames(design)=levels(datTraits$subtype)moduleColors <- labels2colors(net$colors)MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenesMEs = orderMEs(MEs0)moduleTraitCor = cor(MEs, design , use = "p")moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)sizeGrWindow(10,6)textMatrix = paste(signif(moduleTraitCor, 2), "
(",                   signif(moduleTraitPvalue, 1), ")", sep = "");dim(textMatrix) = dim(moduleTraitCor)par(mar = c(6, 8.5, 3, 3));# Display the correlation values within a heatmap plotlabeledHeatmap(Matrix = moduleTraitCor,               xLabels = names(hclust_table[,1:3]),               yLabels = names(MEs),               ySymbols = names(MEs),               colorLabels = FALSE,               colors = greenWhiteRed(50),               textMatrix = textMatrix,               setStdMargins = FALSE,               cex.text = 0.5,               zlim = c(-1,1),               main = paste("Module-trait relationships"))

image

MEs = moduleEigengenes(datExpr, moduleColors)$eigengenesMET = orderMEs(MEs)sizeGrWindow(5,7.5);par(cex = 0.9)plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8,                       xLabelsAngle= 90)# Plot the dendrogramsizeGrWindow(6,6);par(cex = 1.0)## 模块的聚类图plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0),                      plotHeatmaps = FALSE)# Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)par(cex = 1.0)## 性状与模块热图plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2),                      plotDendrograms = FALSE, xLabelsAngle = 90)

image

第八步 感兴趣模块的具体基因分析

# names (colors) of the modulesmodNames = substring(names(MEs), 3)geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));## 算出每个模块跟基因的皮尔森相关系数矩阵## MEs是每个模块在每个样本里面的值## datExpr是每个基因在每个样本的表达量MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));names(geneModuleMembership) = paste("MM", modNames, sep="");names(MMPvalue) = paste("p.MM", modNames, sep="");geneTraitSignificance = as.data.frame(cor(datExpr,use = "p"));GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));module = "turquoise"column = match(module, modNames);moduleGenes = moduleColors==module;sizeGrWindow(7, 7);par(mfrow = c(1,1));verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),                   abs(geneTraitSignificance[moduleGenes, 1]),                   xlab = paste("Module Membership in", module, "module"),                   ylab = "Gene significance for STATUS",                   main = paste("Module membership vs. gene significance
"),                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)module = "brown"column = match(module, modNames);moduleGenes = moduleColors==module;sizeGrWindow(7, 7);par(mfrow = c(1,1));verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),                   abs(geneTraitSignificance[moduleGenes, 1]),                   xlab = paste("Module Membership in", module, "module"),                   ylab = "Gene significance for STATUS",                   main = paste("Module membership vs. gene significance
"),                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

image