|
样本聚类：flashClust包提供了一种比hclust()较快的聚类方法，比如sampleTree = flashClust(dist(datExpr0), method =
"average")。也可以直接使用plotClusterTreeSamples()直接查看样本聚类情况，如下图，y表示样本的一个分类，比如正常组和对照组。
注意：datExpr0中行代表样本，而列代表基因。
 0. 数据前处理
WGCNA包对芯片数据有特殊的处理习惯，即将芯片数据矩阵的行表示为样本，列表示基因。所以N×M矩阵，表示N个样本，M个基因。正因如此，如果使用函数dist(exprData)可以得到按照样本（即“按行聚类”）的聚类结果。
1. 选择何时“软阀值（soft thresholding power）”beta
使用函数pickSoftThreshold(datExpr,
powerVector = powers, networkType =
'unsigned',
这个函数返回一个列表，第一项是powerEstimate是估计的最优power；第二项是fitIndices是详细的矩阵数据，其中第五列是mean.k表示平均“连接度（connectivity）”。连接度是指某个基因与该网络中其他基因的相关程度，常为“相关度”之和。
2. 构建网络
使用函数blockwiseModules()自动构建不同module。这个函数有大量的参数可供调试，比较重要的有：maxBlockSize默认为5000，表示在这个数值内的基因将整体被计算，如果调大需要更多的内存；power是软阀值；deepSplit是分割灵敏度，取值在0到4（最灵敏）之间，默认是2；minModuleSize设定每个modle的对小容量（即节点数），可以调大；numericLabels默认为FALSE返回颜色，设定为TRUE则返回数字；mergeCutHeight是在计算完所有modules后，将特征量高度相似的modules进行和并，和并的标准就是所有小与mergeCutHeight数值，默认为0.15，可以调大（即减少module）；pamRespectsDendro逻辑值，默认为TRUE，可以设定为FALSE。
返回结果是一个列表，第一个colors表示基因被分为不同的module，0表示没有任何module接受，可以使用table()函数查看，如果返回的是数字，可以使用labels2colors()转换为对应颜色；MEs是一个关于modules的特征量矩阵，行数等于筛选的modules数，列数等于样本数；dendrograms表示在一个block中所有基因的进化树图
也可以分步构建网络。1. 计算每个基因之间的“相关度（adjacency）矩阵”，使用函数adjacency()，返回N×N的矩阵，N表示所有基因数。同时，可以计算每个基因的连接度（即colSums()），之后使用scaleFreePlot()查看连接度，并且判断网络是否属于“无标度拓扑网络(scale
free topology)”。因为无标度网络的连接度有幂律分布，因此scaleFreePlot()会呈现线性。2. 计算TOM(Topological Overlap
Matrix)，使用函数TOM=TOMsimilarity()计算相似矩阵，使用TOMdist()计算不相似矩阵（即1-TOM）。3.
构建基因进化树，使用函数flashClust()，使用plot()函数查看进化树，hang参数调整分枝度，hang越小，树枝越少；进化树的每一个树枝代表一个基因。4.使用dynamicTreeCut包修剪树枝，使用函数cutreeDynamic()，参数与blockwiseModules()相似。也可以使用cutreeStaticColor()修剪树枝。两种方法各有利弊，可以根据生物学意义选择。
5. 和并相似的共表达网络，这个过程相对于把上一步筛选得到的modules再进行聚类，使用函数moduleEigengenes()，确定修剪高度后，使用mergeCloseModules()函数和并modules。
使用函数plotDendroAndColors()查看系统发生树；blockGenes是一个block中所有的基因。dendro表示进化树；而colors表示modules的颜色，如果需要表示不同的modules对比，只需要把colors进行和并cbind()即可。这将有利于展示不同的module划分方法，并且有助于选择modules。
3. 结合外部生物学参数选择模型
通常我们会得到多个module，这时需要结合其他生物学信息选择合适的模型。“基因显著参数（Gene
Significance, GS）”为非负数字信息，比如基于样本的临床数据（clinical
trails）和基于每个基因的
-log(p-value)等。
3.1 结合样本信息，比如样本的重量、长度、程度等信息
样本信息应该是一个非负矩阵，行数等于样本数。具体步骤：1.
计算模型特征矩阵和样本信息矩阵的相关度，使用函数cor()，方法可以选择Pearson。矩阵中的数值范围是(-1,1)，越接近于-1表示负相关，越接近于1表示正相关，说接近于0表示相关性最弱。2.
计算相关度的p值，使用函数corPvalueStudent(cor,
nSamples),p值的越小表示相关性越显著。之后，可以使用labeledHeatmap()将模型和上述相关矩阵可视化，比如下图。从下图可以看出与weight相关的modules中，棕色（相关性最强），其次是红色、肉色和蓝色。所以，如果以weight为指标，我们可以选择这四个modules作为下一步研究对象。
注意：
第一步，计算模型与样本信息相关矩阵时候，模型矩阵与样本信息矩阵的行应该相同。比如，模型矩阵是K×P（K个样本，P个模型），样本信息矩阵则为K×Q（K个样本，Q个样本信息）。最后得到一个P×Q的矩阵。
 在选定module后，还可以看这个module的特征向量与样本信息内，在考虑到每个基因层面上的相关性。步骤如下：1.
计算全表达谱矩阵与特征向量矩阵相关性。2. 计算全基因表达谱与样本信息相关性。3. 使用verboseScatterplot()绘制相关性图，如下。
 还可以将样本信息与module特征向量矩阵结合在一起（直接cbind()在一起即可），使用plotEigengeneNetworks()查看，如下图。也可以控制plotHeatmaps参数和plotDendrograms参数分别输出上下两个图。
 查看module与每个基因GS之间的相关性，可以使用函数plotMEpairs()，或者psych包的pairs.panels()函数，以PerformanceAnalytics包的chart.Correlation()函数，如下图。
   同时，可以绘制一个module的基因表达量与该module特征向量图，联合使用plotME()和barplot()函数。  3.2 结合每个基因信息，比如p值、q值等
同样，也可以结合每个基因信息筛选合适的module。比如，经过多重假设检验计算得到的p值（可以转换为-log(p-value)）、矫正后的q值、在某条生物学通路中（GO分类、细胞位置、转录因子调控等）的基因（0表示存在，1表示即不存在）等等。使用plotModuleSignificance(GS,
modulecolor)函数作出条形图，这个条形图的纵坐标的实质就是，基因在每个模型中对应GS的均值。下图表示了每个基因的p值经过处理后，对应module的GS均值。在这种情况下，我们可以优先考虑研究绿色的module。
 与此同时，还可以使用得到网络根据GS筛选基因，networkScreening()筛选样本信息相关的基因，networkScreeningGS()筛选基因信息相关的基因。
 3.3 研究每个module的生物学意义
在得到每个module后，相当于将基因做了分类（对于每个module，类似于找到的一部分基因）。那么，就可以使用类似的GO、生物学通路等等分析，进而解释每一个module的生物学意义。
4. 可视化网络
4.1 查看TOM图
TOM图是指基因和筛选模型的可视化表示。生成TOM图的步骤：1. 生成全基因不相似TOM矩阵，比如1-TOMsimilarityFromExpr(datExpr, power =
6)，可以把得到的不相关矩阵加幂，这样绘制的TOM图色彩差异会比较明显。2. 使用TOMplot()绘制TOM图，如下图。图上和图左是全基因系统发育树，不同颜色亮带表示不同的module，每一个亮点表示每个基因与其他基因的相关性强弱（越亮表示相关性越强）。
查看生成的网络：可以使用VisANT：http://visant./查看筛选的网络。先生成网络拓扑结构，使用函数TOMsimilarityFromExpr()；之后，筛选特定的module，形成一个M×M的矩阵，表示这个网络中有M个节点（node）；最后，使用exportNetworkToVisANT()函数生成visant的输入类型文件，即choose(M,2)个“边（edge）”。也可以使用exportNetworkToCytoscape
|
|
|
