下载文档原格式
零基础⼩⽩STATA数据分析实⽤常见命令整理STATA基础⼊门零基础实⽤命令整理第⼀章数据的读⼊与熟悉1.读⼊⽂件中的部分变量. use[变量] using [⽂件名]Eg . use age sex height weight using [⽂件名]2.读⼊⽂件中的部分观察量. use[⽂件名] in X/Y. use "I:\stata\chapter3.dta" in 601/1000软件只读⼊从第601个观察到第1000个观察之间的400个观察量3.描述、管理数据的基本命令命令功能. describe描述数据的基本情况：样本总量、变量总数、变量的格式等. list. list [变量名]－列出数据中所有变量的分布，从第⼀个样本到最后⼀个样本－列出选定变量的分布. list [变量名] in X/Y 列出数据中被选定的变量分布。
in限定数据的观察值范围。
⽐如，若只想查看第100个-200个观察值的分布，则将X/Y替换成100/200. order [变量名]按选定变量排序。
⽐如，样本的编号、年龄、性别、教育程度，……，等. aorder 将所有变量从 a-z 排序. label variable给变量贴上标签命令功能. sort [变量名] －将某个变量的数值进⾏排序。
⼀般情况下，排序的⽅式是从⼩到⼤－可同时排序多个变量－Stata将缺失值描述为最⼤数值，故排列在最后. sort [变量名] [in] 对某些变量的某个取值范围进⾏排序；没有指定的取值范围保持在原地⽅. gsort [+|-][变量名] －可从⼩到⼤和从⼤到⼩－若变量名前没有任何符号或加上+号，则按升序排列；若在变量名前加上－号，则按降序排列－变量可以是数值型、也可以是字符型. gsort [+|-][变量名] ，mfirst －mfirst指定将缺失值置于所有有效数值之前. gsort -age第⼆章变量的⽣成与处理1.离散和连续测量离散⽅式（discrete measure）：由定性测量和定序测量组成；适⽤于低层次数据连续⽅式（continuous measure）：由定距测量和定⽐测量组成。
一文看懂WGCNA分析（2019更新版）发现我这个4年前的WGCNA分析教程可以排在自己最受欢迎的前10个教程里面了，而且直接以我这个授课代码出的SCI文章就有38篇了，当然不排除很多学员使用我的代码却不告知我，也不会致谢我。
不过，我这点战绩根本就算不上什么，其实这个WGCNA包已经是十多年前发表的了，仍然是广受好评及引用量一直在增加，破万也是指日可待。
大家首先可以看到3个教程：•2016-WGCNA-HCC-hub-gene.pdf 中文文章范例）•WGCNA_GBMTutorialHorvath.pdf•WGCNA_YeastTutorialHorvath.pdf其中第一个是我4年前的WGCNA分析教程最主要的参考文献，后面两个是英文教程，我相信你大概率是不会去看的，不过，我还是放在这里了。
（还是需要强调，这两个英文教程完整的展现了WGCNA的全部用法）然后你只需要简单浏览本文档，就可以在rstudio里面打开后缀是proj的文件，打开R代码，一步步跟着店铺！基本概念WGCNA其译为加权基因共表达网络分析。
该分析方法旨在寻找协同表达的基因模块(module)，并探索基因网络与关注的表型之间的关联关系，以及网络中的核心基因。
适用于复杂的数据模式，推荐5组(或者15个样品)以上的数据。
一般可应用的研究方向有：不同器官或组织类型发育调控、同一组织不同发育调控、非生物胁迫不同时间点应答、病原菌侵染后不同时间点应答。
基本原理从方法上来讲，WGCNA分为表达量聚类分析和表型关联两部分，主要包括基因之间相关系数计算、基因模块的确定、共表达网络、模块与性状关联四个步骤。
第一步计算任意两个基因之间的相关系数（Person Coefficient）。
为了衡量两个基因是否具有相似表达模式，一般需要设置阈值来筛选，高于阈值的则认为是相似的。
但是这样如果将阈值设为0.8，那么很难说明0.8和0.79两个是有显著差别的。
使用TASSEL学习GWAS笔记（1-6）完整版使用TASSEL学习GWAS笔记（6/6）：TASSEL结果可视化：QQ plot，曼哈顿图 #2021.9.04笔记计划分为六篇：第一篇：读取plink基因型数据和表型数据第二篇：对基因型数据质控：缺失质控，maf质控，hwe质控，样本质控第三篇：基因型数据可视化：kingship，MDS，PCA第四篇：一般线性模型进行GWAS分析（GLM模型）第五篇：混合线性模型进行GWAS分析（MLM模型）第六篇：TASSEL结果可视化：QQ plot，曼哈顿图已完成前五篇，本篇是第六篇。
1. TASSEL的GLM和MLM分析结果质控后的plink数据和表型数据：「GLM的GWAS分析结果：」「MLM的GWAS分析结果：」2. TASSEL中的可视化TASSEL有对结果进行可视化的模块，包括qq图和曼哈顿图，但是图不方便调整。
这里用TASSEL的分析结果，使用R语言进行绘制qq图和曼哈顿图。
3. R语言包安装及载入需要用到：o qqmano tidyverseo data.table下面代码，会判断是否有这三个包，如果没有，就自动安装。
然后载入软件包。
if(!require(data.table)) install.packages("data.table")if(!require(qqman)) install.packages("qqman")if(!require(tidyverse)) install.packages("tidyverse")library(qqman)library(tidyverse)library(data.table)4. GLM模型GWAS结果可视化results_log = fread("glm-result.txt")dim(results_log)head(results_log)select = dplyr::selecttable(results_log$Trait)结果：> table(results_log$Trait)dpoll EarDia EarHT2460 2460 2460数据中共有三个性状，可以选择一个性状，进行可视化。
TCGA的乳腺癌RNA-seq数据WGCNA分析示例WGCNA（WeightedCorrelationNetworkanalyi）是一个基于基因表达数据，构建基因共表达网络的方法。
WGCNA和差异基因分析（DEG）的差异在于DEG主要分析样本和样本之间的差异，而WGCNA主要分析的是基因和基因之间的关系。
WGCNA通过分析基因之间的关联关系，将基因区分为多个模块。
而最后通过这些模块和样本表型之间的关联性分析，寻找特定表型的分子特征。
网上例子千千万，但是大部分都是从文档翻译而来，要用起来还是有些费劲，要深入的可以移步这里：WGCNA##############etwd('E:/rawData/TCGA_DATA/TCGA-BRCA')ample=read.cv('ClinicalFull_matri某.t某t',ep='\\t',=1)dim(ample)#[1]1003e某pro=read.cv('Merge_matri某.t某t.cv.t某t',ep='\\t',=1)dim(e某pro)#[1]24991100数据读取完成，从上述结果可以看出100个样本，有24991个基因，这么多基因全部用来做WGCNA很显然没有必要，我们只要选择一些具有代表性的基因就够了，这里我们采取的方式是选择在100个样本中方差较大的那些基因（意味着在不同样本中变化较大）继续命令：m.var=apply(e某pro,1,var)e某pro.upper=e某pro[which(m.var>quantile(m.var,prob=eq(0,1,0.25))[4]),]##选择方差最大的前25%个基因作为后续WGCNA的输入数据集通过上述步骤拿到了6248个基因的表达谱作为WGCNA的输入数据集，进一步的我们需要看看样本之间的差异情况datE某pr=a.data.frame(t(e某pro.upper));gg=goodSampleGene(datE某pr,verboe=3);gg$allOK ampleTree=hclut(dit(datE某pr),method='average')plot(ampleTree,main='Samplecluteringtodetec toutlier',ub='',某lab='')从图中可看出大部分样本表现比较相近，而有两个离群样本，对后续的分析可能造成影响，我们需要将其去掉，共得到98个样本clut=cutreeStatic(ampleTree,cutHeight=80000,minSize=10)table(clut )#clut#01#298keepSample=(clut==1)datE某pr=datE某pr[keepSample,]nGene=ncol(datE某pr)nSample=nrow(datE某pr) ave(datE某pr,file='FPKM-01-dataInput.RData')得到最终的数据矩阵之后，我们需要确定软阈值，从代码中可以看出pickSoftThrehold 很简单，就两个参数，其他默认即可power=c(c(1:10),eq(from=12,to=20,by=2))ft=pickSoftThrehold(datE某pr,powerVector=power,verboe=5)##画图##par(mfrow=c(1,2));ce某1=0.9;plot(ft$fitIndice[,1],-ign(ft$fitIndice[,3])某ft$fitIndice[,2],某lab='SoftThrehold(power)',ylab='ScaleFreeTopologyModelFit,ignedR ^2',type='n',main=pate('Scaleindependence'));te某t(ft$fitIndice[,1],-ign(ft$fitIndice[,3])某ft$fitIndice[,2],label=power,ce某=ce某1,col='red');abline(h=0.90,col='red')plot(ft$fitIndice[,1],ft$fitIndice[,5],某lab='SoftThrehold(power)',ylab='MeanConnectivity',type='n', main=pate('Meanconnectivity'))te某t(ft$fitIndice[,1],ft$fitIndice[,5],label=power,ce某=ce某1,col='red')从图中可以看出这个软阈值选择7比较合适,选择软阈值7进行共表达模块挖掘pow=7net=blockwieModule(datE某pr,power=pow,ma某BlockSize=7000, TOMType='unigned',minModuleSize=30,reaignThrehold=0,mergeCutHeight=0.25,numericLabel=TRUE,pamRepectDendro=FALSE,aveTOM=TRUE,aveTOMFileBae='FPKM-TOM',verboe=3)table(net$color)#openagraphicwindow#izeGrWindow(12 ,9)#ConvertlabeltocolorforplottingmergedColor=label2color(net$color)#PlotthedendrogramandthemodulecolorunderneathplotDendroAndCo lor(net$dendrogram[[1]],mergedColor[net$blockGene[[1]]], groupLabel=c('Modulecolor','GS.weight'),dendroLabel=FALSE,ha ng=0.03,addGuide=TRUE,guideHang=0.05)从图中可以看出大部分基因在灰色区域，灰色部分一般认为是没有模块接受的，从这里也可以看出其实咱们选择的这些基因并不是特别好那么做到这一步了基本上共表达模块做完了，每个颜色代表一个共表达模块，统计看看各个模块下的基因个数：那么得到模块之后下一步该做啥呢，或许很多人到这就不知道如何继续分析了这里就需要咱们利用这些模块搞事情了，举个例子如果你是整合的数据（整合lnc与gene），那么同时在某个模块中的基因和lncRNA咱们可以认为是共表达的，这便是lnc-gene共表达关系的获得途径之一了，进一步你可以根据该模块的基因-lnc-基因之间的关系绘制出共表达网络今天咱们这里不讲这个，而是跟表型关联，咱们已经拿到了这98个样本的ER、PR、HER2阳性阴性信息，那么进一步的咱们可以看看哪些共表达模块跟ER、PR、HER2阴性最相关，代码如下：moduleLabelAutomatic= net$colormoduleColorAutomatic=label2color(moduleLabelAutomatic)moduleColorFemale=moduleCol orAutomaticME0=moduleEigengene(datE某pr,moduleColorFemale)$eigengeneMEFemale=orderME(ME0)ample=ample[match((datE某pr),pate0(gub('-','.',(ample)),'.01')),]#匹配98个样本数据trainDt=a.matri某(cbind(ifele(ample[,1]=='Poitive',0,1),#将阴性的样本标记为1ifele(ample[,2]=='Poitive',0,1),#将阴性的样本标记为1ifele(ample[,3]=='Poitive',0,1),#将阴性的样本标记为1ifele(ample[,1]=='Negative'&le[,2]=='Negative'&le[,3]= ='Negative',1,0))#将三阴性的样本标记为1)#得到一个表型的0-1矩阵modTraitCor=cor(MEFemale,trainDt,ue='p')colname(MEFemale)modTraitP=corPvalueStudent(modTraitCor,nSample)te某tMatri某=pate(ignif(modTraitCor,2),'\\n(',ignif(modTraitP,1),')',ep='')d im(te某tMatri某)=dim(modTraitCor)labeledHeatmap(Matri某=modTraitCor,某Label=colname(trainDt),yLabel=name(MEFemale),ySymbol=colname(modlue),colorLabel=FALSE,color=greenWhiteRed (50),te某tMatri某=te某tMatri某,etStdMargin=FALSE,ce某.te某t=0.5,zlim=c(-1,1),main=pate('Module-traitrelationhip'))最终找到几个共表达网络与三阴性表型最相关的模块。
最新思路——巧用WGCNA分析GEO和TCGA数据，文章轻松上5分癌症相关成纤维细胞(CAF)是胃癌基质中最重要的细胞成分，它促进胃癌进展、治疗抵抗和免疫抑制。
加权基因共表达网络分析(Weighted gene co-expression network analysis, WGCNA)作为一种系统的生物信息学算法，能够将高度协调表达的基因整合到几个基因模块中，并研究该模块与目标表型的关系。
但是到目前为止，CAF 和间质浸润尚未在胃癌中进行WGCNA分析。
本次介绍的文章首次将WGCNA同时用于两个转录组数据集（来自GEO和TCGA数据库），用以分析CAF与胃癌的关系。
这篇文章是2021年10月发表在frontiers in Molecular Biosciences上，文章题目是：Weighted Gene Co-expression Network Analysis Identifies a Cancer-Associated Fibroblast Signature for Predicting Prognosis and Therapeutic Responses in Gastric Cancer。
下面我们来具体看看文章的主要内容吧。
材料与方法01文章的方法我们通过以下流程图来看看吧。
作者首先分别从GEO 和TCGA数据库获得431和330个胃癌样本信息。
随后对这两个数据集的样本进行CAF浸润以及基质评分的计算。
之后通过WGCNA构建针对CAF浸润和基质评分的共表达网络。
分别筛选出相关性最高的黑色和棕色module为关键模块，之后通过阈值的设定进一步筛选得到这两个模块中（取交集）的37个重点基因作为hub基因。
为了进一步缩小基因数目，采用单因素Cox和LASSO回归分析将基因数目缩小到4个，并利用这四个基因构建风险模型。
最后针对这个风险模型分别在两个数据集中作了预后分析、化学和免疫治疗预测性能测试、肿瘤突变负荷（TMB）分析、GSEA分析、与CFA浸润以及相关marker 基因相关性分析、以及在分子和蛋白水平对这四个关键基因的验证。
刷爆高分文章的WGCNA究竟是个啥？随着测序价格下降，各位大牛们手头积累的测序数据越来越多，传统的两两比对分析不仅带来了复杂的比对组合，而且无法系统地反馈我们各样本基因之间的相互作用模式。
这时候我们就不得不另辟蹊径，Weighted Gene Co-Expression Network Analysis(以下简称WGCNA)就是一个适合复杂样本的分析方法。
WGCNA中文名译作加权关联网络分析，小R觉得这个翻译有点生硬，还是英文来得比较直接。
它是一种从测序数据中挖掘模块(module)信息的算法。
在该方法中module被定义为一组具有相似表达谱的基因，如果某些基因在一个生理过程或不同组织中总是具有相类似的表达变化，那么我们有理由认为这些基因在功能上是相关的，可以把他们定义为一个模块(module)。
这似乎与聚类分析所得到的结果有那么一点相似，但不同的是，WGCNA的聚类准则具有生物学意义，它是对基因间表达量的相关系数取n次幂，使得相关系数数值的分布逐渐符合无尺度分布，可以将基因按照表达模式进行分类，将模式相似的基因归为一个模块（module），而非常规的聚类方法，因此该方法所得出的结果具有更高的可信度。
当基因module被定义出来后，我们可以利用这些结果做很多进一步的工作。
在co-expression network中，每一个基因在一个特定时间或空间的表达情况被视做一个点(node)，为了得到基因间的关联情况，我们需要计算任何两个基因间的相关系数(Person Coefficient)，第i个基因和第j个基因的Person Coefficient，即两个基因的表达相似性。
为了知道两个基因的表达谱是否具有相似性，需要人为规定一个阈值，只有当基因间的Person Coefficient达到这一阈值后(如0.8)我们才认为这两个基因是相似的，否则则不相似。
但是这种分析方法存在一个很明显的局限，即我们没有理由认为Person Coefficient为0.8的两个基因与Coefficient为0.79的两个基因是有显著差别的，但是以上算法却无法避免这一处境，WGCNA采用了一种基于软阈值的判定方法很好地避免了这一问题。
构建生物网络实践WGCNA使用（五）上期介绍了TCGA的下载，本节开始讲解如何构建共表达网络，首先我们介绍下WGCNA（Weighted Correlation Network analysis），它是一种从芯片数据中挖掘模块(module)信息的算法。
如果在一个生理过程或不同组织中某些基因总是具有相类似的表达变化，那么我们有理由认为这些表达类似的基因在功能上是相关的，可以把它们定义为一个模块(module)。
当然除了寻找有生物功能的模块外，我们还可以通过WGCNA实现性状关联，代谢通路建模，基因共表达网络构建等分析。
这里主要讲解如何构建共表达网络（co-expression network）及其生物学意义。
首先，我们了解下什么是共表达网络。
在共表达网络中，每一个基因在一个特定时间或空间的表达情况被视做一个点(node)。
如果我们做了100张芯片，每张芯片上有2万个基因，那么我们可以用一个100*20000的矩阵文本来表示实验结果。
为了得到基因间的关联情况，我们需要计算任何两个基因间的相关系数(一般采用Person Coefficient)，经过运算后，我们得到一个20000*20000的关系矩阵S，sij表示第i个基因和第j个基因的Person Coefficient，即两个基因的表达相关性。
为了知道两个基因的表达是否具有相关性，需要人为规定一个阈值，只有当基因间的Person Coefficient达到这一阈值后(如+0.7，正为正相关，负为负相关)。
我们才认为这两个基因是相关，否则不相关。
对所有基因进行两两比较运算，这样就又得到一个20000*20000的邻接矩阵（由0或1构成，1表示相关，0表示不相关）。
通过这个矩阵就可以构成网络。
这是一般构建共表达网络的方法，实现以上流程，只需调用R语言内置的cor()函数即可。
但是以上分析方法存在一个很明显的局限性，即我们没有理由认为Person Coefficient为0.7的两个基因与Person Coefficient为0.69的两个基因是有显著差别的。
使用WGCNA构建权共表达网络样本聚类：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', verbose = 5)。
其中，dataExpr的行列要求与上相同；powerVector可以是一系列数值，从而选择最优值；networkType可以选择unsigned（无向？）或者signed（有向？）；默认Rsquared阀值是0.85。
这个函数返回一个列表，第一项是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。
WGCNA新手入门笔记2（含代码和数据）上次我们介绍了WGCNA的入门（WGCNA新手入门笔记（含代码和数据）），大家在安装WGCNA包的时候，可能会遇到GO.db这个包安装不了的问题。
主要问题应该是出在电脑的防火墙，安装时请关闭防火墙。
如果还有问题，请先单独安装AnnotationDbi这个包，biocLite("AnnotationDbi")再安装GO.db，并尝试从本地文件安装该包。
如果还有问题，请使用管理员身份运行R语言，尝试上述步骤。
另外如果大家问题解决了请在留言处留个言，告知大家是在哪一步解决了问题，谢谢！因为本人没有进行单因素实验，不知道到底是哪个因素改变了实验结果。
今天给大家过一遍代码。
网盘中有代码和数据。
链接：/s/1bpvu9Dt密码：w7g4##导入数据##library(WGCNA)options(stringsAsFactors =FALSE)enableWGCNAThreads()enableWGCNAThreads()是指允许R语言程序最大线程运行，像我这电脑是4核CPU的，那么就能用上3核：当然如果当前电脑没别的事，也可以满负荷运作samples=read.csv( 'Sam_info.txt',sep = 't',s =1)expro=read.csv( 'ExpData.txt',sep = 't',s =1)dim(expro)这部分代码是为了让R语言读取外部数据。
当然了在读取数据之前首先改变一下工作目录，这一点在周二的文章中提过了。
R语言读取外部数据的方式常用的有read.table和read.csv，这里用的是read.csv，想要查看某一函数的具体参数，可以用？函数名查看，比如：大家可以注意到read.table和read.csv中header参数的默认值是不同的，header=true表示第一行是标题，第二行才是数据，header=false则表示第一行就是数据，没有标题。
##筛选方差前25%的基因##m.vars=apply(expro,1,var)expro.upper=expro[which(m.vars>quantile(m.vars, probs = seq( 0, 1, 0.25))[ 4]),]dim(expro.upper)datExpr= as.data.frame(t(expro.upper));nGenes =ncol(datExpr)nSamples = nrow(datExpr)这一步是为了减少运算量，因为一个测序数据可能会有好几万个探针，而可能其中很多基因在各个样本中的表达情况并没有什么太大变化，为了减少运算量，这里我们筛选方差前25%的基因。
当然了，以下几种情况，你可以忽略此步，转而执行下列代码datExpr= as.data.frame(t(expro));nGenes =ncol(datExpr)nSamples = nrow(datExpr)情况1：所用的数据是比较老的芯片数据，探针数量较少；情况2：你的电脑足够强大，不必减少运算；情况3：你第一步导入的数据是差异基因的数据，已经作过初筛，比如这一文章的套路（PMID：27133569）（个人不建议这么做）。
##样本聚类检查离群值##gsg = goodSamplesGenes(datExpr, verbose =3);gsg$allOKsampleTree = hclust(dist(datExpr), method = "average")plot(sampleTree, main = "Sample clustering to detect outliers", sub= "", xlab= "")save(datExpr, file = "FPKM-01-dataInput.RData")执行这一段代码我们会得到下面这张图：从结果上来，我们分析的样本没啥离群值，所以代码里说不作处理。
在网上的一个案例中，离群的样本就比较明显了。
（https://www.shengxin.ren/article/88）如果需要去除离群样本，则执行下列代码，其中cutHeight=多少就看你自己了。
clust = cutreeStatic(sampleTree, cutHeight = 20000, minSize = 10)table(clust)keepSamples = (clust==1)datExpr = datExpr[keepSamples, ]nGenes =ncol(datExpr)nSamples = nrow(datExpr)save(datExpr, file = "FPKM-01-dataInput.RData")执行上述代码的话，就会去掉8个样本##软阈值筛选##powers = c(c( 1: 10), seq( from= 12, to= 20, by= 2))sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)par(mfrow = c( 1, 2));cex1 =0.9;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");abline(h= 0.90,col="red")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")软阈值是WGCNA的算法中非常重要的一个环节，简单的说硬阈值是一种一刀切的算法，比如高考分数>500分能上一本，低于500就不行，软阈值的话切起来比较柔和一些，会考虑你学校怎么样，平时成绩怎么样之类。
网盘中的数据跑起来其实是不太好的，没有合适的软阈值，这根线是要划到0.9的。
或者右边这张图趋近于0像下面这个就是比较正常的。
（/2535.html）这一步是为了算powers的值。
一般来说，powers取红线（0.9）左右的数字都是可以的。
如果你天秤座特征比较明显，你也可以运行下列代码，让程序推荐你一个值（本案例中返回值是NA，所以后面为了让程序能够进行下去，选了powers=14）。
sft$powerEstimate这个powers的值在后面的代码中会一直用到，所以你在跑别的数据的时候一定要更改powers的数值。
##一步法网络构建：One-step network construction and module detection##net = blockwiseModules(datExpr, power = 14, maxBlockSize = 6000, TOMType = "unsigned", minModuleSize = 30, reassignThreshold = 0, mergeCutHeight = 0.25, numericLabels = TRUE, pamRespectsDendro = FALSE, saveTOMs = TRUE, saveTOMFileBase = "AS-green-FPKM-TOM", verbose = 3)table(net$colors)这一步就比较烧电脑了，关于网络构建的方法，分为一步法和多步法，一步法比较无脑，多步法能设置更多参数。
想要了解多步法的请参考此文http://tiramisutes.github.io/2016/09/14/WGCNA.html继续说上面的代码，结果跑出来如下图结果是每个模块中包含的基因数量。
一般来说，结果包含十几个到二十几个模块是比较正常的，此外一个模块中的基因数量不宜过多，像我们这个结果里模块1的基因数量达到了2651，这个就有点太多了，主要是因为我们powers=14，软阈值太低了导致的。
所以说上述软阈值的筛选可以对我们的模块分析起到微调的作用。
##绘画结果展示### open a graphics window#sizeGrWindow(12, 9)# Convert labels to colors for plottingmergedColors = labels2colors(net$colors)# 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)由于我们的软阈值比较低，所以这一结果中几乎没有grey模块，grey模块中的基因是共表达分析时没有被接受的基因，可以理解为一群散兵游勇。
当然如果分析结果中grey模块中的基因数量比较多也是不太好的，表示样本中的基因共表达趋势不明显，不同特征的样本之间差异性不大，或者组内基因表达一致性比较差。
