单细胞WGCNA

加权基因共表达网络分析 (WGCNA, Weighted correlation network analysis)是用来描述不同样品（单细胞中为cell-barcode）之间基因关联模式的系统生物学方法，可以用来鉴定高度协同变化的基因集，并根据基因集的内连性和基因集与表型之间的关联鉴定marker gene 或治疗靶点。
参考链接：https://mp.weixin.qq.com/s/z0Kr9ZbHSeRBWnSs0bJS1Q
本文仅做个人学习使用。

WGCNA相关概念

Co-expression network(共表达网络)：共表达网络定义为无向的、加权的基因网络。这样一个网络的节点对应于基因，基因之间的边代表基因表达量的相关性，加权是将相关性的绝对值提高到幂β≥1(软阈值)，加权基因共表达网络的构建以牺牲低相关性为代价，强调高相关性。
Module(模块) ：模块是高度互连的基因簇。模块对应于正相关的基因。这里的加权的网络就等于邻接矩阵。通过幂邻接转换，就强化了高相关性基因的关系，弱化了相关性基因的关系。
Connectivity(连接度)：对于每个基因，连接性(也称为度)被定义为与其他基因的连接强度之和：在共表达网络中，连接度衡量一个基因与所有其他网络基因的相关性。
Intramodular connectivity(模块内连接度)：模块内链接度衡量给定基因相对于特定模块的基因的连接或共表达程度。模块内连接度可以做为Module membership的度量。
Module eigengene E：给定模块的第一主成分，代表整个模块的基因表达谱
Module Membership(MM)：对于每个基因，我们通过将其基因表达谱与模块的Module eigengene相关性来定义Module Membership。
hub gene ：高度连接基因的缩写，根据定义，它是共表达网络模块内具有高连接度的基因。
Gene significance(GS) ：模块的显著性（module significance,Ms）：定义为模块包含的所有基因显著性性的平均值，然后比较MS，一般MS越高，说明这个模块与疾病之间的关联度越高。

官方教程

构建基因共表达网络：使用加权的表达相关性。
识别基因集：基于加权相关性，进行层级聚类分析，并根据设定标准切分聚类结果，获得不同的基因模块，用聚类树的分枝和不同颜色表示。
如果有表型信息，计算基因模块与表型的相关性，鉴定性状相关的模块。
研究模块之间的关系，从系统层面查看不同模块的互作网络。
从关键模块中选择感兴趣的驱动基因，或根据模块中已知基因的功能推测未知基因的功能。
导出TOM矩阵，绘制相关性网络图。

流程详细说明

Step 1 Loading expression data and cleaning
导入数据(datExpr0)然后进行数据处理，清除缺失值和离群值，生成可用于下一步分析的表达量数据(datExpr)。其中，goodSamplesGenes用于发现缺失值，hclust用于发现离群值。
step 2 Automatic network construction
这个过程分两步: (1)选择合适的软阈值； (2)构建表达网络。
(1)　设置一个power区间，一般为c(c(1:10), seq(from = 12, to=20, by=2))。利用pickSoftThreshold()函数对datExpr在该power区间内筛选出合适的阈值。我一般选择函数估测的阈值(powerEstimate)作为最优值。

powers = c(c(1:10), seq(from = 12, to=20, by=2))
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5,
                        networkType = "signed")

(2)　blockwiseModules()函数构建网络。需要注意的是，参数networkType 有”unsigned”、”signed”等选择，”signed”与”unsigned”方法的不同会导致阈值和网络的不同，此处作者建议选择”signed”。不过即使该步与之后TOMsimilarityFromExpr()函数均选择”signed”，最终输出用于Cytoscape可视化的网络时仍显示”undirected”。
该步结束后，datExpr中的全部基因会被划分到n个module中。其中module grey包含的是未被分配的基因，后续研究中可以不关注该module。

net = blockwiseModules(datExpr, power = sft$powerEstimate,
                       TOMType = "signed", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       networkType = "signed",
                       saveTOMFileBase = "TOM",
                       verbose = 3)
moduleColors = labels2colors(net$colors)

step 3 基因网络可视化
moduleEigengenes()绘制n-1个mudule间的关系(1表示module grey)
步骤：
(1) moduleEigengenes(datExpr, moduleColors)$eigengenes生成样本的特征基因矩阵。
(2) 对该矩阵用orderMEs()排序后，moduleEigengenes()生成图像。

TOMplot()绘制TOM图
步骤：
(1) 生成全基因不相似TOM矩阵，1-TOMsimilarityFromExpr(datExpr, powerEstimate)，将得到矩阵加幂(default)，使色彩差异更明显。
(2)TOMplot()对不相似TOM矩阵进行可视化。

step 4 特征基因矩阵与性状间的关系
该步骤计算排序后的特征基因矩阵(MET)与性状数据(Traits)间的相关性。性状数据可以是一列或多列，但必须是数字矩阵。
(1)　cor(MET,Traits)计算二者之间的相关性，结果存为moduleTraitCor。
(2)　corPvalueStudent()对moduleTraitCor添加P值。
(3)　labeledHeatmap()对相关性进行可视化。
(4)　根据相关性结果，选择与所研究性状相关性高的module为待研究mudule。

step 5 基因与性状、基因模块间的关系该步骤计算表达量矩阵(datExpr)与性状数据(Traits)、基因模块间(MET)的相关性。性状数据可以是一列或多列，但必须是数字矩阵。

part A 基因与性状间的关系
(1)　cor(datExpr,Traits)计算二者之间的相关性，结果存为geneTraitSignificance。
(2)　corPvalueStudent()对geneTraitSignificance添加P值。
part B 基因与模块间的关系
(1)　cor(datExpr,MET)计算二者之间的相关性，结果存为geneModuleMembership。
(2)　corPvalueStudent()对geneModuleMembership添加P值。
Note: 该步骤完成后，可以就step 4(4)中的待研究module，观察module的特征向量与性状数据，在基因层面上的相关性。 verboseScatterplot(abs(geneModuleMembership),abs(geneTraitSignificance))实现该步骤的可视化。相关性系数越高越好，因为与性状高度相关的基因，也是与性状相关的模型的关键基因。

step 6 hub gene
hub genes指模块中连通性（connectivity）较高的基因。高连通性的Hub基因通常为调控因子（调控网络中处于偏上游的位置）。
一般计算出KME值，并利用|kME| >=阈值（0.8）来筛选出hub gene。
(1)　connet = abs(cor(datExpr,use=”p”))^6计算基因间的相关性。intramodularConnectivity(connet, moduleColors)计算模块内连接度，即节点与同一模块内其他节点的连接度，结果记为Alldegrees1。
(2)　计算基因-性状显著性与模块内连接度之间的关系。

colorlevels=unique(moduleColors)
par(mfrow=c(2,as.integer(0.5+length(colorlevels)/2)))
par(mar = c(4,5,3,1))
for (i in c(1:length(colorlevels))) 
{
  whichmodule=colorlevels[[i]]; 
  restrict1 = (moduleColors==whichmodule);
  verboseScatterplot(Alldegrees1$kWithin[restrict1], 
                     geneTraitSignificance[restrict1,1], 
                     col=moduleColors[restrict1],
                     main=whichmodule, 
                     xlab = "Connectivity", ylab = "Gene Significance", abline = TRUE)
}

(3)　signedKME()函数计算KME，根据阈值筛选出hub gene。

datKME=signedKME(datExpr, MET, outputColumnName="MM.")
FilterGenes = abs(geneTraitSignificance[,1])> .8 & abs(datKME$MM.blue)>.8
table(FilterGenes)
hubgene <- dimnames(data.frame(datExpr))[[2]][FilterGenes]

scWGCNA

Github中创建的单细胞WGCNA分析 https://github.com/milescsmith/scWGCNA

scWGCNA输入数据

需要的是均一化之后的数据，在Seurat对象中是pbmc_small@assays$RNA@data中的数据，表型数据是pbmc_small@meta.data。
其中每个细胞被当做一个样本。但是这也是一个挑战：A）数据纬度高了计算量大；B)纬度高数据稀疏，相关性差，找不到明显的模块。
Weighted gene co-expression network analysis (WGCNA) was performed with functions in the WGCNA R package. To attenuate the effects of noise and outliers, the analyses were performed on pseudocells. calculated as averages of 4-10 cells randomly chosen within each cluster.
也就是每个cluster随机选一部分细胞构成Pseudocell（局部bulk的方法）进行分析。
单细胞转录组数据的一个特点就是纬度高，数据稀疏，除了要考虑细胞的特殊处理之外，还可以过滤基因，如只用高変基因（FindVariableFeatures()）。这里请注意，也是不推荐只选用某一群的差异基因做的，因为某一群的差异基因，已经是一个明显的模块了，这样做很可能只得到很少的模块。

基因过滤，可以挑选一部分基因做
细胞过滤，A选择某一群细胞看某一群的基因表达模块；B整个样本做；C如果细胞数据比较离散可以考虑我们上面提到的构造Pseudocell再做
数据一般使用均一化之后的

scWGCNA代码

# R版本最好3.6.3以上
setwd("D:\\6月\\单细胞\\sraData\\Seurat\\filtered_gene_bc_matrices\\hg19")
library(WGCNA)
library(Seurat)
library(tidyverse)
library(reshape2)
library(stringr)
library(patchwork)
library(dplyr)

# load the data
pbmc.data <- read10x('./')
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
head(pbmc@meta.data)

# Seurat常规分析流程
## QC
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)  # 绘图
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

## Normalizing
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

## Identification of highly variable features (feature selection)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(pbmc), 10)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

## Scaling the data
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

## linear dimensional reduction线性降维
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
DimPlot(pbmc, reduction = "pca")
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)

## 挑选维度
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:15)
ElbowPlot(pbmc)

## 聚类
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
head(Idents(pbmc), 5)

## non-linear dimensional reduction非线性降维
pbmc <- RunUMAP(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap")
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")

## 差异基因分析
# find all markers of cluster 1
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)
# find markers for every cluster compared to all remaining cells, report only the positive ones
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
cluster1.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
# you can plot raw counts as well
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()

## 类群注释
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono", 
    "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
saveRDS(pbmc, file = "../output/pbmc3k_final.rds")

# 加载rds文件
pbmc <- readRDS(../output/pbmc_tutorial.rds)
pbmc
head(pbmc@meta.data)

# 细胞数较多，且单细胞表达矩阵稀疏，先讲细胞筛选一下
# 采用pseudocell概念
datadf <- as.matrix(pbmc@assays$RNA@data)
idd1 <- pbmc@meta.data
Inter.id1<-cbind(rownames(idd1),idd1$seurat_clusters)
rownames(Inter.id1)<-rownames(idd1)
colnames(Inter.id1)<-c("CellID","Celltype")
Inter.id1<-as.data.frame(Inter.id1)
head(Inter.id1)
Inter1<-datadf[,Inter.id1$CellID]
Inter2<-as.matrix(Inter1)
Inter2[1:4,1:4]

pseudocell.size = 10
new_ids_list1 = list()
length(levels(Inter.id1$Celltype))
for (i in 1:length(levels(Inter.id1$Celltype))) {
  cluster_id = levels(Inter.id1$Celltype)[i]
  cluster_cells <- rownames(Inter.id1[Inter.id1$Celltype == cluster_id,])
  cluster_size <- length(cluster_cells)     
  pseudo_ids <- floor(seq_along(cluster_cells)/pseudocell.size)
  pseudo_ids <- paste0(cluster_id, "_Cell", pseudo_ids)
  names(pseudo_ids) <- sample(cluster_cells)    
  new_ids_list1[[i]] <- pseudo_ids      
}

new_ids <- unlist(new_ids_list1)
new_ids <- as.data.frame(new_ids)
head(new_ids)
new_ids_length <- table(new_ids)
new_ids_length
new_colnames <- rownames(new_ids)

# 筛选细胞
gc()
colnames(datadf)  
all.data<-datadf[,as.character(new_colnames)] ###add
all.data <- t(all.data)###add
new.data<-aggregate(list(all.data[,1:length(all.data[1,])]),
                    list(name=new_ids[,1]),FUN=mean)
rownames(new.data)<-new.data$name
new.data<-new.data[,-1]
new_ids_length<-as.matrix(new_ids_length)##
short<-which(new_ids_length<10)##
new_good_ids<-as.matrix(new_ids_length[-short,])##
result<-t(new.data)[,rownames(new_good_ids)]
dim(result)

# 基因过滤
pbmc <- FindVariableFeatures(pbmc,nfeatures = 5000)
colnames(result)[grepl("[12]_Cel",colnames(result))]
Cluster1 <- result[intersect(Seurat::VariableFeatures(pbmc),rownames(result)),]

# WGCNA参数设置
type = "unsigned"  # 官方推荐 "signed" 或 "signed hybrid"
corType = "pearson" # 相关性计算  官方推荐 biweight mid-correlation & bicor  corType: pearson or bicor 
corFnc = ifelse(corType=="pearson", cor, bicor)
corFnc
maxPOutliers = ifelse(corType=="pearson",1,0.05) # 对二元变量，如样本性状信息计算相关性时， # 或基因表达严重依赖于疾病状态时，需设置下面参数
# 关联样品性状的二元变量时，设置
robustY = ifelse(corType=="pearson",T,F)
dataExpr  <- as.matrix(Cluster1)

# 重新再进行一次筛选
## 筛选中位绝对偏差前75%的基因，至少MAD大于0.01
## 筛选后会降低运算量，也会失去部分信息
## 也可不做筛选，使MAD大于0即可
m.mad <- apply(dataExpr,1,mad)
dataExprVar <- dataExpr[which(m.mad > 
                                max(quantile(m.mad, probs=seq(0, 1, 0.25))[2],0.01)),]
## 转换为样品在行，基因在列的矩阵
dataExpr <- as.data.frame(t(dataExprVar))
dim(dataExpr)
head(dataExpr)[,1:8]
## 检测缺失值
gsg = goodSamplesGenes(dataExpr, verbose = 3)
gsg$allOK
gsg$goodSamples

if (!gsg$allOK){
  # Optionally, print the gene and sample names that were removed:
  if (sum(!gsg$goodGenes)>0) 
    printFlush(paste("Removing genes:", 
                     paste(names(dataExpr)[!gsg$goodGenes], collapse = ",")));
  if (sum(!gsg$goodSamples)>0) 
    printFlush(paste("Removing samples:", 
                     paste(rownames(dataExpr)[!gsg$goodSamples], collapse = ",")));
  # Remove the offending genes and samples from the data:
 dataExpr = dataExpr[gsg$goodSamples, gsg$goodGenes]
}

nGenes = ncol(dataExpr)
nSamples = nrow(dataExpr)
dim(dataExpr)
# 在筛选的时候，有些基因是要保留的，哪些基因呢？就是那些你关注的基因

## 查看是否有离群样品
sampleTree = hclust(dist(dataExpr), method = "average")
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="")

powers = c(c(1:10), seq(from = 12, to=30, by=2))
# 需要重点理解pickSoftThreshold函数及其返回的对象，最佳的beta值就是sft$powerEstimate
sft = pickSoftThreshold(dataExpr, powerVector=powers, 
                        networkType="signed", verbose=5)
par(mfrow = c(1,2))
cex1 = 0.9
# 横轴是Soft threshold (power)，纵轴是无标度网络的评估参数，数值越高，
# 网络越符合无标度特征 (non-scale)
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")
# 筛选标准。R-square=0.85
abline(h=0.85,col="red")

# Soft threshold与平均连通性
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")
power = sft$powerEstimate
softPower  = power
softPower

参数beta取值默认是1到30，上述图形的横轴均代表权重参数β，左图纵轴代表对应的网络中log(k)与log(p(k))相关系数的平方。相关系数的平方越高，说明该网络越逼近无网路尺度的分布。右图的纵轴代表对应的基因模块中所有基因邻接函数的均值。最佳的beta值就是sftpowerEstimate¨G10G¨G11G参数beta取值默认是1到30，上述图形的横轴均代表权重参数β，左图纵轴代表对应的网络中log(k)与log(p(k))相关系数的平方。相关系数的平方越高，说明该网络越逼近无网路尺度的分布。右图的纵轴代表对应的基因模块中所有基因邻接函数的均值。最佳的beta值就是sftpowerEstimate，已经被保存到变量了，不需要知道具体是什么，后面的代码都用这个即可，在本例子里面是5。

# 无向网络在power小于15或有向网络power小于30内，没有一个power值可以使
# 无标度网络图谱结构R^2达到0.8，平均连接度较高如在100以上，可能是由于
# 部分样品与其他样品差别太大。这可能由批次效应、样品异质性或实验条件对
# 表达影响太大等造成。可以通过绘制样品聚类查看分组信息和有无异常样品。
# 如果这确实是由有意义的生物变化引起的，也可以使用下面的经验power值。
if (is.na(power)){
  power = ifelse(nSamples<20, ifelse(type == "unsigned", 9, 18),
                 ifelse(nSamples<30, ifelse(type == "unsigned", 8, 16),
                        ifelse(nSamples<40, ifelse(type == "unsigned", 7, 14),
                               ifelse(type == "unsigned", 6, 12))       
                 )
  )
}
power

#一步法网络构建：One-step network construction and module detection##
# power: 上一步计算的软阈值
# maxBlockSize: 计算机能处理的最大模块的基因数量 (默认5000)；
#  4G内存电脑可处理8000-10000个，16G内存电脑可以处理2万个，32G内存电脑可
#  以处理3万个
#  计算资源允许的情况下最好放在一个block里面。
# corType: pearson or bicor
# numericLabels: 返回数字而不是颜色作为模块的名字，后面可以再转换为颜色
# saveTOMs：最耗费时间的计算，存储起来，供后续使用
# mergeCutHeight: 合并模块的阈值，越大模块越少
# type = unsigned
cor <- WGCNA::cor
net = blockwiseModules(dataExpr, power = power, maxBlockSize = nGenes,#nGenes
                       TOMType = "unsigned", minModuleSize = 10,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs=TRUE, corType = corType, 
                       maxPOutliers=maxPOutliers, loadTOMs=TRUE,
                       saveTOMFileBase = paste0("dataExpr", ".tom"),
                       verbose = 3)

# 对象探索
table(net$colors)
net$unmergedColors
head(net$MEs)
net$goodSamples
net$goodGenes
net$TOMFiles
net$blockGenes
net$blocks
net$MEsOK
head(net$MEs)
table(net$colors)

## 灰色的为**未分类**到模块的基因。灰色太多一定程度显示基因选择的问题
# Convert labels to colors for plotting
moduleLabels = net$colors
moduleColors = labels2colors(moduleLabels)
moduleColors
# Plot the dendrogram and the module colors underneath
# 如果对结果不满意，还可以recutBlockwiseTrees，节省计算时间
plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

# module eigengene, 可以绘制线图，作为每个模块的基因表达趋势的展示
MEs = net$MEs
### 不需要重新计算，改下列名字就好
### 官方教程是重新计算的，起始可以不用这么麻烦
MEs_col = MEs
colnames(MEs_col) = paste0("ME", labels2colors(
  as.numeric(str_replace_all(colnames(MEs),"ME",""))))
MEs_col = orderMEs(MEs_col)
# 根据基因间表达量进行聚类所得到的各模块间的相关性图
# marDendro/marHeatmap 设置下、左、上、右的边距
head(MEs_col)
?plotEigengeneNetworks
plotEigengeneNetworks(MEs, "Eigengene adjacency heatmap", 
                      marDendro = c(3,3,2,4),
                      marHeatmap = c(3,4,2,2),
                      plotDendrograms = T,
                      xLabelsAngle = 90)

# 计算每个模块的特征向量基因，为某一特定模块第一主成分基因E。代表了该模块内基因表达的整体水平
MEList = moduleEigengenes(dataExpr, colors = dynamicColors)
MEs = MEList$eigengenes
# 计算根据模块特征向量基因计算模块相异度：
MEDiss = 1-cor(MEs);
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average");
# Plot the result

plotEigengeneNetworks(MEs, 
                      "Eigengene adjacency heatmap", 
                      marHeatmap = c(3,4,2,2), 
                      plotDendrograms = FALSE, 
                      xLabelsAngle = 90) 

# 画出指定模块表达量的热图
which.module="turquoise"; 
ME=mergedMEs[, paste("ME",which.module, sep="")]
par(mfrow=c(2,1), mar=c(0,4.1,4,2.05))
plotMat(t(scale(dataExpr[,moduleColors==which.module ]) ),
        nrgcols=30,rlabels=F,rcols=which.module,
        main=which.module, cex.main=2)
par(mar=c(2,2.3,0.5,0.8))
barplot(ME, col=which.module, main="", cex.main=2,
        ylab="eigengene expression",xlab="array sample")

# 所有基因模块关系
load(net$TOMFiles, verbose=T)
## Loading objects:
##   TOM
TOM <- as.matrix(TOM)
TOM[1:4,1:4]
#dim(TOM2)
dissTOM = 1-TOM
# Transform dissTOM with a power to make moderately strong 
# connections more visible in the heatmap
plotTOM = dissTOM^7
# Set diagonal to NA for a nicer plot
diag(plotTOM) = NA
# Call the plot function
table(moduleColors)
# 这一部分特别耗时，行列同时做层级聚类
TOMplot(plotTOM, net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]], 
        main = "Network heatmap plot, all genes")

# 由于之前更改了细胞名称，所以需要构造表型数据
Cluster1[1:4,1:4]
mypbmc<- CreateSeuratObject(Cluster1)
mypbmc
mypbmc[["percent.mt"]] <- PercentageFeatureSet(mypbmc, pattern = "^CD")
mypbmc <- FindVariableFeatures(mypbmc, selection.method = "vst", nfeatures = 2000)
mypbmc %>% NormalizeData( normalization.method = "LogNormalize", scale.factor = 10000)%>%
  FindVariableFeatures( selection.method = "vst", nfeatures = 2000) %>%
  ScaleData(features=VariableFeatures(mypbmc),vars.to.regress = "percent.mt")  %>%
  RunPCA(features = VariableFeatures(object = mypbmc)) %>%
  FindNeighbors( dims = 1:10) %>%
  FindClusters( resolution = 0.5) %>%
  BuildClusterTree() %>%
  RunUMAP( dims = 1:10)  -> mypbmc
head(mypbmc@meta.data)
head(mypbmc@meta.data)

# 通过模块与各种表型的相关系数，可以很清楚的挑选自己感兴趣的模块进行下游分析了。
moduleTraitCor_noFP <- cor(mergedMEs, mypbmc@meta.data, use = "p");
moduleTraitPvalue_noFP = corPvalueStudent(moduleTraitCor_noFP, nSamples); 
textMatrix_noFP <- paste(signif(moduleTraitCor_noFP, 2), "\n(", signif(moduleTraitPvalue_noFP, 1), ")", sep = ""); 
par(mar = c(10, 8.5, 3, 3)); 
labeledHeatmap(Matrix = moduleTraitCor_noFP, 
               xLabels = names(mypbmc@meta.data), 
               yLabels = names(mergedMEs), 
               ySymbols = names(mergedMEs), 
               colorLabels = FALSE, 
               colors = blueWhiteRed(50), 
               textMatrix = textMatrix_noFP,
               setStdMargins = FALSE, 
               cex.text = 0.65, 
               zlim = c(-1,1), 
               main = paste("Module-trait relationships"))

# 根据性状与模块特征向量基因的相关性及pvalue来挖掘与性状相关的模块
library(pheatmap)
cor_ADR <- signif(WGCNA::cor(mypbmc@meta.data,mergedMEs,use="p",method="pearson"),5)
p.values <- corPvalueStudent(cor_ADR,nSamples=nrow(mypbmc@meta.data))
pheatmap(cor_ADR,display_numbers = matrix(ifelse(p.values <= 0.01, "**", ifelse(p.values<= 0.05 ,"*"," ")), nrow(p.values)),fontsize=18)

# 根据基因网络显著性，也就是性状与每个基因表达量相关性在各个模块的均值作为该性状在该模块的显著性，显著性最大的那个模块与该性状最相关
GS1 <- as.numeric(WGCNA::cor(mypbmc@meta.data[,3],dataExpr,use="p",method="pearson"))
# 显著性是绝对值：
GeneSignificance <- abs(GS1)
length(GeneSignificance)
length(mergedColors)
mypbmc
dim(dataExpr)
dim(Cluster1)
dim(mypbmc@meta.data)
# 获得该性状在每个模块中的显著性：
ModuleSignificance <- tapply(GeneSignificance,mergedColors,mean,na.rm=T)
ModuleSignificance

# 寻找与该性状相关的枢纽基因（hub genes）,首先计算基因的内部连接度和模块身份，内部连接度衡量的是基因在模块内部的地位，而模块身份表明基因属于哪个模块。
# 计算每个基因模块内部连接度，也就是基因直接两两加权相关性。
ADJ1=abs(cor(dataExpr,use="p"))^softPower 
# 根据上面结果和基因所属模块信息获得连接度：
# 整体连接度 kTotal，模块内部连接度：kWithin，kOut=kTotal-kWithin， kDiff=kIn-kOut=2*kIN-kTotal
Alldegrees1=intramodularConnectivity(ADJ1, moduleColors) 
head(Alldegrees1)

# 注意模块内基于特征向量基因连接度评估模块内其他基因：de ne a module eigengene-based connectivity measure for each gene as the correlation between a the gene expression and the module eigengene
# 如 brown 模块内：kM Ebrown(i) = cor(xi, MEbrown) ， xi is the gene expression pro le of gene i and M Ebrown is the module eigengene of the brown module
# 而 module membership 与内部连接度不同。MM 衡量了基因在全局网络中的位置。
datKME=signedKME(dataExpr, MEs, outputColumnName="MM.")
datKME[1:4,1:4]

# 注意模块内基于特征向量基因连接度评估模块内其他基因：de ne a module eigengene-based connectivity measure for each gene as the correlation between a the gene expression and the module eigengene
# 如 brown 模块内：kM Ebrown(i) = cor(xi, MEbrown) ， xi is the gene expression pro le of gene i and M Ebrown is the module eigengene of the brown module
# 而 module membership 与内部连接度不同。MM 衡量了基因在全局网络中的位置。
datKME=signedKME(dataExpr, MEs, outputColumnName="MM.")
datKME[1:4,1:4]

# 选择特定模块的基因
table(moduleColors)
module = "yellow";
# Select module probes
probes = colnames(dataExpr) ## 我们例子里面的probe就是基因名
inModule = (moduleColors==module);
modProbes = probes[inModule]; 
modProbes

# 导出作图文件，主要模块里面的基因直接的相互作用关系信息可以导出到cytoscape,VisANT等网络可视化软件。
## 也是提取指定模块的基因名
# Select the corresponding Topological Overlap
modTOM = TOM[inModule, inModule];
dimnames(modTOM) = list(modProbes, modProbes)
vis = exportNetworkToVisANT(modTOM,
                            file = paste("VisANTInput-", module, ".txt", sep=""),
                            weighted = TRUE,
                            threshold = 0)
?exportNetworkToCytoscape
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.005,
  nodeNames = modProbes, 
  nodeAttr = moduleColors[inModule]
);

WGCNA四步骤详解

确定软阈值，构建邻接矩阵

https://mp.weixin.qq.com/s/rOr8nVEYBY3P4Z9CLii67g
无尺度分布，随机网络是可以用一个值来衡量大多数节点之间的距离，或者是6或者是60，也就是可以用一个尺度去衡量他，可以称为尺度网络。而有hub的网络，没办法来衡量两个节点之间的距离，叫作无尺度分布，该分布又叫做幂律分布，也就是二八定律，长尾定律都是幂律分布的口头化呈现。
而蛋白和蛋白之间的互作也是符合无尺度分布的，WGCNA就是让基因间的联系符合无尺度分布。

手动计算软阈值，挑选合适的软阈值是为了让基因间的相关性符合无尺度网络分布。

# 手动计算软阈值过程
mypick <- function(powerVector,datExpr){
  power <- powerVector
  cor<-stats::cor
  ADJ=abs(cor(datExpr,use = 'p'))^power
  k=apply(ADJ,2,sum) -1
  cut1=cut(k,10)
  binned.k=tapply(k,cut1,mean)
  freq1=tapply(k,cut1,length)/length(k)
  xx= as.vector(log10(binned.k))
  lm1=lm(as.numeric(log10(freq1+.000000001))~ xx )
  return(data.frame(Power=power,
                    SFT.R.sq=as.character(round(summary(lm1)$adj.r.squared,2)),
                    slope=round(lm1$coefficients[[2]],2),
                    mean.k=mean(k)))
}
mypick(10,datExpr)
powers = c(c(1:10), seq(from = 12, to=20, by=2))
do.call(rbind,lapply(powers,mypick,datExpr))

sizeGrWindow(9, 5)
par(mfrow = c(1,2))
cex1 = 0.85
# 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="red")
# 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")

# 官方脚本计算软阈值
powers = c(c(1:10), seq(from = 12, to=20, by=2))
sft = pickSoftThreshold(datExpr, powerVector = powers)
# 返回结果为一个列表
# 第一个：确定的软阈值
sft$powerEstimate
# 第二个：计算得到的结果表格
sft$fitIndices

基因相关性中存在unsigned和signed的，指的是相关性正负的问题。作者推荐signed，因为更符合真实情况，但是不绝对，如果是单细胞数据，用unsigned比较合适。

把邻接矩阵变成拓扑重叠矩阵

## 加载R包
library(WGCNA)
## 确定软阈值
softPower = 6
## 得到邻接矩阵(就相当于相关性分析)
adjacency = adjacency(datExpr, power = softPower)
## 转换为拓扑矩阵
TOM = TOMsimilarity(adjacency)
dissTOM = 1-TOM
geneTree = hclust(as.dist(dissTOM), method = "average")
plot(geneTree,labels = FALSE, hang = 0.04)

仅仅通过表达的相关性矩阵，不足于反应体内真实的情况。
假如基因1和基因2的表达相关性是0.8，基因2和基因3的表达相关性也是0.8,是不能说这两个关系是相等的，因为基因1和基因2之间还有很多共同的小伙伴跟着两个基因建立联系。也就是说，我们评价两个基因相关性的时候，不能仅仅看这个两个基因的表达相关性，还要考虑跟其他基因之间的互作。
WGCNA相关性 TOM矩阵计算

根据拓扑重叠矩阵识别基因模块

WGCNA中识别模块是使用层次聚类。

geneTree = hclust(as.dist(dissTOM), method = "average")
# 其中methods有很多种，比如average，single，complete
# 在计算两个cluster距离的时候，single选择的离得最近的点，complete选择的是最远的点
# 然后average是计算每个点然后取平均值。
# 而WGCNA中这里选择的是average。
sizeGrWindow(12,9)
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",labels = FALSE, hang = 0.04);

# 动态切合算法
minModuleSize = 30
dynamicMods = cutreeDynamic(dendro = geneTree, 
                            distM = dissTOM,
                            deepSplit = 2, 
                            pamStage = TRUE,
                            pamRespectsDendro = FALSE,
                            minClusterSize = minModuleSize)
# minModuleSize规定最小的模块不小于30
# dendro = geneTree,  distM = dissTOM 输入的是前面的树和TOM矩阵
# deepSplit 是切割强度，越大得到的模块越多
# pamStage = TRUE, pamRespectsDendro = FALSE,
这是对PAM方法的限制，pamStage = TRUE是使用PAM方法，pamRespectsDendro = FALSE，是说PAM不会考虑树的感受

table(dynamicMods)
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)

# 合并模块（模块太多了，还可以把相近的合并）
merge = mergeCloseModules(datExpr, 
                          dynamicColors, 
                          cutHeight = 0.25)

eigengene意义

eigengene指每个模块基因和样本矩阵进行PCA分析后的第一主成分。

moduleEigengenes(datExpr, moduleColors)
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes # 行是样本，列是模块的eigengene

# 具体计算过程
table(moduleColors)
datablack <- datExpr[,moduleColors == "black"]
### 缺失值填充
### 要求就是行是基因，列是样本
datModule = t(as.matrix(datablack))
datModule = impute::impute.knn(datModule, k = min(10, nrow(datModule)-1))
## 提取数据，行是样本，列是基因
datModule_knn = t(datModule$data)
pc <- prcomp(datModule_knn)
pc$x[, 1] # 提取第一主成分

# WGCNA在进行eigengene分析时，采用的是奇异值分解
datModule_scale= scale(datModule_knn)
svd = svd(datModule_scale)
svd$u[,1]

性状关联。模块里面的基因，经过压缩，变成了一个eigengene，这样跟性状信息可以求相关性。有了相关性，我们就可以得到性状和模块的相关性热图，这也是WGCNA里面最重要的一张图。
hubgene，eigengene的存在，让找每个模块中的是hub gene变得简单。hubgene指那个模块中跟eigengene 相关性最强的基因。
迭代WGCNA。eigengene的出现，我们筛选模块里面高价值的基因就有了标准，我们可以把模块中的每个基因都跟eigengene计算相关性，把小于0.8的全部剔除(人为定义) 每个模块都这么做，现在就到了一个新的行是基因，列是样本的矩阵。重新来做一次WGCNA，再来筛选一次，直到无法筛选基因为止。此时，就得到了都比较纯净模块。

单细胞转录组WGCNA应该如何做？