利用RNA-Seq数据构建人类共表达网络13

1、学士学位论文论文题目：利用RNA-Seq数据构建人类共表达网络 作 者： 黄立波导 师： 庞尔丽 副教授系别年级： 生命科学学院 2010级学 号： 201011202950学科专业： 生物科学完成日期： 2014年5月北京师范大学教务处北京师范大学学士学位论文（设计）诚信承诺书本人郑重声明： 所呈交的学士学位论文（设计），是本人在导师的指导下，独立进行研究工作所取得的成果。除文中已经注明引用的内容外，本论文不含任何其他个人或集体已经发表或撰写过的作品成果。对本文的研究做出重要贡献的个人和集体，均已在文中以明确方式标明。本人完全意识到本声明的法律结果由本人承担。本人签名： 年 月 日北京师范大

2、学学士学位论文（设计）使用授权的说明本人完全了解北京师范大学有关收集、保留和使用学士学位论文（设计）的规定，即：本科生在校攻读学位期间论文（设计）工作的知识产权单位属北京师范大学。学校有权保留并向国家有关部门或机构送交论文的复印件和电子版，允许学位论文（设计）被查阅和借阅；学校可以公布学位论文的全部或部分内容，可以采用影印、缩印或扫描等复制手段保存、汇编学位论文。保密的学位论文在解密后遵守此规定。本论文（是、否）保密论文。保密论文在 年解密后适用本授权书。本人签名： 年 月 日 目录摘要IABSTRACTII前言11.背景综述21.1基因共表达网络及构建方法21.1.1基因和基因表达21.1.

3、2复杂网络和基因共表达网络及其构建方法21.2 转录组与RNA-Seq技术51.2.1转录组51.2.2早期研究转录组的基本方法61.2.3全转录组鸟枪法测序(RNA-Seq)61.2.4 RNA-seq的主要用途81.3 共表达网络可视化81.3.1 Cytoscape:网络可视化工具81.3.2 WGCNA生成节点和边的纯文本文件81.3.3将边文件导入Cytoscape生成网络图92.材料与方法102.1 计算环境、操作环境与研究流程简介102.2数据来源112.2.1 人类参考基因组文件112.2.2 人类基因组注释文件112.2.3人类RNA-Seq数据112.3 RNA-Seq数据

4、处理112.3.1 将sra数据转换为fastq数据112.3.2 TopHat：拼接RNA-Seq112.3.3 Cufflinks：组装转录本122.3.4 根据注释文件提取编码蛋白的表达量数据122.3.5根据表达量筛选基因122.4 使用WGCNA构建基因共表达网络132.4.1表达量矩阵导入与数据预处理132.4.2 网络构建与模块检测162.4.3 将基因网络文件导出到Cytoscape203.结果与讨论223.1模块个数及各模块的网络情况介绍223.2根据基因网络文件在Cytoscape中作网络可视化243.2.1 三个模块网络图的可视化243.2.2 单个模块网络图的可视化24

5、4.总结与展望264.1 研究过程中存在的主要问题264.2 基因共表达网络的应用展望26参考文献27附录一 筛选编码蛋白基因的perl代码28附录二 WGCNA中使用的R代码30致谢37正文图表目录图 1复杂网络图例3图 2无尺度网络与随机网络的对比4图 3 RNA-Seq的测序步骤7图 4 研究流程图10图 5 检测离群样本的层次聚类15图 6 为值的选取分析网络拓扑18图 7 所有模块的层次聚类图19图 8 模块检测结果22图 9 通过聚类分析寻找离群样本23图 10 三个模块的网络图（相关系数0.68）24图 11 对单个模块基因网络的可视化（相关系数0.94）25图 12 处于网络中

6、心位置的gene的id25表 1现实中的无尺度网络5表 2三种转录组研究方法比较974利用RNA-Seq数据构建人类共表达网络摘要基因共表达是指基因组中某些基因具有相似表达谱的现象，它们可能受到相似的调控，使其表达量的变化趋势相同。基因共表达网络是一种无尺度网络，该网络中的节点代表基因，基因之间的边是由两两相关的基因的表达量确定的，具有共表达关系的基因处于同一个基因共表达子网络之中。构建基因共表达网络，需要基因的表达量数据。RNA-Seq是基于第二代测序技术的全转录组测序技术，利用RNA-Seq数据能计算出各个基因的表达量。得到表达量数据之后，便可以通过一个名为WGCNA的R软件包构建基因共表

7、达网络。关键词:基因共表达网络，人类，RNA-Seq，WGCNA，R软件Using RNA-Seq Data to Construct Human Co-expression NetworkABSTRACTGene co-expression refers to the phenomenon that genes has similar expression profile in genome, they may be under similar regulation, which makes their expression tends to consistent. Co-expressio

8、n networks are undirected, weighted gene networks, the nodes of such a network correspond to genes in expression profiles, and edges between genes are determined by the pairwise correlations between gene expressions. Genes that has co-expression relationship are in the same sub-network. To construct

9、 gene co-expression network, gene expression data is needed. RNA-Seq is “Whole Transcriptome Shotgun Sequencing” based on the next generation sequencing, by using RNA-Seq data can calculate the expression of each gene. After obtaining expression data, then we can construct gene co-expression network

10、 via an R package which is called WGCNA.KEY WORDS: gene co-expression networks, human, RNA-Seq, WGCNA, R software37前言人类基因组大约有两万多个基因，在单个细胞中表达的基因通常仅有几百到几千个，而且很多基因只在特定组织或发育阶段表达。在这些表达的基因中，有一些基因的表达谱相似，也就是说表达量的变化趋于一致，这些基因很可能受到相同的调控1，这种现象就是基因共表达。基因共表达网络正在越来越多地被用于探索基因的系统级别功能。这种网络的建设从概念上来讲是简单直观的：网络的节点代表基因，如果

11、相关的基因表达谱相似，那么代表基因的节点就被边连接起来，同时可以给边权重2。网络提供了一个系统观察节点之间相互作用的方法，它给我们提供了一个在系统水平上研究基因之间关系的平台。1.背景综述1.1基因共表达网络及构建方法1.1.1基因和基因表达基因是现代生物学的最基本的概念之一。随着生物学的发展，我们对基因的认识也在逐步深入。最初，遗传学的奠基人孟德尔在研究豌豆的过程中发现了孟德尔分离定律和自由组合定律，并提出了“遗传因子”的概念。孟德尔指出生物的每一个性状都是通过“遗传因子”来传递，“遗传因子”是独立的遗传单位，这样就把可观察的遗传性状和控制它的“遗传因子”区分开来。1909年，丹麦遗传学家约

12、翰逊提出“基因”的概念，替代了孟德尔假定的“遗传因子”。1926年，摩尔根在其出版的基因论中提出基因不仅是决定性状的功能单位，而且是一个突变单位和交换单位，这进一步丰富了基因的概念。随着科学技术的发展，生物体的主要遗传物质DNA及其独特的双螺旋结构被发现，基因断裂、基因重叠理论被提出，内含子、外显子、转座子、启动子以及假基因等被陆续发现，所有的相关研究成果都使人们加深了对基因的认识。基因表达指基因在生物体内的转录、剪接、翻译以及转变成具有生物活性的蛋白质分子的所有过程。在DNA被转录翻译的过程中，一种特殊的RNA起到了媒介作用，我们把这种RNA称作mRNA（message RNA，信使RNA）

13、。蛋白质是基因表达的终产物，蛋白质参与了生物体内种类繁多的生命过程，有的蛋白质参与构建生物机体，有的参与催化过程，有的参与调控过程，有的参与运输等等。总而言之，蛋白质是生命活动的主要承担者。因此，研究基因的表达情况，在各类生物学研究中是非常重要的。由于通过测量蛋白质的丰度来研究基因的表达情况需要复杂的技术和昂贵的价格，我们更多的选择通过测量mRNA的丰度来反映基因的表达情况。1.1.2复杂网络和基因共表达网络及其构建方法复杂网络基于网络的方法在许多领域都很有用，如：基因共表达网络，蛋白质相互作用网络，细胞相互作用网络，万维网和社交网络等。事实表明自然界中存在的大量复杂系统均可通过

14、网络加以描述。网络由节点与节点之间的边组成，其中节点表示真实系统中的个体，而边表示个体之间的关系。复杂网络是点与点之间连接关系较复杂一类网络的总称。图 1是一个随机生成的BA（BarabásiAlbert (BA) model，即无尺度网络模型）模型复杂网络。网络的复杂性体现在以下三个方面：a. 结构复杂性：网络的连接结构错综复杂，而且节点之间的连接可能具有不同的权重和方向。b. 节点复杂性：网络中的每个节点都可能具有复杂的演化行为。而且网络中可能存在不同类型的节点，如生化网络中的基质和酶、神经网络中不同的神经元等。c. 各种复杂因素的影响和作用：实际的复杂网络会受到各种复杂因素的影

15、响和作用，例如神经元被同时激活，其连接就加强。目前复杂网络的研究内容主要有以下几个方面：a. 揭示复杂网络结构的统计性质，发现度量这些性质的方法。b. 建立合适的网络模型帮助人们理解网络统计性质的意义和产生机理。c. 基于单个网络节点的性质和整个网络的拓扑结构性质分析与预测网络的行为。d. 提出改善现有网络性能和设计新的网络的有效方法。图 1复杂网络图例基因共表达网络基因共表达网络是无尺度（scale-free）的权重基因网络3。无尺度特性，又称作无标度特性，是指网络的度分布满足幂律分布。所谓一个网络的度分布，是当随机地从网络中抽取一个节点时，与这个节点相连的节点数（叫做这个节点

16、的度）d的概率分布。比如说对一个n个节点组成的完全图（所有节点之间都连有边的图），度分布是：d = n - 1 的概率是1，其余的都是0。无尺度网络的度分布满足幂律分布，也就是说d = k 的概率正比于k 的某个幂次（一般是负的）：Pd=kk- 式 （1-1）幂律分布这一特性，正说明了无尺度网络的度分布与一般随机网络的同。随机网络的度分布属于正态分布，因此有一个特征度数，即大部分节点的度数都接近它。无尺度网络的度分布是呈集散分布：大部分的节点只有比较少的连接，而少数节点有大量的连接。由于不存在特征度数，因此得名“无尺度”。如图 2所示，(a)为随机网络，该网络中大部分节点都连出2到3条边，0条

17、与一条边的和四条边的都很少。(b)为无尺度网络，大部分节点连一条边，少数节点（红色）连有大量边。图 2无尺度网络与随机网络的对比不少现实的网络结构属于无尺度网络，或者有无尺度特性。比如：表 1现实中的无尺度网络网络节点连接电影演员网络演员出演同一部电影因特网路由器物理连接蛋白质相互作用网络蛋白质蛋白质之间的相互作用关系美国飞机航班网络机场飞机航线和这些无尺度网络类似，基因共表达网络中的节点是少数“关键基因”，这些“关键基因”和许多基因的表达谱相似 。基因共表达网络的构建方法在构建基因共表达网络时，我们要借助R4软件及一个名为WGCNA3, 5第三方R包。R软件是一个很方便进行数据

18、处理、计算和图像展示的软件集成套件6。R是当前最受欢迎的数据分析和可视化平台之一，并且它开源、拥有Windows、Mac OS和Linux的版本。WGCNA（weighted gene co-expression network analysis）的全称是加权基因共表达网络分析，该算法是一种构建基因共表达网络的典型系统生物学算法，该算法基于高通量的基因信使RNA（mRNA）表达芯片数据，被广泛应用于生物信息和医学领域。WGCNA算法首先假定基因网络服从无尺度分布，并定义基因共表达相关矩阵、基因网络形成的邻接函数，然后计算不同节点的相异系数，并据此构建分层聚类树（hierarchical clu

19、stering tree），该聚类数的不同分支代表不同的基因模块（module），模块内基因共表达程度高，而分属不同模块的基因共表达程度低7。通过探索模块与特定表型或疾病的关联，可以鉴定引发疾病的靶基因、预测调控通路等。1.2 转录组与RNA-Seq技术 1.2.1转录组转录组是特定细胞或组织在特定时间或状态下转录出来的所有RNA的集合。通过对转录组的研究可以揭示生物体的基因表达、研究结构变异及发现新基因等。转录组分析的研究方法、研究平台发生着日新月异的变化, 同时生物信息学分析的内容也在随之逐渐完善。广义上讲，转录组是生物体细胞或组织在特定状态下所转录出来的所有RNA的总和。RNA包括编码蛋

20、白的RNA（即mRNA）和非编码蛋白的RNA（即ncRNA，如rRNA,tRNA,microRNA等），而mRNA和基因的表达量息息相关，所有狭义上的转录组通常指所有mRNA的总和，而不包含ncRNA。与基因组不同，转录组更具有时间性和空间性。例如，人体大部分细胞具有一模一样的基因，而即使同一细胞在不同生长时期及环境下，其基因表达情况也是不完全相同的。所有，转录组通常反映的是特定条件下活跃表达的基因。1.2.2早期研究转录组的基本方法早期的转录组研究方法主要有：a. 基于杂交技术，如cDNA芯片和寡聚核苷酸芯片；b. 基于第一代测序技术，如基于Sanger测序的SAGE (Serial Ana

21、lysis of Gene Expression)和MPSS(Massively Parallel Signature Sequencing)、全长cDNA文库和EST文库的测序分析。1.2.3全转录组鸟枪法测序(RNA-Seq)RNA测序（RNA Sequencing，简称RNA-Seq，也被称为全转录组鸟枪法测序：Whole Transcriptome Shotgun Sequencing8，WTSS）是基于第二代测序技术的转录组学研究方法：首先提取生物样品的全部转录的RNA，然后反转录为c-DNA后进行二代高通量测序，在此基础上进行片段的重叠 组装，从而可得到一个个的转录本，进而可以形成

22、对该生物样品当前发育状态的基因表达状况的全局了解（globle）。进一步说，若和下一阶段的生物样品的 RNA-Seq转录组进行比较，则可以得到全部的（在转录层面）基因表达的上调及下调-这就形成了表达谱，针对关键基因则可以形成你要想要的 pathway的构建。图 3显示的是RNA-Seq的测序步骤。图 3 RNA-Seq的测序步骤三种主要的转录组研究方法对比，见表 2。其中RNA-Seq主要具有以下优势：a 通量高。运用第二代测序平台可得到几个到几百亿个碱基序列，可以达到覆盖整个基因组和转录组的要求；b 灵敏度高。可以检测细胞中少至几个拷贝的稀有转录本；c 分辨率高。RNA-seq的分辨率能达到

23、单个碱基，准确度好，同时不存在传统微阵列杂交的荧光模拟信号带来的交叉反应和背景噪声；d 不受限制性。可以对任意物种进行全转录组分析，无需预先设计特异性探针，能够直接对任何物种进行转录组分析。同时能够检测未知基因、发现新的转录本，并准确地识别可变剪切位点及SNP、UTR区域。表 2三种转录组研究方法比较9技术芯片SAGE/MPSS/Cdna/ESTRNA-Seq原理杂交Sanger测序高通量测序信号荧光模拟信号数字化信号数字化信号分辨率数个-10bp单碱基单碱基通量高低高背景噪声高低低分析成本高高低起始RNA用量多多少同时映射转录区域和基因表达是有限的基因表达是能够区分不同的亚型有限是是能够区别

24、等位基因有限是是1.2.4 RNA-seq的主要用途RNA-Seq技术能够在单核苷酸水平对特定物种的整体转录活动进行检测，从而全面快速地获得该物种在某一状态下的几乎所有转录本信息。由于转录组测序可以得到全部RNA转录本的丰度信息，加之准确度高，使得它具有非常广泛的应用领域。如：a 检测新的转录本10, 11，包括未知转录本和稀有转录本；b 基因转录水平研究12，如基因表达量、不同样本间差异表达；c 非编码区域功能研究，如microRNA13、非编码长RNA（lncRNA）14、RNA编辑15；d 转录本结构变异研究，如可变剪接、基因融合；e 开发SNPs和SSR等。1.3 共表达网络可视化1.

25、3.1 Cytoscape:网络可视化工具Cytoscape是一个开源软件项目，它使用高通量的表达数据和其他分子状态到一个概念框架来整合生物分子相互作用网络16。生物网络的计算机辅助模型是系统生物学的基石。1.3.2 WGCNA生成节点和边的纯文本文件将网络导入Cytoscape之前，需要将网络数据导入至两个纯文本文件中，便于Cytoscape读取数据。第一个文件包含的是nodes数据，该文件包含三列：第一列是节点名；第二列是节点别名（这是一个可选列，若在WGCNA生成文件过程中提供注释文件，则显示各节点的别名；若不提供注释文件，那么此列全部显示为NA，及无值）；第三列是节点属性，以模块的颜色

26、名为值。第二个文件包含的是edges数据，该文件包含六列：第一列和第二列是有共表达关系的两个基因的基因名；类似的，第五列和第六列是这两个基因别名；第三列是关系的权重（weight）；第四列是网络关系的属性，由于基因共表达网络是无尺度网络，所有该列的值全为undirected。1.3.3将边文件导入Cytoscape生成网络图在上一步中我们获得了两个纯文本文件，我们只需要第二个edges文件就能构建网络图。导入该文件后，在软件的导入设置中，将第一列设置为fromNode，第二列设置为toNode，最后把第三列设为网络关系属性，完成设置，便可生成网络图了。2.材料与方法2.1 计算环境、操作环境与

27、研究流程简介数据处理和分析在北京师范大学计算分子生物学实验室的大型HP计算机集群上完成。HP高性能计算集群系统配置了20个计算节点，其中专配有两台32核心、256G内存的胖节点及50T的存储。在作网络图是，我使用的是HP工作站，CPU为2GHz*8，内存为32G，使用系统为CentOS6.2。在我使用的服务器（本地ip为55）上，我的工作目录为/home/huanglb。本次的研究流程如图 4所示。图 4 研究流程图在图 4中，蓝色矩形图标代表数据，其他颜色的的椭圆形图标代表处理数据的软件。在流程图的最上面是原始的三大类数据，在处理后我们依次得到了基因表达量数据、基因网络

28、数据并最终实现基因网络可视化。2.2数据来源2.2.1 人类参考基因组文件本研究使用的人类基因组数据序列版本号是GRCH37，发布日期是2009年2月，下载于Ensemble(/index.html)。2.2.2 人类基因组注释文件本研究使用的人类基因组文件gencode.v17.gtf下载于GENCODE网（/），版本号是Gencode17，文件格式为gtf，发布日期为2013-06-17。2.2.3人类RNA-Seq数据在NCBI（美国国家生物信息中心）数据库中下载所需的十六个组织的RNA-Se

29、q数据。（/geo/query/acc.cgi?acc=GSE12946）。该数据是Nature杂志上发表的文章Alternative Isoform Regulation in Human Tissue Transcriptomes中使用的实验数据17。该实验数据中包含人类16个组织（heart, breast, testes, BT474, MAQC UHR, MB435, lymph node, T47D, HME, MCF-7, adipose, colon, liver, brain, MAQC human cell lines,

30、skeletal muscle）的RNA-Seq数据。2.3 RNA-Seq数据处理2.3.1 将sra数据转换为fastq数据从NCBI上下载的数据是sra格式的，我们需要把数据的格式转换为fastq格式，我们使用的工具是sratoolkit，软件版本为sratoolkit.2.1.16-centos_linux64 (/Traces/sra/?view=software)。2.3.2 TopHat：拼接RNA-SeqTopHat是一种高效的序列拼接软件，可以将RNA-Seq短片段拼接到人类参考基因组上。我们使用的TopHat(http:/t

31、/)软件版本为tophat-2.0.8b.Linux_x86_64。在使用该软件时，我们使用的参数均为默认参数。2.3.3 Cufflinks：组装转录本Cufflinks能将已拼接到基因组上的RNA-Seq短序列组装起来。我们使用的版本是cufflinks_2.0.2 (/)。Cufflinks将TopHat中比对到人类参考基因组上的序列组装成转录本，并通过计算转录本的相对丰度得出基因的表达量即FPKM(Fragments Per Kilobase of transcript per Million

32、mapped reads，每1百万个拼接上的短序列中拼接到外显子上的每一千个碱基上的短序列数)值。在这一步结束之后，我们就可以得到十六个组织的基因表达量数据。我们把十六个组装的表达量数据整合到一个名为fpkm.txt文件中。2.3.4 根据注释文件提取编码蛋白的表达量数据在我们得到的表达量数据中，除了有编码蛋白的基因，还有假基因和线粒体基因等，而我们需要的是编码蛋白的基因，那么我们编写一个perl脚本（主要代码见附录一，前期shell处理过程略），根据gencode.v17.gtf文件中的gene_type属性来筛选数据，只把编码蛋白基因的表达量数据留下了，并将文件另存为fpkm.new，为t

33、ab分隔符文件。处理好的文件中含有20345个基因，整个文件有20345行，16列。2.3.5根据表达量筛选基因在WGCNA中计算基因相关性的时候不能含零值，否则会产生NaN（not a number，非数），比如除0就会产生NaN。所有我们把含零的基因都删除。最后，剩下15749个基因。最终文件中有15749行，16列。该文件就可以拿来作共表达网络分析和可视化了。2.4 使用WGCNA构建基因共表达网络2.4.1表达量矩阵导入与数据预处理 导入表达量数据在打开R软件（版本为3.03）之后，先设置好工作路径，加载相应的包和数据，并设置好相关参数：#设置工作目录setwd(&quo

34、t;/home/huanglb");#查看工作目录getwd();#加载WGCNA包library(WGCNA);#设置参数，这个参数非常重要options(stringsAsFactors = FALSE);#读取表达量数据。我数据文件名为fpkm.new，是一个tab分隔符文件#该数据是一个16*15749的表达量矩阵，将其导入为R中名为fpkm的数据框中fpkm <- read.table("fpkm.new",header=TRUE,sep="t");#查看数据框的维数和列名dim(fpkm);names(fpkm);在可视化环境

35、中，可以使用fix(fpkm)来查看数据框。表达量矩阵中的16代表十六个组织，即基因表达的十六个样本；表达量矩阵中的15749代表从20345个人类基因中筛选出来的基因。为了后续的计算，需要把导入的表达量矩阵作转置：#将表达量矩阵转置为16*15749矩阵datExpr0 = as.data.frame(t(fpkm, -c(1);names(datExpr0) = fpkm$tracking_id;rownames(datExpr0) = names(fpkm)-c(1);在转置后的矩阵中，每一列为相应基因在不同组织中的表达量。那么我们就得到了15749个基因的表达量列向量。

36、检测数据中有过多缺失值和离群的样本我们先检测有过多缺失值的基因和样本：gsg = goodSamplesGenes(datExpr0, verbose = 3);gsg$allOK如果最后一行代码返回TURE，那么所有的基因都符合要求。如果返回FALSE，那么我们就从基因和样本中移除违规基因：if (!gsg$allOK)#以下打印出被移除的基因和样本名是可选的if (sum(!gsg$goodGenes)>0)printFlush(paste("Removing genes:", paste(names(datExpr0)!gsg$goodGenes, collap

37、se = ", ");if (sum(!gsg$goodSamples)>0)printFlush(paste("Removing samples:", paste(rownames(datExpr0)!gsg$goodSamples, collapse = ", ");#从数据中移除违规基因:datExpr0 = datExpr0gsg$goodSamples, gsg$goodGenes在该过程中，程序对我的数据返回了TRUE，那么就可以跳过以上移除违规基因的代码，直接进入下一步。接下来我们聚类样本来看看是否有明显的离群样本

38、。我们使用提供快速进行层次聚类的flashClust函数。#作样本树图pdf(file = "Plots/fpkm-01.pdf", width = 12, height = 9);par(cex = 0.6);par(mar = c(0,4,2,0)plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,cex.axis = 1.5, cex.main = 2)图 5 检测离群样本

39、的层次聚类该图所显示的样本从左至右依次为：heart, breast, testes, BT474, MAQC UHR, MB435, lymph node, T47D, HME, MCF-7, adipose, colon, liver, brain, MAQC human cell lines, skeletal muscle。从图中我们可以很明显地找出一个离群样本，即heart样本，我们将其移除：如图 5，通过选择一个高度来移除样本，我们选择的是200000，并画一条红线。#在剪切位置画一条红色的线abline(h = 200000,col="red")dev.off

40、()#确定红线之下样本的聚类clust = cutreeStatic(sampleTree, cutHeight = 200000, minSize = 10)table(clust)#返回的结果为：#0 1#1 15#在clust 1 含有我们想要的聚类keepSamples = (clust=1)datExpr = datExpr0keepSamples, nGenes = ncol(datExpr)nSamples = nrow(datExpr)现在变量fpkm中包含的表达量数据可以用于网络分析了。现在我们将变量fpkm储存起来，就可以在下次想使用的时候可以之间载入并调用，不必再次计算。

41、save(datExpr,file="fpkm-01-dataInput.RData")2.4.2 网络构建与模块检测 初步：设置R会话在这一步开始的时候我们打开了一个新的R会话。接下来加载WGCNA包，设置基本的参数并加载上一步保存的数据。接下来的代码会在拥有多核心的计算机上使用多线程计算。这种多线程计算在终端和R软件原生的GUI(Graphical User Interface)上工作良好，但是目前在RStudio和其他的第三方R环境中不支持这种多线程计算。如果你在使用第三方R环境，请在以下代码中不要调用enableWGCNAThreads()函数。#设置

42、工作目录setwd("/home/huanglb");#查看工作目录getwd();#加载WGCNA包library(WGCNA);#设置参数，这个参数非常重要options(stringsAsFactors = FALSE);#允许WGCNA使用多线程计算。这在多核心计算机上会提高计算速度enableWGCNAThreads()#加载第一部分保存的数据lnames = load(file = "fpkm-01-dataInput.RData");#在变量lnames中包含加载的变量的名字lnames 选择邻接矩阵权重参数：分析网络拓扑构建

43、一个权重基因网络需要选择邻接矩阵权重参数，选择一个合适的值可以使网络更符合无尺度网络。在下面的代码中将使用pickSoftThreshold来分析网络拓扑并帮助使用者选择合适的值。# 给出候选的值powers = c(c(1:10), seq(from = 12, to=20, by=2)# 调用网络拓扑分析函数sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)# 对结果作图pdf(file="Plots/fpkm-02-1.pdf",width=9,height=5)par(mfrow =

44、 c(1,2);cex1 = 0.8;plot(sft$fitIndices,1, -sign(sft$fitIndices,3)*sft$fitIndices,2,xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R2",type="n",main = paste("Scale independence");text(sft$fitIndices,1, -sign(sft$fitIndices,3)*sft$fitIn

45、dices,2,labels=powers,cex=cex1,col="red");# 在图中高度为0.8的位置画线，这个值对应的是相关系数的平方R2abline(h=0.8,col="red")plot(sft$fitIndices,1, sft$fitIndices,5,xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",main = paste("Mean connectivity")text

46、(sft$fitIndices,1, sft$fitIndices,5, labels=powers, cex=cex1,col="red")dev.off()计算结果如图 6所示，在左图中，纵轴表示相关系数的平方，该值取值越高，网络越接近无尺度分布。这里我们取该值为0.8，即红线所示位置。从图中可以看出，当power取6的时候即将接近0.8。在接下来的计算中，我们对，即代码中对应的power值取6。右图的纵轴代表对应的基因模块中所有基因邻接函数的均值，即平均连通性。图 6 为值的选取分析网络拓扑 网络构建和模块检测现在，我们只需要调用一个函数就能构建基因网络

47、和识别模块：net = blockwiseModules(datExpr, power = 6,TOMType = "unsigned", minModuleSize = 30,reassignThreshold = 0, mergeCutHeight = 0.25,numericLabels = TRUE, pamRespectsDendro = FALSE,saveTOMs = TRUE,saveTOMFileBase = "fpkmTOM",verbose = 3)以上代码的计算时间稍长。在计算结束之后，我们可以看看程序检测出了多少个模块：结果显示

48、总共有25个模块，编号为1-25，编号正下方对应的数字是每个模块中含有的基因数。其中1号模块最大，有2025个基因，25号模块最小，含有33个模块。编号0中包含的是不属于任何模块的基因，有141个。接下来，我们对模块作层次聚类树状图，并用不同的颜色表示不同的模块。# 将每个模块编号对应上一个模块颜色以便作图pdf(file="Plots/fpkm-02-2.pdf",width=12,height=9)mergedColors = labels2colors(net$colors)# 作层次图并在下方显示相应颜色plotDendroAndColors(net$dendrog

49、rams1, mergedColorsnet$blockGenes1,"Module colors",dendroLabels = FALSE, hang = 0.03,addGuide = TRUE, guideHang = 0.05)dev.off()作图结果如图 7所示。图 7 所有模块的层次聚类图接下来，将这次计算的得到的变量保存，以便下次使用。moduleLabels = net$colorsmoduleColors = labels2colors(net$colors)MEs = net$MEs;geneTree = net$dendrograms1;save(

50、net, MEs, moduleLabels, moduleColors, geneTree,net,file = "fpkm-02-networkConstruction-auto.RData")2.4.3 将基因网络文件导出到CytoscapeCytoscape允许用户导入节点文件和边文件，并设置线条粗细和节点颜色。在WGCNA中，我们一次可以导出一个模块的点和边文件，也可以导出多个模块的点和边文件。在接下来的代码中，我们将所有编号的点和边文件导出到一个文件中。#载入第一部分保存的变量：表达量矩阵lnames = load(file = "fpkm-01-da

51、taInput.RData");#查看载入的变量名lnames#载入第一部分保存的变量：网络数据lnames = load(file = "fpkm-02-networkConstruction-auto.RData");lnames#计算拓扑重叠。这个计算所需时间比较长TOM = TOMsimilarityFromExpr(datExpr, power = 6);#选择模块modules = c("black","blue","brown","cyan","darkgree

52、n","darkgrey","darkred","darkturquoise","green","greenyellow","grey","grey60","lightcyan","lightgreen","lightyellow" ,"magenta","midnightblue","orange","pink&qu

53、ot;,"purple","red","royalblue","salmon","tan","turquoise","yellow");#选择模块探测probes = names(datExpr)inModule = is.finite(match(moduleColors, modules);modProbes = probesinModule;#选择对应的拓扑重叠modTOM = TOMinModule, inModule;dimnames(modTO

54、M) = list(modProbes, modProbes)# 导出能被Cytoscape读入的点和边文件cyt = exportNetworkToCytoscape(modTOM,edgeFile = paste("CytoscapeInput-edges-all-070-", paste(modules, collapse="-"), ".txt", sep=""),nodeFile = paste("CytoscapeInput-nodes-all-070-", paste(module

55、s, collapse="-"), ".txt", sep=""),weighted = TRUE,threshold = 0.70,nodeNames = modProbes,nodeAttr = moduleColorsinModule);#将本次计算的变量保存，以便下次直接使用save(TOM,file="fpkm-06-dataInput.RData")在以上代码中，TOMsimilarityFromExpr函数计算需要较长的时间。需要注意的是，在生成点和边文件的时候，我们选择了拓扑重叠阈值为0.7，该值对

56、应的是相关系数为0.94（0.7=0.94，其中的值在之前确定为6）左右的基因对。代码执行完成之后，将生成所有网络的点和边的txt文件。3.结果与讨论3.1模块个数及各模块的网络情况介绍使用在WGCNA导入和初步处理数据的时候，我们发现心脏数据明显离群，如图 9所示，作聚类分析可以看出它和其他十五个组织分属两类。我们在这一步把heart数据删除了，所以后继的所有计算都是由剩下的十五个人类组织完成的。在第二步网络构建和模块检测中，我们根据分析网络拓扑来确定邻接矩阵权重参数的值为6，据此将基因分成了25个模块（如图 8），编号为0的基因不属于任何模块。在图中我们可以看出模块的大小差异很大，最大的模

57、块包含2025个基因，而最小的模块含33个基因。在设置模块参数时，我们将模块包含的最少基因数设置为30。图 8 模块检测结果图 9 通过聚类分析寻找离群样本3.2根据基因网络文件在Cytoscape中作网络可视化3.2.1 三个模块网络图的可视化在WGCNA作分层聚类时，我们将基因分成25个模块。这25个模块中最大的含有两千多个基因，最小的只有33个基因。理论上我们可以对25个模块做全局网络图，用25种不同的颜色表示不同模块中的基因。但是由于机器内存有限（32G），无法将如此大的数据导入Cytoscape，我们采取了折衷的办法，对23个模块作网络图，如图 10。图 10 三个模块的网络图（相关系数0.68）在该图中，主要有两个模块，另外一个较小的模块未能显示出来。由于网络过密集，我们无法获得更多信息，所以接下来会单独分析一个模块3.2.2 单个模块网络图的可视化我们对图 10中左边的网络单独拿出来实现了可视化，如图 11所