QQ登录

只需一步,快速开始

使用微信账号登录

查看: 2181|回复: 2

TCGA数据库癌症WGCNA分析

  [复制链接]

该用户从未签到

4

主题

16

帖子

68

积分

注册会员

Rank: 2

积分
68
发表于 2017-9-14 05:47:05 | 显示全部楼层 |阅读模式
生信自学课堂
首先我是学习了生信自学网提供的课程,得到了表达矩阵,然后用表达矩阵做的WGCNA分析,当然网上有些什么“小工具”,如果你想法论文,还是去官网下载数据吧,又一次一个师兄,答辩被问到,当场蒙圈,答辩直接被问到,你这博士毕业,用的就是这乱七八糟的软件,你就不会去官网下载数据?
有这一次血的教训,告诫大家学习和发论文千万别混淆,有些东西就是上不了台面,科研不是儿戏,还得严谨。
现在可以上分析了
  1. setwd('E:/biowolf/TCGA/TCGA_BRCA')
  2. samples=read.csv('BRCA_matrix.txt',sep = '\t',row.names = 1)
  3. dim(samples)
  4. #[1] 100   3
  5. expro=read.csv('merge_matrix.txt.cv.txt',sep = '\t',row.names = 1)
  6. dim(expro)
  7. #[1] 24991   100
复制代码

数据读取完成,从上述结果可以看出100个样本,有24991个基因,这么多基因全部用来做WGCNA很显然没有必要,我们只要选择一些具有代表性的基因就够了,这里我们采取的方式是选择在100个样本中方差较大的那些基因(意味着在不同样本中变化较大)

继续命令:

  1. m.vars=apply(expro,1,var)
  2. expro.upper=expro[which(m.vars>quantile(m.vars, probs = seq(0, 1, 0.25))[4]),]##选择方差最大的前25%个基因作为后续WGCNA的输入数据集
复制代码

通过上述步骤拿到了6248个基因的表达谱作为WGCNA的输入数据集,进一步的我们需要看看样本之间的差异情况

  1. datExpr=as.data.frame(t(expro.upper));
  2. gsg = goodSamplesGenes(datExpr, verbose = 3);
  3. gsg$allOK
  4. sampleTree = hclust(dist(datExpr), method = 'average')
  5. plot(sampleTree, main = 'Sample clustering to detect outliers'
  6.      , sub='', xlab='')
复制代码

从图中可看出大部分样本表现比较相近,而有两个离群样本,对后续的分析可能造成影响,我们需要将其去掉,共得到98个样本


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复

使用道具 举报

该用户从未签到

4

主题

16

帖子

68

积分

注册会员

Rank: 2

积分
68
 楼主| 发表于 2017-9-14 05:49:15 | 显示全部楼层
生信自学课堂
  1. clust = cutreeStatic(sampleTree, cutHeight = 80000, minSize = 10)
  2. table(clust)
  3. #clust
  4. #0  1
  5. #2 98
  6. keepSamples = (clust==1)
  7. datExpr = datExpr[keepSamples, ]
  8. nGenes = ncol(datExpr)
  9. nSamples = nrow(datExpr)
  10. save(datExpr, file = 'FPKM-01-dataInput.RData')
复制代码
得到最终的数据矩阵之后,我们需要确定软阈值,从代码中可以看出pickSoftThreshold很简单,就两个参数,其他默认即可
  1. powers = c(c(1:10), seq(from = 12, to=20, by=2))
  2. sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
  3. ##画图##
  4. par(mfrow = c(1,2));
  5. cex1 = 0.9;
  6. plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
  7.      xlab='Soft Threshold (power)',ylab='Scale Free Topology Model Fit,signed R^2',type='n',
  8.      main = paste('Scale independence'));
  9. text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
  10.      labels=powers,cex=cex1,col='red');
  11. abline(h=0.90,col='red')
  12. plot(sft$fitIndices[,1], sft$fitIndices[,5],
  13.      xlab='Soft Threshold (power)',ylab='Mean Connectivity', type='n',
  14.      main = paste('Mean connectivity'))
  15. text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col='red')
复制代码

从图中可以看出这个软阈值选择7比较合适,选择软阈值7进行共表达模块挖掘
  1. pow=7
  2. net = blockwiseModules(datExpr, power = pow, maxBlockSize = 7000,
  3.                        TOMType = 'unsigned', minModuleSize = 30,
  4.                        reassignThreshold = 0, mergeCutHeight = 0.25,
  5.                        numericLabels = TRUE, pamRespectsDendro = FALSE,
  6.                        saveTOMs = TRUE,
  7.                        saveTOMFileBase = 'FPKM-TOM',
  8.                        verbose = 3)
  9. table(net$colors)
  10. # open a graphics window
  11. #sizeGrWindow(12, 9)
  12. # Convert labels to colors for plotting
  13. mergedColors = labels2colors(net$colors)
  14. # Plot the dendrogram and the module colors underneath
  15. plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
  16.                     groupLabels = c('Module colors',
  17.                                     'GS.weight'),
  18.                     dendroLabels = FALSE, hang = 0.03,
  19.                     addGuide = TRUE, guideHang = 0.05)
复制代码

从图中可以看出大部分基因在灰色区域,灰色部分一般认为是没有模块接受的,从这里也可以看出其实咱们选择的这些基因并不是特别好

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

该用户从未签到

4

主题

16

帖子

68

积分

注册会员

Rank: 2

积分
68
 楼主| 发表于 2017-9-14 05:51:08 | 显示全部楼层
生信自学课堂
那么做到这一步了基本上共表达模块做完了,每个颜色代表一个共表达模块,统计看看各个模块下的基因个数:

那么得到模块之后下一步该做啥呢,或许很多人到这就不知道如何继续分析了

这里就需要咱们利用这些模块搞事情了,举个例子

如果你是整合的数据(整合lnc与gene),那么同时在某个模块中的基因和lncRNA咱们可以认为是 共表达的,这便是lnc-gene共表达关系的获得途径之一了,进一步你可以根据该模块的基因-lnc-基因之间的关系绘制出共表达网络

今天咱们这里不讲这个,而是跟表型关联,咱们已经拿到了这98个样本的ER、PR、HER2阳性阴性信息,那么进一步的咱们可以看看哪些共表达模块跟ER、PR、HER2阴性最相关,代码如下:

  1. moduleLabelsAutomatic = net$colors
  2. moduleColorsAutomatic = labels2colors(moduleLabelsAutomatic)
  3. moduleColorsFemale = moduleColorsAutomatic
  4. MEs0 = moduleEigengenes(datExpr, moduleColorsFemale)$eigengenes
  5. MEsFemale = orderMEs(MEs0)
  6. samples=samples[match(row.names(datExpr),paste0(gsub('-','.',row.names(samples)),'.01')),]#匹配98个样本数据
  7. trainDt=as.matrix(cbind(ifelse(samples[,1]=='Positive',0,1),#将阴性的样本标记为1
  8.       ifelse(samples[,2]=='Positive',0,1),#将阴性的样本标记为1
  9.       ifelse(samples[,3]=='Positive',0,1),#将阴性的样本标记为1
  10.       ifelse(samples[,1]=='Negative'&samples[,2]=='Negative'&samples[,3]=='Negative',1,0))#将三阴性的样本标记为1
  11. )#得到一个表型的0-1矩阵
  12. modTraitCor = cor(MEsFemale, trainDt, use = 'p')
  13. colnames(MEsFemale)
  14. modTraitP = corPvalueStudent(modTraitCor, nSamples)
  15. textMatrix = paste(signif(modTraitCor, 2), '\n(', signif(modTraitP, 1), ')', sep = '')
  16. dim(textMatrix) = dim(modTraitCor)
  17. labeledHeatmap(Matrix = modTraitCor, xLabels = colnames(trainDt), yLabels = names(MEsFemale),
  18.                ySymbols = colnames(modlues), colorLabels = FALSE, colors = greenWhiteRed(50),
  19.                textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.5, zlim = c(-1,1)
  20.                , main = paste('Module-trait relationships'))
复制代码

最终找到几个共表达网络与三阴性表型最相关的模块。

到了这一步其实还没完,咱们已经拿到了最相关的模块比如上图中的yellow模块,那么yellow模块中的基因,哪些基因最重要,我们该如何获取

首先我们需要计算每个基因与模块的相关关系

  1. modTraitCor = cor(MEsFemale, datExpr, use = 'p')
  2. modTraitP = corPvalueStudent(modTraitCor, nSamples)
  3. corYellow=modTraitCor[which(row.names(modTraitCor)=='MEyellow'),]
  4. head(corYellow[order(-corYellow)])
  5. #RAD51AP1     HDAC2     FOXM1    NCAPD2      TPI1      NOP2
  6. #0.9249354 0.9080872 0.8991733 0.8872607 0.8717050 0.8708449
复制代码

从上诉结果中可以看出RAD51AP1,HDAC2两个基因与yellow相关性最好,也就是说这两个基因是yellow模块的hub-gene

进一步的咱们导出yellow模块的基因共表达关系使用cytoscape进行可视化,代码如下:

  1. TOM = TOMsimilarityFromExpr(datExpr, power = pow);
  2. probes = names(datExpr)
  3. mc='yellow'
  4. mcInds=which(match(moduleColorsAutomatic, gsub('^ME','',mc))==1)
  5. modProbes=probes[mcInds]
  6. modTOM = TOM[mcInds, mcInds];

  7. dimnames(modTOM) = list(modProbes, modProbes)
  8. cyt = exportNetworkToCytoscape(modTOM,
  9.                                  edgeFile = paste('edges-', mc, '.txt', sep=''),
  10.                                  nodeFile = paste('nodes-', mc, '.txt', sep=''),
  11.                                  weighted = TRUE,
  12.                                  threshold = median(modTOM),
  13.                                  nodeNames = modProbes,
  14.                                  #altNodeNames = modGenes,
  15.                                  nodeAttr = moduleColorsAutomatic[mcInds]);
复制代码


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

客服热线
18520221056(微信) 周一至周日:09:00 - 22:00
公司官网:http://www.biowolf.cn

速科生物是一家融生信创新、设计、技术开发、服务为核心的生物公司,生信自学网专注于生信培训周边课程开发和代码设计,坚持为客户打造高品质的精品课程和培训服务。

Powered by 生信自学网 © 2016-2019 江西速科生物

QQ|生信自学论坛 ( 赣ICP备19001400号-1 )

GMT+8, 2019-5-20 08:37 , Processed in 0.186697 second(s), 25 queries .

快速回复 返回顶部 返回列表