hdWGCNA与WGNCA分析的关键点很类似，主要有以下几点：
- 构建无标度网络，识别关键模块基因。这里的无标度网络就像是社会中的人与人之间网络关系，大多数人之间的关系是普通的，但有少部分的人具有很强的"引力"，与他们链接在一起的有一大帮子人。这里的少部分人在WGCNA中就相当于是关键模块基因，是非常重要的节点，如果没有了这个节点这个网络就可能会“瘫痪”。
- 在构建这个无标度网络的时候采用了加权共表达的方式，并且由定义的软阈值去计算加权网络。加权共表达的方式非常好理解，比如人与人之间的交流肯定存在一定的强弱关系，不会只是好与不好的正反两种情况。
- 软阈值的作用是为了更好的放大或者缩小不同节点之间的相关性情况，从而减少在未使用软阈值情况下节点因为稍未达到阈值而被认为“不重要”的情况发生。
hdWGCNA分析流程
1.导入
rm(list = ls())
library(WGCNA)
library(hdWGCNA)
library(tidyverse)
library(cowplot)
library(patchwork)
library(qs)
library(Seurat)
set.seed(12345)
#enableWGCNAThreads(nThreads = 8)
#这里加载的是seurat对象，替换自己的数据即可
scRNA <- qread("./9-CD4+T/CD4+t_final.qs")
#检查一下自己导入进来的数据
Idents(scRNA) <- "celltype"
DimPlot(scRNA,reduction = 'umap',
label = TRUE,pt.size = 0.5) +NoLegend()
2.为WGCNA设置Seurat对象
- WGCNA分析的时候会把信息储存在seurat对象的@misc槽中；
- variable: 使用存储在Seurat对象的VariableFeatures中的基因；
- fraction: 使用在整个数据集或每组细胞中表达的基因，由 group.by 指定;
- custom: 使用在Custom 列表中指定的基因；
- 一个seurat对象可以包含多个hdWGCNA实验对象；
# V5版本需要这行代码，V4不需要
scRNA <- SeuratObject::UpdateSeuratObject(scRNA)
seurat_obj <- SetupForWGCNA(
scRNA,
gene_select = "fraction", # 默认fraction;其他有variable(Seurat对象中存储的基因),custom(自定义)
fraction = 0.05, # fraction of cells that a gene needs to be expressed in order to be included
wgcna_name = "CD4+T" # the name of the hdWGCNA experiment
)
3.构建metacells
metacells是由来自同一个生物样本的、相似细胞组成的小群体聚合而成的。该过程使用k最近邻(k-Nearest Neighbors, KNN)算法来识别相似细胞的群体，然后计算这些细胞的平均或总表达量，从而生成一个metacell基因表达矩阵。与原始表达矩阵相比，metacell表达矩阵的稀疏性大大降低，更适合用于后续分析。
# 各组构建metacell
seurat_obj <- MetacellsByGroups(
seurat_obj = seurat_obj,
group.by = c("celltype", "orig.ident"), #指定seurat_obj@meta.data中要分组的列
reduction = 'harmony', # 选择要执行KNN的降维
k = 25, # 最近邻居参数
max_shared = 10, # 两个metacell之间共享细胞的最大数目
ident.group = 'celltype', # 设置metacell安全对象的标识
min_cells = 100 # 排除数量小于100的细胞亚群
)
# normalize metacell expression matrix:
seurat_obj <- NormalizeMetacells(seurat_obj)
4.共表达网络分析
设置表达式矩阵，使用hdWGCNA对目标细胞亚群进行共表达网络分析
seurat_obj <- SetDatExpr(
seurat_obj,
group_name = "Treg", # the name of the group of interest in the group.by column
group.by='celltype', # the metadata column containing the cell type info. This same column should have also been used in MetacellsByGroups
assay = 'RNA', # using RNA assay
slot = 'data' # using normalized data
)
# # 选择多个group_names
# seurat_obj <- SetDatExpr(
# seurat_obj,
# group_name = c("Treg", "Tm"),
# group.by='celltype'
# )
4.1 选择软阈值
# Test different soft powers:
seurat_obj <- TestSoftPowers(
seurat_obj,
powers = c(seq(1, 10, by = 1), seq(12, 30, by = 2)),
networkType = 'unsigned'# you can also use "unsigned" or "signed hybrid"
)
# plot the results:
plot_list <- PlotSoftPowers(seurat_obj)
# 1 1 0.4193109 -2.829398 0.9708088 551.0532441 532.1816874 1079.59624
# 2 2 0.7722457 -2.846229 0.9935455 90.8116702 80.0299812 331.23789
# 3 3 0.8469701 -2.813601 0.9812133 19.9509707 15.0395890 127.98068
# 4 4 0.8914758 -2.532259 0.9884110 5.4419966 3.2775043 57.54649
# 5 5 0.9024030 -2.199419 0.9834531 1.7820619 0.7925626 28.86319
# 6 6 0.9511793 -1.813044 0.9753076 0.6868952 0.2083341 15.72666
# assemble with patchwork
wrap_plots(plot_list, ncol=2)
# check以下数据
power_table <- GetPowerTable(seurat_obj)
head(power_table)
# Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
# 1 1 0.02536182 3.273051 0.9541434 4370.9149 4379.0629 4736.8927
# 2 2 0.11091306 -3.571441 0.8008960 2322.5480 2286.2454 2871.2953
# 3 3 0.50454728 -4.960822 0.8035027 1286.6453 1241.8414 1898.5501
# 4 4 0.79569568 -4.812735 0.9183803 740.0525 697.1193 1338.0185
# 5 5 0.86641323 -4.110731 0.9517671 440.6141 402.5530 985.0984
# 6 6 0.88593187 -3.582879 0.9624951 270.9020 237.8831 750.2825
WGCNA和hdWGCNA的一般原则是选择使尺度自由拓扑模型拟合度(Scale Free Topology Model Fit)大于或等于0.8的最低软阈值(soft power threshold)。在构建网络时，如果用户没有提供软阈值，程序会自动选择一个软阈值。
4.2 构建共表达网络
# 如果没有指定软阈值，construcNetwork会自动指定软阈值
# construct co-expression network:
seurat_obj <- ConstructNetwork(
seurat_obj,
soft_power = 4, # 自定义了4，如果是自动选择的话可能会是3
tom_outdir = "TOM",
tom_name = 'Treg', # name of the topoligical overlap matrix written to disk
overwrite_tom = TRUE# 允许覆盖已存在的同名文件
)
# 可视化WGCNA树状图
# “灰色”模块由那些未被归入任何共表达模块的基因组成。对于所有下游分析和解释，应忽略灰色模块。
PlotDendrogram(seurat_obj, main='Treg hdWGCNA Dendrogram')
# 可选：检查拓扑重叠矩阵(topoligcal overlap matrix，TOM)
# TOM <- GetTOM(seurat_obj)
# TOM
4.3 计算模块特征基因
模块特征基因(Module Eigengenes，MEs)是用于总结整个共表达模块基因表达谱的常用指标。模块特征基因是通过在每个模块的基因表达矩阵子集上执行主成分(PCA)分析来计算的。这些PCA矩阵的第一个主成分就是MEs。此外对MEs应用Harmony批量校正，从而得到harmony后的模块特征基因(hMEs）。
# 需要先运行ScaleData，否则harmony会报错:
# seurat_obj <- ScaleData(seurat_obj, features=VariableFeatures(seurat_obj))
# 计算完整单细胞数据集中的所有MEs
seurat_obj <- ModuleEigengenes(
seurat_obj,
group.by.vars="orig.ident"
)
# 协调模特征基因:
# 允许用户对MEs应用Harmony批量校正，生成协调模块特征基因(hMEs)
hMEs <- GetMEs(seurat_obj)
# module eigengenes:
#MEs <- GetMEs(seurat_obj, harmonized=FALSE)
4.4计算模块连接性
在共表达网络分析中，通常希望关注“枢纽基因”，即在每个模块内高度连接的基因。因此，希望确定每个基因的基于特征基因(eigengene)的连接性，也称为kME。hdWGCNA提供了ModuleConnectivity 函数，用于在完整的单细胞数据集(而不是metacell数据集)中计算基因的kME值。这个函数本质上是计算基因与模块特征基因之间的成对相关性。虽然可以在整个数据集中计算所有细胞的kME，但建议在之前用于运行ConstructNetwork的细胞类型或分组中计算kME，这里选择了Treg细胞。
# 计算基于特征基因的连接性(kME)：
# 关注枢纽基因
seurat_obj <- ModuleConnectivity(
seurat_obj,
group.by = 'celltype',
group_name = 'Treg'
)
# 模块重命名
seurat_obj <- ResetModuleNames(
seurat_obj,
new_name = "Treg_NEW"
)
# 绘制每个模块按kME排序的基因
p <- PlotKMEs(seurat_obj, ncol=4)
p
4.5获取模块内部信息
# 获取模块内部信息:
# 个人认为这一部分很关键，毕竟咱们就是想要得到不同模块的基因
# 这一步去除了不需要的灰色模块基因
modules <- GetModules(seurat_obj) %>%
subset(module != 'grey')
# 显示前六列:
head(modules[,1:6])
# gene_name module color kME_grey kME_Treg_NEW1 kME_Treg_NEW2
# ISG15 ISG15 Treg_NEW1 turquoise 0.09485063 0.31618006 0.2177907
# TNFRSF18 TNFRSF18 Treg_NEW1 turquoise 0.12119087 0.39886246 0.4605542
# TNFRSF4 TNFRSF4 Treg_NEW1 turquoise 0.08844463 0.35922337 0.3728684
# SDF4 SDF4 Treg_NEW2 black 0.11518097 0.11212155 0.1883993
# B3GALT6 B3GALT6 Treg_NEW3 purple 0.03314139 0.08610811 0.1067775
# AURKAIP1 AURKAIP1 Treg_NEW1 turquoise 0.09062613 0.26244827 0.1252306
# 得到枢纽基因
# 可以提取按kME排序的前N个枢纽基因的表格,这里选择了10
hub_df <- GetHubGenes(seurat_obj, n_hubs = 10)
head(hub_df)
# gene_name module kME
# 1 GAPDH Treg_NEW1 0.6160237
# 2 S100A4 Treg_NEW1 0.5886924
# 3 MYL6 Treg_NEW1 0.5558792
# 4 TMSB10 Treg_NEW1 0.5371290
# 5 IL32 Treg_NEW1 0.5161320
# 6 ARPC1B Treg_NEW1 0.5138853
# 保存数据
qsave(seurat_obj, 'hdWGCNA_object.qs')
4.6计算hub基因siganture得分
# 计算每个模块前25个枢纽基因的kME得分
# 使用UCell方法
library(UCell)
seurat_obj <- ModuleExprScore(
seurat_obj,
n_genes = 25,
method='UCell' # Seurat方法(AddModuleScore)
)
5.可视化
模块特征图
# 每个模块制作hMEs的特征图
plot_list <- ModuleFeaturePlot(
seurat_obj,
features='hMEs', # plot the hMEs
order=TRUE # order so the points with highest hMEs are on top
)
# stitch together with patchwork
wrap_plots(plot_list, ncol=4)
特征图是常用的可视化方法，可以直接在降维空间中展示感兴趣的特性。hdWGCNA包含一个ModuleFeaturePlot函数，用于为每个共表达模块构建特征图，并以每个模块唯一分配的颜色进行着色。
相同函数绘制hub基因特征得分
# 每个模块制作hub scores的特征图
plot_list <- ModuleFeaturePlot(
seurat_obj,
features='scores', # plot the hub gene scores
order='shuffle', # order so cells are shuffled
ucell = TRUE # depending on Seurat vs UCell for gene scoring
)
# stitch together with patchwork
wrap_plots(plot_list, ncol=4)
跟上面的结果相比hub基因更加均匀分布在细胞中
每个模块在不同细胞亚群中的情况
# 每个模块在不同样本中的情况
seurat_obj$cluster <- do.call(rbind, strsplit(as.character(seurat_obj$orig.ident), ' '))[,1]
ModuleRadarPlot(
seurat_obj,
group.by = 'cluster',
barcodes = seurat_obj@meta.data %>%
subset(celltype == 'Treg') %>%
rownames(),
axis.label.size=4,
grid.label.size=4
)
相关性分析
# 查看模块相关图
ModuleCorrelogram(seurat_obj)
气泡图
# get hMEs from seurat object
MEs <- GetMEs(seurat_obj, harmonized=TRUE)
modules <- GetModules(seurat_obj)
mods <- levels(modules$module); mods <- mods[mods != 'grey']
# add hMEs to Seurat meta-data:
seurat_obj@meta.data <- cbind(seurat_obj@meta.data, MEs)
# plot with Seurat's DotPlot function
p <- DotPlot(seurat_obj, features=mods, group.by = 'celltype')
# flip the x/y axes, rotate the axis labels, and change color scheme:
p <- p +
RotatedAxis() +
scale_color_gradient2(high='red', mid='grey95', low='blue')
# plot output
p
单模块的网络图
# 使用ModuleNetworkPlot可视化每个模块前50(数值可自定)的hub gene
ModuleNetworkPlot(
seurat_obj,
outdir='ModuleNetworks', # new folder name
n_inner = 20, # number of genes in inner ring
n_outer = 30, # number of genes in outer ring
n_conns = Inf, # show all of the connections
plot_size=c(10,10), # larger plotting area
vertex.label.cex=1 # font size
)
会导出一个文件，文件中有每个模块前50的枢纽基因
结合hub基因的网络图
options(future.globals.maxSize = 5 * 1024^3) # 5GB
# hubgene network(基因数可自定)
HubGeneNetworkPlot(
seurat_obj,
n_hubs = 2,
n_other=2,
edge_prop = 0.75,
mods = 'all'
)
# 可以选择模块数
g <- HubGeneNetworkPlot(seurat_obj, return_graph=TRUE)
# get the list of modules:
modules <- GetModules(seurat_obj)
mods <- levels(modules$module); mods <- mods[mods != 'grey']
# hubgene network
HubGeneNetworkPlot(
seurat_obj,
n_hubs = 2,
n_other= 2,
edge_prop = 0.75,
mods = mods[1:5] # only select 5 modules
)
使用HubGeneNetworkPlot函数将所有模块合并在一起制作网络图，此函数可以制定展示特定数量的中心基因和其他基因。
可以减少模块的数量UMAP共表达网络
seurat_obj <- RunModuleUMAP(
seurat_obj,
n_hubs = 10, # number of hub genes to include for the UMAP embedding
n_neighbors=15, # neighbors parameter for UMAP
min_dist=0.1# min distance between points in UMAP space
)
# get the hub gene UMAP table from the seurat object
umap_df <- GetModuleUMAP(seurat_obj)
# plot with ggplot
ggplot(umap_df, aes(x=UMAP1, y=UMAP2)) +
geom_point(
color=umap_df$color, # color each point by WGCNA module
size=umap_df$kME*2# size of each point based on intramodular connectivity
) +
umap_theme()
ModuleUMAPPlot(
seurat_obj,
edge.alpha=0.25,
sample_edges=TRUE,
edge_prop=0.1, # proportion of edges to sample (20% here)
label_hubs=2 ,# how many hub genes to plot per module?
keep_grey_edges=FALSE
)
上面分析更加侧重枢纽基因和模块的可视化，而这个方法则把UMAP应用到了hdWGCNA网络中使得图形更加美观。此外，hdWGCNA还提供了更多丰富的可视化教程，其他部分可以点击参考资料中的链接自行查看。
6.富集分析
library(Seurat)
library(tidyverse)
library(cowplot)
library(patchwork)
library(WGCNA)
library(hdWGCNA)
library(enrichR)
library(GeneOverlap)
library(qs)
#dir.create("14-hdWGCNA")
#setwd("14-hdWGCNA")
seurat_obj <- qread("hdWGCNA_object.qs")
# 定义enrichr databases
dbs <- c('GO_Biological_Process_2023',
'GO_Cellular_Component_2023',
'GO_Molecular_Function_2023')
# 富集分析
seurat_obj <- RunEnrichr(
seurat_obj,
dbs=dbs,
max_genes = 100# use max_genes = Inf to choose all genes
)
# 检索输出表
enrich_df <- GetEnrichrTable(seurat_obj)
# 查看结果
head(enrich_df)
# make GO term plots:
EnrichrBarPlot(
seurat_obj,
outdir = "enrichr_plots", # name of output directory
n_terms = 10, # number of enriched terms to show (sometimes more are shown if there are ties)
plot_size = c(5,7), # width, height of the output .pdfs
logscale=TRUE# do you want to show the enrichment as a log scale?
)
# enrichr dotplot
EnrichrDotPlot(
seurat_obj,
mods = "all", # use all modules (default)
database = "GO_Biological_Process_2023", # this must match one of the dbs used previously
n_terms=2, # number of terms per module
term_size=8, # font size for the terms
p_adj = FALSE# show the p-val or adjusted p-val?
) + scale_color_stepsn(colors=rev(viridis::magma(256)))
会输出一个文件夹，里面有很多结果GSEA
https:///Enrichr/#libraries 下载数据
library(fgsea)
# load the GO Biological Pathways file (downloaded from EnrichR website)
pathways <- fgsea::gmtPathways('GO_Biological_Process_2023.txt')
# optionally, remove the GO term ID from the pathway names to make the downstream plots look cleaner
names(pathways) <- stringr::str_replace(names(pathways), " \\s*\\([^\\)]+\\)", "")
# get the modules table and remove grey genes
modules <- GetModules(seurat_obj) %>% subset(module != 'grey')
# rank by Treg_NEW1 genes only by kME
cur_mod <- 'Treg_NEW1'
modules <- GetModules(seurat_obj) %>% subset(module == cur_mod)
cur_genes <- modules[,(c('gene_name', 'module', paste0('kME_', cur_mod)))]
ranks <- cur_genes$kME; names(ranks) <- cur_genes$gene_name
ranks <- ranks[order(ranks)]
# run fgsea to compute enrichments
gsea_df2 <- fgsea::fgsea(
pathways = pathways,
stats = ranks,
minSize = 3,
maxSize = 500
)
# 可视化
top_pathways <- gsea_df2 %>%
subset(pval < 0.05) %>%
slice_max(order_by=NES, n=25) %>%
.$pathway
plotGseaTable(
pathways[top_pathways],
ranks,
gsea_df2,
gseaParam=0.5,
colwidths = c(10, 4, 1, 1, 1)
)
# name of the pathway to plot
selected_pathway <- 'Cellular Respiration'
plotEnrichment(
pathways[[selected_pathway]],
ranks
) + labs(title=selected_pathway)
这里是选择了特定模块的基因进行的分析
7.转录因子分析
library(Seurat)
library(tidyverse)
library(cowplot)
library(patchwork)
library(magrittr)
library(WGCNA)
library(hdWGCNA)
library(igraph)
library(JASPAR2020)
library(JASPAR2024)
library(motifmatchr)
library(TFBSTools)
library(EnsDb.Hsapiens.v86)
library(BSgenome.Hsapiens.UCSC.hg38)
library(GenomicRanges)
library(xgboost)
library(JASPAR2024)
library(RSQLite)
library(EnsDb.Hsapiens.v86)
library(qs)
#dir.create("14-hdWGCNA")
#setwd("14-hdWGCNA")
seurat_obj <- qread("hdWGCNA_object.qs")
# JASPAR 2020
pfm_core <- TFBSTools::getMatrixSet(
x = JASPAR2020,
opts = list(collection = "CORE",
tax_group = 'vertebrates',
all_versions = FALSE)
)
# JASPAR 2024 (not used for this tutorial)
# JASPAR2024 <- JASPAR2024()
# sq24 <- RSQLite::dbConnect(RSQLite::SQLite(), db(JASPAR2024))
# pfm_core <- TFBSTools::getMatrixSet(
# x = sq24,
# opts = list(collection = "CORE", tax_group = 'vertebrates', all_versions = FALSE)
# )
# 进行motif分析
seurat_obj <- MotifScan(
seurat_obj,
species_genome = 'hg38',
pfm = pfm_core,
EnsDb = EnsDb.Hsapiens.v86
)
# 获取motif df:
motif_df <- GetMotifs(seurat_obj)
# 保留所有TFs, 并去除灰色模块基因
tf_genes <- unique(motif_df$gene_name)
modules <- GetModules(seurat_obj)
nongrey_genes <- subset(modules, module != 'grey') %>% .$gene_name
genes_use <- c(tf_genes, nongrey_genes)
# update the gene list and re-run SetDatExpr
seurat_obj <- SetWGCNAGenes(seurat_obj, genes_use)
seurat_obj <- SetDatExpr(seurat_obj, group.by = 'celltype', group_name='Treg')
# define model params:
model_params <- list(
objective = 'reg:squarederror',
max_depth = 1,
eta = 0.1,
nthread=16,
alpha=0.5
)
# 构建转录因子网络
seurat_obj <- ConstructTFNetwork(seurat_obj, model_params)
results <- GetTFNetwork(seurat_obj)
head(results)
# tf gene Gain Cover Frequency Cor
# 1 ZKSCAN1 FOXD1 0.11708119 0.04695019 0.04695019 -0.19391353
# 2 NFIL3 FOXD1 0.10931756 0.05379589 0.05379589 0.19738680
# 3 ZNF652 FOXD1 0.09789632 0.07897635 0.07897635 -0.17072678
# 4 NR4A1 FOXD1 0.09624640 0.04498028 0.04498028 0.18091337
# 5 ZNF24 FOXD1 0.05250378 0.02255133 0.02255133 -0.09566174
# 6 NFKB2 FOXD1 0.05148296 0.06481095 0.06481095 0.16717336
# 策略“A”为每个基因选择前10个TF
# 策略“B”选择每个转录因子的顶级基因
# 策略“C”保留所有高于一定调控分数的TF-基因对
seurat_obj <- AssignTFRegulons(
seurat_obj,
strategy = "A", # 还有B和C策略
reg_thresh = 0.01,
n_tfs = 10
)
# 可视化
# 根据基因表达分为与TF正相关(右侧)或负相关(左侧)的目标基因
p1 <- RegulonBarPlot(seurat_obj, selected_tf='ZNF652')
p2 <- RegulonBarPlot(seurat_obj, selected_tf='NFKB2', cutoff=0.15)
p1 | p2
可以探索与转录因子相关的不同下游靶基因。
获取正/负共表达的TF-基因对
# 正向regulons
seurat_obj <- RegulonScores(
seurat_obj,
target_type = 'positive',
ncores=8
)
# 负向regulons
seurat_obj <- RegulonScores(
seurat_obj,
target_type = 'negative',
cor_thresh = -0.05,
ncores=8
)
# 获取数据结果:
pos_regulon_scores <- GetRegulonScores(seurat_obj, target_type='positive')
neg_regulon_scores <- GetRegulonScores(seurat_obj, target_type='negative')
# 选择感兴趣的TF
cur_tf <- 'FOXP3'
# 把regulon分数添加到Seurat metadata
seurat_obj$pos_regulon_score <- pos_regulon_scores[,cur_tf]
seurat_obj$neg_regulon_score <- neg_regulon_scores[,cur_tf]
# plot using FeaturePlot
p1 <- FeaturePlot(seurat_obj, feature=cur_tf) + umap_theme()
p2 <- FeaturePlot(seurat_obj, feature='pos_regulon_score', cols=c('lightgrey', 'red')) + umap_theme()
p3 <- FeaturePlot(seurat_obj, feature='neg_regulon_score', cols=c('lightgrey', 'seagreen')) + umap_theme()
p1 | p2 | p3
FOXP3主要在Treg细胞中表达，其正相关得分在Treg细胞中也表现为较高的表达水平。
TFNetworkPlot
需要较高配置电脑才能绘制该图
# select TF of interest
cur_tf <- 'FOXP3'
# plot with default settings
p <- TFNetworkPlot(seurat_obj, selected_tfs=cur_tf)
p
# plot the FOXP3 network with primary, secondary, and tertiary targets
p1 <- TFNetworkPlot(seurat_obj, selected_tfs=cur_tf, depth=1, no_labels=TRUE)
p2 <- TFNetworkPlot(seurat_obj, selected_tfs=cur_tf, depth=2, no_labels=TRUE)
p3 <- TFNetworkPlot(seurat_obj, selected_tfs=cur_tf, depth=3, no_labels=TRUE)
p1 | p2 | p3
此外，若数据中含有很多临床信息可以进行模块-性状相关性分析,也可以进行PPI网络分析。
本次分析完成了hdWGCNA的完整实践流程。基于hdWGCNA分析，研究者可以识别与感兴趣细胞亚群相关的模块基因，并进一步围绕这些模块基因开展富集分析、转录因子分析等多种下游探索。
参考资料：
- hdWGCNA identifies co-expression networks in high-dimensional transcriptomics data. Cell Rep Methods. 2023 Jun 12;3(6):100498.
- hdWGCNA：https://smorabit./hdWGCNA/
- hdWGCNA的Network可视化流程：https://smorabit./hdWGCNA/articles/network_visualizations.html
- hdWGCNA的其他可视化流程：https://smorabit./hdWGCNA/articles/hdWGCNA.html
- 生信技能树：https://mp.weixin.qq.com/s/OBvS0I7IUuwcChGoaJVgtQ https://mp.weixin.qq.com/s/3w2e6LwRowMsm2ZzyMAvdA
- 生信菜鸟团：https://mp.weixin.qq.com/s/pLIGZhopfi-NCNtSuvwOBg
- 生信星球：https://mp.weixin.qq.com/s/VlACgl6xlQ9YON4153Rp9g
- 生信小白要知道：https://mp.weixin.qq.com/s/fCvLizKQNWDQKeWBuSG3UQ
