单细胞实战之单细胞hdWGCNA分析——入门到进阶(高级篇3）

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表达矩阵的稀疏性大大降低，更适合用于后续分析。

# 各组构建metacellseurat_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 patchworkwrap_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))# 计算完整单细胞数据集中的所有MEsseurat_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个枢纽基因的表格,这里选择了10hub_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 patchworkwrap_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 patchworkwrap_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 objectMEs <- 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 functionp <- 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 outputp

图片

单模块的网络图

# 使用ModuleNetworkPlot可视化每个模块前50(数值可自定)的hub geneModuleNetworkPlot(    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 networkHubGeneNetworkPlot(  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 objectumap_df <- GetModuleUMAP(seurat_obj)# plot with ggplotggplot(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 databasesdbs <- 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 dotplotEnrichrDotPlot(  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://maayanlab.cloud/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 cleanernames(pathways) <- stringr::str_replace(names(pathways), " \\s*\\([^\\)]+\\)", "")# get the modules table and remove grey genesmodules <- GetModules(seurat_obj) %>% subset(module != 'grey')# rank by Treg_NEW1 genes only by kMEcur_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_nameranks <- ranks[order(ranks)]# run fgsea to compute enrichmentsgsea_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) %>%     .$pathwayplotGseaTable(    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 2020pfm_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_namegenes_use <- c(tf_genes, nongrey_genes)# update the gene list and re-run SetDatExprseurat_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-基因对

# 正向regulonsseurat_obj <- RegulonScores(    seurat_obj,    target_type = 'positive',    ncores=8)# 负向regulonsseurat_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')# 选择感兴趣的TFcur_tf <- 'FOXP3'# 把regulon分数添加到Seurat metadataseurat_obj$pos_regulon_score <- pos_regulon_scores[,cur_tf]seurat_obj$neg_regulon_score <- neg_regulon_scores[,cur_tf]# plot using FeaturePlotp1 <- 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 interestcur_tf <- 'FOXP3'# plot with default settingsp <- TFNetworkPlot(seurat_obj, selected_tfs=cur_tf)p# plot the FOXP3 network with primary, secondary, and tertiary targetsp1 <- 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.github.io/hdWGCNA/hdWGCNA的Network可视化流程：https://smorabit.github.io/hdWGCNA/articles/network_visualizations.htmlhdWGCNA的其他可视化流程：https://smorabit.github.io/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 本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。