目录

1 背景 BackgroundTypes of biological networks

Motivation for using co-expression networks

Network inference and reverse engineering

Basic graph terminology and data structures

Steps for building a co-expression network

Optimizing parameters for network construction

2 共表达流程 TutorialPreparing RNA-seq data for network construction

Building a co-expression network

Detecting co-expression modules

Annotating a co-expression network

Visualizing network正文开始

1.1 Types of biological networks

生物网络可以包含不同的数据类型，用点(node)和边(edge)区分。常见的网络类型：蛋白互作(PPI)

表示蛋白之间物理联系，它们几乎占据了细胞生物过程的中心位置。蛋白作为点，用无向的线连接

代谢网络

主要表示生化反应，有助于生物生长、繁殖、维持结构。点是代谢产物，并用有向的箭头表示代谢过程或特定反应的调节作用

基因互作

不同的点表示不同基因，描述它们功能相关性；可以根据基因的背景知识来推断线的方向

基因/转录调控

表示基因表达是如何被调控的；点是基因或转录因子，它们之间的关系也是定向，例如Reactome、KEGG等数据库中表示基因调节的关系

细胞信号

点表示通路中的物质，如蛋白、核酸或其他代谢物

各种通路

1.2 Motivation for using co-expression networks

为何研究共表达

看图说话：某个细胞受到刺激1，也许它的A通路就会上调表达，B通路下调，结果可能比刺激前还要理想；

受到另一种刺激2后，A通路下调，B通路上调，那么可能就比较糟糕

通过共表达网络，就可以探索A、B通路是如何被调控的，以及背后基因的相互关系；另外，互作的基因一般都参与同样的生物途径

一般来讲，探索基因表达数据的标准流程是这样：Differnentail expression analysis

想了解转录水平不同处理的差别背后的原因

Gene enrichment analysis (GO/KEGG)

得到差异基因，但是它们是干嘛的不知道，需要注释一下

但是有个弊端，它只能两两比较(如：感染与未感染)，然后得到的结果也只是知道哪些上调哪些下调，是一个宏观的结论

使用Co-expression network 共表达网络可以分析多个处理的基因表达数据(例如：不同时间段处理)，还能推断未知基因产物的功能、检测sub-groups

1.3 Network inference and reverse engineering

利用网络进行推断：可以使用表达量数据、已知的转录因子、ChIP-ChIP或ChIP-seq、时间序列等，因为网络是有向、交叉 的，所以可以判断许多的关系信息

说到网络，就要看一下有向和无向网络：install.packages("igraph")

library(igraph)

set.seed(1)# 先构建一个有向网络临近矩阵# create a 5x5 adjacency matrixadj

plot(g)#再构建一个无向网络adj[upper.tri(adj)]

plot(g)#再构建加权网络 weighted networkset.seed(1)

adj

adj[upper.tri(adj, diag = TRUE)]

plot(g, edge.width = E(g)$weight)

三种网络

构建共表达网络的关键步骤：数据预处理 Data pre-processing数据转换、过滤或者均一化都会影响下游的分析

选择数据: 如果想了解整体调控关系，可以选择全部样本，因为它们中间会存在不同的扰动，更能体现共表达；

也可以选择和某种表型相关的样本

过滤低表达基因 Filter low count genes

过滤低相关性/非差异基因 Filter low-variance/ non-DE genes ：这是构建Robust性能更强的网络的关键一步

数据转换 Log2-CPM：将RNA数据转换成更像芯片数据的类型

标准化 Normalization

临近矩阵构建 Adjacency matrix constructionstep1: 得到相似性矩阵

既然目标是构建共表达网络，那么就要找到这个“共”所在，即找到基因表达量的相似性，因此需要用到一些寻找基因间相似性的算法

基因共表达分析中最常用的两种相关系数

皮尔森 Pearson correlation，是线性相关系数，反映两个变量线性相关程度。如果两个基因的表达量呈线性关系(数学上，线性相关指的是直线相关，指数、幂函数、正弦函数等曲线相关不属于线性相关)，那么两个基因表达量的就有显著的皮尔森相关系性。Pearson相关系数适用条件为两个变量间有线性关系、变量是连续变量、变量均符合正态分布。同一量纲数据可以选择Pearson，例如mRNA表达量数据

但在生物体内的许多调控关系，例如转录因子、小干扰RNA与靶基因，可能都是非线性关系，于是斯皮尔曼系数登场

斯皮尔曼 Spearman correlation ，是针对不同量纲计算的，比如两个通路看着相似，但其实单位不同，无法用pearson直接统计。无论两个变量的数据如何变化，符合什么样的分布，只关心每个数值在变量内的排列顺序。如果两个变量的对应值，在各组内的排序顺位是相同或类似的，则具有显著的相关性。

相关性|r|表明两变量间相关的程度，r为正表示正相关，为负表示负相关

step2 ：根据相关性分组

step3 ： 将“假相关“去除方法一：利用sigmoid 转换：1/(1+e-x )

方法二：利用power转换xn方法一：不管正负，按照相关性绝对值分组

方法二：只将正相关的基因聚在一起

检测共表达模块 Network module detection利用聚类树(反映数据的关系强弱)

聚类树检测

纵坐标反映了相关性，数值越低越相似，其实也可以看作是相关性的反义指标，因此蓝色的那一坨就表现特别相关下面主要介绍前期数据处理部分，数据处理好了，后面真正的WGCNA才能做好

2.1 Preparing RNA-seq data for network construction

2.2 Prepare data and package拿到一个数据，第一件事就是配置包library('gplots')

library('ggplot2')

library('knitr')

library('limma')

library('reshape2')

library('RColorBrewer')

library('WGCNA')

然后读取表型和表达矩阵samples

raw_counts

dim(raw_counts)

大概的样子是这样：

两个必备数据

加载注释信息library('Homo.sapiens') #加载物种注释，这里以人为例keytypes(Homo.sapiens) #查看注释包的内容，里面包含了各种的ID以及基因位置信息# 比如要匹配Ensembl的原始ID到gene ID，并且显示染色体信息以及基因位置gene_ids

select(Homo.sapiens, keytype='ENSEMBL', keys=gene_ids,

columns=c('ALIAS', 'TXCHROM', 'TXSTART', 'TXEND'))# 就会得到类似下面的数据## 'select()' returned 1:many mapping between keys and columns##            ENSEMBL    ALIAS TXCHROM  TXSTART    TXEND## 1  ENSG00000121410      A1B   chr19 58858172 58864865## 2  ENSG00000121410      A1B   chr19 58859832 58874214## 3  ENSG00000121410      ABG   chr19 58858172 58864865## 4  ENSG00000121410      ABG   chr19 58859832 58874214## 5  ENSG00000121410      GAB   chr19 58858172 58864865

Sample check

分析之前，看看样本质量如何，可以用热图或聚类# 用热图# add a colorbar along the heatmap with sample condition{

num_conditions

pal

cond_colors

heatmap.2(cor(raw_counts), RowSideColors=cond_colors,

trace='none', main='Sample correlations (raw)')

}#用聚类{

sampleTree = hclust(dist(cor(raw_counts)), method = "average")

pdf(file = "pre-sampleClustering.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)

dev.off()

}# 如果发现跑偏的样本，除掉它if(F){

abline(h = 20, col = "red") #画一条辅助线,h的值自定义

# 比如这里设置把高于20的切除

clust = cutreeStatic(sampleTree, cutHeight = 20, minSize = 10)

table(clust) # 0代表切除的，1代表保留的

keepSamples = (clust==1)

datExpr = datExpr0[keepSamples, ]

}

聚类+热图Low-count filtering

对样本质量满意后，我们只保留表达量丰富的基因，去除低表达量或者为0的基因# 基因在每个样本中平均表达量为1就要被过滤low_count_mask

filt_raw_counts

sprintf("Removing %d low-count genes (%d remaining).",           sum(low_count_mask), sum(!low_count_mask))#结果会返回类似：## [1] "Removing 6154 low-count genes (14802 remaining)."

Log-transforming data# 注意log使用要将真数+1 (真数就是这里的filt_raw_counts)log_counts

colnames(x) = c('gene_id', 'sample', 'value')

ggplot(x, aes(x=value, color=sample)) + geom_density()#再画个热图heatmap.2(cor(log_counts), RowSideColors=cond_colors,

trace='none', main='Sample correlations (log2-transformed)')

过滤、转换前后对比过滤并且log转换后，确实数据开始变好了

2.3 Remove non differentially-expressed genes

对于多个分组信息，需要生成几组两两组合的差异比较矩阵(取决于表型数据中的因子信息)；并且方差不显著的基因就要去除# 首先，去除方差为0的基因，因为这些基因没有表现出任何差别，还有可能对下面构建模型产生误导log_counts  0,]# 构建design矩阵为差异分析作准备mod

fit

comparisons

for (i in 1:nrow(condition_pairs)) {

comparisons[[i]]

}

# 设置一个储存差异基因的向量sig_genes

contrast_formula

contrast_mat

contrast_fit

eb

sig_genes

rownames(topTable(eb, number=Inf, p.value=0.005)))

}# 最后把差异不显著的基因去除，留下DEGslog_counts

作者：刘小泽

链接：https://www.jianshu.com/p/b05145d0020a