lncRNA实战项目-第六步-WGCNA相关性分析

WGCNA将lncRNA分成18个模块(3635个lncRNA),空间模块中lncRNA表达呈现明显的组织区域特异性,如:CB (M1, 794个lncRNAs),DG/CA1 (M2, 443个lncRNAs), CA1 (M4, 369个lncRNAs),neocortex (M7, 123个lncRNAs)和OC (M10,57个lncRNAs)。时间模块中lncRNA表达与年龄有关,而与组织区域不明显;性别模块中lncRNA表达与性别和年龄都相关。每个模块就必须做pathway/go等数据库的注释分析咯!


资料收集:

google搜索或在生信技能树和生信菜鸟团搜索WGCNA ,能找到很多教程,下面列出几个中文教程和英文教程,强烈推荐中文教程1和英文教程3。

  1. 一文学会WGCNA分析
  2. 学习WGCNA总结
  3. Tutorials for the WGCNA package
  • WGCNA Background and glossary
  • Data input and cleaning
  • Network construction and module detection
  • Relating modules to external information and identifying important genes
  • Interfacing network analysis with other data such as functional annotation and gene ontology
  • Network visualization using WGCNA functions
  • Exporting a gene network to external visualization software

背景知识

基本概念:

WGCNA(weighted correlation network analysis)加权基因共表达网络分析, 用于提取与性状或临床特征相关的基因模块,解析基础代谢途径,转录调控途径、翻译水平调控等生物学过程。WGCNA适合于复杂的数据模式,推荐5组以上的数据,如: 不同器官、组织类型发育调控; 同一组织不同时期发育调控; 非生物胁迫不同时间点应答; 病原物侵染后不同时间点应答。

基本步骤:

WGCNA分为表达量聚类分析和表型关联两部分,具体步骤包括基因之间相关系数的计算,共表达网络的构建,筛选特定模块,模块与性状关联,核心基因的筛选。

Overview of a typical WGCNA analysis.

术语:

Co-expression weighted network: 是一个无向有权重(undirected, weighted)的网络。“无权重(unweighted network)”,基因与基因之间的相关度只能是0或者1,0表示两个基因没有联系,而1表示有。“有权重(weighted network)”基因之间不仅仅是相关与否,还记录着它们的相关性数值,数值就是基因之间的联系的权重(相关性)。 Module:(模块)指表达模式相似的基因聚为一类,这样的一类基因称为模块。 Connectivity(连接度):连接度是指某个基因与该网络中其他基因的相关程度,常为“相关度”之和。 Eigengene(eigen- +‎ gene):基因和样本构成的矩阵。 Hub gene: 模块中的高度相关的核心基因。 Gene signicance ,GS:基因显著性参数,为非负数字信息,比如基于样本的临床数据(clinical trails)和基于每个基因的-log(p-value)等。


分析流程

WGCNA输入文件需要一个表达矩阵,最好是RPKM或其他归一化好的表达量,还需要一个矩阵关于临床信息或者其它关于样本属性的信息。

STEP1: 输入数据的准备

表达矩阵可以从作者GitHub上下载(https://github.com/DChenABLife/RhesusLncRNA),我这里只尝试对lncRNA的表达矩阵(https://github.com/DChenABLife/RhesusLncRNA/blob/master/All_sample_LncRNA_exp_RPKM.xlsx)做后续分析。下载的表达矩阵文件是Excel格式的,需要转为csv格式方便后续用R处理,可以直接打开这个excel文件,然后另存为csv格式即可。

读入原始表达数据

原始数据包含64个样本,9904个lncRNA表达量,这时的矩阵行为基因,列为样本信息,其中第一列是lncRNA ID,(这里的lncRNA ID是cufflinks 组装时给的自由的ID, 是需要和已有的ID对应, 对于新的转录本再通过nr/nt等数据库注释), 第66列是作者给出的注释(我查了几个注释有的也查不出来是什么意思)。

setwd("G:/My_exercise/WGCNA/")
lncRNAexpr <- read.csv("All_sample_LncRNA_exp_RPKM.csv",sep=",",header = T)
head(lncRNAexpr)
dim(lncRNAexpr)
#[1] 9904   66

rawdata

重命名数据列表,行名和列名

##去掉Annotation这列
fpkm <- lncRNAexpr[,-66]
head(fpkm)
##重命名行名和列名
rownames(fpkm)=fpkm[,1]
fpkm=fpkm[,-1]
fpkm[1:4,1:4]

fpkm

准备表型信息

这里有64个样本,包含猕猴脑不同空间区域,不同发育时期,以及性别,因为每个样本都交叉包含着三种不同的信息,如果选择全部表型信息,我试了试,后续的模块和性状完全看不清关系,所以我这里仅选择脑不同区域的表型信息,包括CB、DG、PFC、PCC、CA1、OC、PC、TC。

##Sample Info
subname=sapply(colnames(fpkm),function(x) strsplit(x,"_")[[1]][1])
datTraits = data.frame(gsm=names(fpkm),
                       subtype=subname)
rownames(datTraits)=datTraits[,1]
head(datTraits)

sample.Info

下载并载入WGCNA包

source("http://bioconductor.org/biocLite.R")
#biocLite(c("AnnotationDbi", "impute", "GO.db", "preprocessCore")) 
##如果上面的包已经下载过了,就不用下载
biocLite("WGCNA")
library(WGCNA)

行列转置

WGCNA针对的是基因进行聚类,而一般我们的聚类是针对样本用hclust即可,也就是说要转置为行表示样本, 列表示基因。

RNAseq_voom <- fpkm 
WGCNA_matrix = t(RNAseq_voom[order(apply(RNAseq_voom,1,mad), decreasing = T)[1:5000],])
datExpr <- WGCNA_matrix  ## top 5000 mad genes
datExpr[1:4,1:4]

datExpr

确定临床表型与样本名字

sampleNames = rownames(datExpr);
traitRows = match(sampleNames, datTraits$gsm)  
rownames(datTraits) = datTraits[traitRows, 1]

datExpr和datTraits准备好后,接下来就是构建基因网络,鉴定模块。网络构建有三种方法:1)一步法构建网络;2)多步法构建网络;3)block-wise构建网络(主要针对大数据集)。下面的介绍的步骤是一步法构建网络。

STEP2:确定最佳soft-thresholding power

选择合适的“软阀值(soft thresholding power)”beta, 用到的函数是pickSoftThreshold,pickSoftThreshold(datExpr, powerVector = powers, verbose = 5),powerVector可以是一系列数值,从而选择最优值。这个函数返回一个列表,第一项是powerEstimate是估计的最优power;第二项是fitIndices是详细的矩阵数据,其中第五列是mean.k表示平均“连接度(connectivity)”。

Constructing a weighted gene network entails the choice of the soft thresholding power to which co-expression similarity is raised to calculate adjacency.

# Choose a set of soft-thresholding powers
powers = c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
# Plot the results:
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.90,col="green")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
# 选择beta值
best_beta=sft$powerEstimate
#> best_beta
[1] 3

最佳beta值是3。

Soft Threshold


STEP3: 一步法构建共表达矩阵

一步法构建网络,使用函数blockwiseModules(), 这个函数包含很多参数,其中power=sft$powerEstimate=3 即上一步得到的最佳软阈值;maxBlockSize 默认为5000, 表示在这个数值内的基因将整体被计算,如果调大需要更多的内存;numerricLabels 默认为返回颜色,设置为TRUE则返回数字;mergeCutHeight是合并模块阈值的一个参数。

net = blockwiseModules(datExpr, power = sft$powerEstimate,
                       maxBlockSize = 6000,TOMType = "unsigned", 
                       minModuleSize = 30,reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "AS-green-FPKM-TOM",
                       verbose = 3)

STEP4:模块鉴定及可视化

模块鉴定

上一步的返回结果是一个列表,可以用table()函数查看,0表示没有任何module接受。table(net$colors) 可以看总共有多少模块,每个模块的大小,这里共有9个模块,从1-9每个模块的大小是递减的,从2254-115,0表示这些基因不在所有模块内。

table(net$colors)

可视化

dendrograms表示在一个block中所有基因的进化树图, 使用函数plotDendroAndColors()查看系统发生树;blockGenes是一个block中所有的基因。

# Convert labels to colors for plotting
mergedColors = labels2colors(net$colors)
table(mergedColors)
# Plot the dendrogram and the module colors underneath
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
## assign all of the gene to their corresponding module 
## hclust for the genes.

cluster Dendrogram

#明确样本数和基因数
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.lab = 2, 
     cex.axis = 1, cex.main = 1,cex.lab=1)
## 如果这个时候样本是有性状,或者临床表型的,可以加进去看看是否聚类合理
#针对前面构造的样品矩阵添加对应颜色
sample_colors <- numbers2colors(as.numeric(factor(datTraits$subtype)), 
                                colors = c("grey","blue","red","green"),signed = FALSE)
## 这个给样品添加对应颜色的代码需要自行修改以适应自己的数据分析项目。
#  sample_colors <- numbers2colors( datTraits ,signed = FALSE)
## 如果样品有多种分类情况,而且 datTraits 里面都是分类信息,那么可以直接用上面代码,
##当然,这样给的颜色不明显,意义不大。
#构造10个样品的系统聚类树及性状热图
par(mar = c(1,4,3,1),cex=0.8)
plotDendroAndColors(datExpr_tree, sample_colors,
                    groupLabels = colnames(sample),
                    cex.dendroLabels = 0.8,
                    marAll = c(1, 4, 3, 1),
                    cex.rowText = 0.01,
                    main = "Sample dendrogram and trait heatmap")

Sample dendrogram and trait heatmap

STEP5:模块和性状的关系

MEs是一个关于modules的特征量矩阵,行数等于筛选的modules数,列数等于样本数;

design=model.matrix(~0+ datTraits$subtype)
colnames(design)=levels(datTraits$subtype)
moduleColors <- labels2colors(net$colors)
# Recalculate MEs with color labels
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0); ##不同颜色的模块的ME值矩阵(样本vs模块)
moduleTraitCor = cor(MEs, design , use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

sizeGrWindow(10,6)
# Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(design),
               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"))

图中第二列第五行,即CB/turquoise相关性最强(0.97),pvalue=1e-41,后续分析可以挑选这个模块。 每一列对应的样本特征可以从design这里查看。

STEP6:感兴趣性状的模块的具体基因分析

选定CB/turquoise这个模块后,接下来看看模块与样本特征和基因的相关性,1. 首先计算全模块与基因的相关性矩阵,2. 再计算性状与基因的相关性矩阵,3. 使用verboseScatterplot()绘制相关性图:

首先计算模块与基因的相关性矩阵
# names (colors) of the modules
modNames = 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="");
再计算性状与基因的相关性矩阵
## 只有连续型性状才能只有计算
## 这里把是否属于 CB 表型这个变量用0,1进行数值化。
CB = as.data.frame(design[,2]);
names(CB) = "CB"
geneTraitSignificance = as.data.frame(cor(datExpr, CB, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(CB), sep="");
names(GSPvalue) = paste("p.GS.", names(CB), sep="")
最后把两个相关性矩阵联合起来,指定感兴趣模块进行分析
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 CB",
                   main = paste("Module membership vs. gene significance\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

上图可以看出基因跟其对应的性状高度相关,可以导出做个GO/KEGG注释,看看这些基因的具体功能。

STEP7:网络的可视化

Topological Overlap Matrix (TOM)图是指基因和筛选模型的可视化表示。WGCNA认为基因之间的简单的相关性不足以计算共表达,所以它利用邻近矩阵,又计算了一个新的邻近矩阵。一般来说,TOM就是WGCNA分析的最终结果,后续的只是对TOM的下游注释。

生成TOM图的步骤:1. 生成全基因不相似TOM矩阵,比如dissTOM=1-TOMsimilarityFromExpr(datExpr, power = 6),可以把得到的不相关矩阵加幂,这样绘制的TOM图色彩差异会比较明显。2. 使用TOMplot()绘制TOM图。

#首先针对所有基因画热图
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
geneTree = net$dendrograms[[1]]; 
#生成全基因不相似TOM矩阵
dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6); 
plotTOM = dissTOM^7; 
diag(plotTOM) = NA; 
#TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")

#然后随机选取部分基因作图
nSelect = 400
# For reproducibility, we set the random seed
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
# There’s no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select];
# Open a graphical window
sizeGrWindow(9,9)
# Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
# the color palette; setting the diagonal to NA also improves the clarity of the plot
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")

图上和图左是全基因系统发育树,不同颜色亮带表示不同的module,每一个亮点表示每个基因与其他基因的相关性强弱(越亮表示相关性越强)。

#最后画模块和性状的关系
# Recalculate module eigengenes
MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
## 只有连续型性状才能只有计算
## 这里把是否属于 Luminal 表型这个变量用0,1进行数值化。
CB = as.data.frame(design[,2]);
names(CB) = "CB"
# Add the weight to existing module eigengenes
MET = orderMEs(cbind(MEs, CB))
# Plot the relationships among the eigengenes and the trait
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 dendrogram
sizeGrWindow(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)

STEP8:提取指定模块的基因名

提取基因信息,可以做GO/KEGG等分析,进而解释这些module的生物学意义。这里的lncRNA ID转换着有点麻烦,这一步先略过,之后再看看。

# Select module
module = "turquoise";
# Select module probes
probes = colnames(datExpr) ## 我们例子里面的probe就是基因名
inModule = (moduleColors==module);
modProbes = probes[inModule];

STEP9: 模块的导出

主要模块里面的基因直接的相互作用关系信息可以导出到cytoscape,VisANT等网络可视化软件。

# Recalculate topological overlap
TOM = TOMsimilarityFromExpr(datExpr, power = 6); 
# Select module
module = "turquoise";
# Select module probes
probes = colnames(datExpr) ## 我们例子里面的probe就是基因名
inModule = (moduleColors==module);
modProbes = probes[inModule]; 
## 也是提取指定模块的基因名
# Select the corresponding Topological Overlap
modTOM = TOM[inModule, inModule];
dimnames(modTOM) = list(modProbes, modProbes)
首先是导出到VisANT
vis = exportNetworkToVisANT(modTOM,
                            file = paste("VisANTInput-", module, ".txt", sep=""),
                            weighted = TRUE,
                            threshold = 0)

也可以输出模块里的部分基因:

nTop = 30;
IMConn = softConnectivity(datExpr[, modProbes]);
top = (rank(-IMConn) <= nTop)
vis = exportNetworkToVisANT(modTOM[top, top],
       file = paste("VisANTInput-", module, "-top30.txt", sep=""),
       weighted = TRUE,
       threshold = 0,
       probeToGene = data.frame(annot$substanceBXH, annot$gene_symbol) )
然后是导出到cytoscape
cyt = exportNetworkToCytoscape(
  modTOM,
  edgeFile = paste("CytoscapeInput-edges-", paste(module, collapse="-"), ".txt", sep=""),
  nodeFile = paste("CytoscapeInput-nodes-", paste(module, collapse="-"), ".txt", sep=""),
  weighted = TRUE,
  threshold = 0.02,
  nodeNames = modProbes, 
  nodeAttr = moduleColors[inModule]
  );

编辑:jimmy

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2018-02-17

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏数说工作室

文本分析 | 常用距离/相似度 一览

这个系列打算以文本相似度为切入点,逐步介绍一些文本分析的干货,包括分词、词频、词频向量、TF-IDF、文本匹配等等。 第一篇中,介绍了文本相似度是干什么的; 第...

3654
来自专栏ATYUN订阅号

【教程】使用TensorFlow对象检测接口标注数据集

当为机器学习对象检测和识别模型构建数据集时,为数据集中的所有图像生成标注非常耗时。而这些标注是训练和测试模型所必需的,并且标注必须是准确的。因此,数据集中的所有...

4467
来自专栏深度学习思考者

DL开源框架Caffe | 模型微调 (finetune)的场景、问题、技巧以及解决方案

前言 什么是模型的微调?   使用别人训练好的网络模型进行训练,前提是必须和别人用同一个网络,因为参数是根据网络而来的。当然最后一层是可以修改的,因为我们...

4406
来自专栏人工智能

机器学习入门之HelloWorld(Tensorflow)

本文介绍一些机器学习的入门知识,从安装环境到跑通机器学习入门程序MNIST demo。

2.2K37
来自专栏人人都是极客

OpenCV和SVM分类器在自动驾驶中的车辆检测

这次文章的车辆检测在车辆感知模块中是非常重要的功能,本节课我们的目标如下: 在标记的图像训练集上进行面向梯度的直方图(HOG)特征提取并训练分类器线性SVM分类...

5376
来自专栏PaddlePaddle

【图像分类】使用经典模型进行图像分类

场景文字识别 图像相比文字能够提供更加生动、容易理解及更具艺术感的信息,是人们转递与交换信息的重要来源。图像分类是根据图像的语义信息对不同类别图像进行区分,是计...

7715
来自专栏专知

【专知-Java Deeplearning4j深度学习教程04】使用CNN进行文本分类:图文+代码

【导读】主题链路知识是我们专知的核心功能之一,为用户提供AI领域系统性的知识学习服务,一站式学习人工智能的知识,包含人工智能( 机器学习、自然语言处理、计算机视...

1.5K6
来自专栏机器学习算法与Python学习

Torch7搭建卷积神经网络详细教程

(如果有好的建议和问题欢迎在留言区指出) 之前的博文,如一文读懂卷积神经网络(CNN)、多层网络与反向传播算法详解、感知机详解、卷积神经网络详解等已经比较详细的...

35114
来自专栏专知

【专知国庆特刊-PyTorch手把手深度学习教程系列01】一文带你入门优雅的PyTorch

【导读】主题链路知识是我们专知的核心功能之一,为用户提供AI领域系统性的知识学习服务,一站式学习人工智能的知识,包含人工智能( 机器学习、自然语言处理、计算机视...

7037
来自专栏AI研习社

利用卷积自编码器对图片进行降噪

前言 这周工作太忙,本来想更把 Attention tranlsation 写出来,但一直抽不出时间,等后面有时间再来写。我们这周来看一个简单的自编码器实战...

3554

扫码关注云+社区