WGCNA复现：小胶质细胞亚群在脑发育时髓鞘形成的作用

最新推荐文章于 2024-10-10 11:25:11 发布

原创

最新推荐文章于 2024-10-10 11:25:11 发布 · 3.8k 阅读

24 ·

CC 4.0 BY-SA版权

文章标签：

#r语言 #大数据 #数据挖掘

本文详细介绍了使用WGCNA分析小胶质细胞在大脑发育中作用的全过程，从数据预处理到模块分析，再到GO富集分析，最后导出感兴趣模块进行可视化。通过实战案例，深入理解WGCNA在生物信息学中的应用。

Hi大家好，我是老米！生信入门一个月，感谢生信技能树平台。

Homework Response to: 重复一篇WGCNA分析的文章。WGCNA学习教程：一文学会WGCNA分析。

复现的文献：A novel microglial subset plays a key role in myelinogenesis in developing brain。EMBO J，10多分的好杂志。

课题背景

小胶质细胞是神经系统中的巨噬细胞，可介导脑稳态和神经炎症。作用至关重要，但在大脑发育中具体机制尚未完全阐明。这篇文章中，通过对比健康成人和炎症激活(构建的实验性自身免疫性髓鞘炎模型EAE)的细胞群，发现新生鼠的小胶质细胞对于髓鞘生成及神经形成起关键作用。其中CD11c +小胶质细胞亚群在发育中大脑的髓鞘区域高表达，介导神经和胶质细胞的存活，迁移和分化。这些细胞是IGF1的主要来源，而选择性耗竭会该亚群可导致性髓鞘形成受损。引用Jimmy的话：

值得一提的是，在单细胞水平研究小神经胶质细胞（Microglia）动态发育和异质性已经有了不少研究。

波士顿儿童医院的研究者们分析了超过76,000个来自于发育、衰老和脑部感染后的小鼠脑部的小胶质细胞，结果表明至少有9种转录特异的小胶质细胞形态，它们可以表达特定的基因集，且位于特定的脑区。发表于免疫学杂志Immunity， doi:10.1016/j.immuni.2018.11.004 (2019).

斯坦福大学医学院的研究者采用高深度scRNA测序揭示了小胶质细胞和脑髓细胞的发育异质性，发表于Neuron，这些细胞取自于胚胎期、出生后早期和成年的小鼠不同脑区。我们发现大部分的成年小胶质细胞表达稳定的基因（homeostatic genes），且不同脑区间没有差异。相反，出生后早期的小胶质细胞异质性更高。 doi:10.1016/j.neuron.2018.12.006 (2019).

德国弗莱堡大学医学院神经病理学研究所的研究者采用单细胞RNA测序揭示小鼠和人的小神经胶质细胞的空间和时间异质性，成果最近以Letter的形式发表于Nature杂志。doi:10.1038/s41586-019-0924-x (2019).

这些都是当前的脑发育研究的热点和前沿！（Jimmy老师，你知识怎么能这么渊博？我这个搞神经发育的都惭愧【捂脸】）

该文章对五个处理组，共17个样本：

orange represents neonatal CD11c+ microglia (n = 4),

green neonatal CD11c microglia (n = 4),

blue EAE CD11c+ microglia (n = 3),

purple EAE CD11c microglia (n = 3),

black adult microglia (n = 3).
其实就是 neonatal(新生的) 和 EAE的Microglia，还有CD11C阳性和阴性，然后和成年小鼠的Microglia进行比较。

的老鼠进行了RNASeq测序，数据在GEO可供下载：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE78809。当然我偷了个懒，该文章Supplemental material提供了整理之后的csv矩阵，大概1万3千个基因。作者进行了WGCNA分析，以此证明小胶质细胞的作用及寻找关键信号通路及基因（嗯，估计这个实验室后续还有大文章出来）。组学及生信分析的篇幅占了整个文章的一半以上。在此，仅重复文章的Figure3。

跌跌撞撞学了WGCNA半个月，大概了解了这个原理，现在还一知半解，对很多细节摸索了很久，头很大。下面进入正题，请大家一起跟着我头大。

1.数据处理

rm(list = ls())
dev.off()
library(WGCNA)
norm.counts <- read.csv("Dataset_EV2.csv",header = T,sep = ",") ##original counts
rownames(norm.counts)<-norm.counts[,1]
norm.counts <- norm.counts[,2:18]
save(norm.counts,file = "norm.counts.Rdata")
## 处理表型信息，就是老鼠的分组信息
a <- colnames(norm.counts)
library(stringr)
str_tmp=as.data.frame(str_split(a,"[_]",simplify = T))
rownames(str_tmp) <- colnames(norm.counts)
colnames(str_tmp) <- c("group","sampleRep")
datTraits <- str_tmp
datTraits$groupNo <- c(rep(1,3),rep(2,4),rep(3,4),rep(4,3),rep(5,3)) 
###数据初步处理完成

## 筛选中位绝对偏差前75%的基因，取MAD大于0.3
## 筛选后会降低运算量，也会失去部分信息
## 也可不做筛选，使MAD大于0即可
norm.counts <- log2(norm.counts)  ##注意我在这里进行了log2，因为我发现如果不log的话，出现了个很有意思的现象，请您往下看
m.mad <- apply(norm.counts,1,mad)
dataExprMad <- norm.counts[which(m.mad > max(quantile(m.mad, probs=seq(0, 1, 0.25))[2],0.3)),] #25%  = 0.22
datExpr0 <- as.data.frame(t(dataExprMad))
##最终取到8222个基因。

###########start handling the data
##看看有没有为空的值，需要剔除。这里下载的都是作者处理好的，基本没问题。
gsg <- goodSamplesGenes(datExpr0,verbose = 3)
gsg$allOK
if (!gsg$allOK){
   
   
  # Optionally, print the gene and sample names that were removed:
  if (sum(!gsg$goodGenes)>0)
    printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes],
                                              collapse = ", ")));
  if (sum(!gsg$goodSamples)>0)
    printFlush(paste("Removing samples:",
                     paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")));
  # Remove the offending genes and samples from the data:
  datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}
gsg <- goodSamplesGenes(datExpr0,verbose = 3)
gsg$allOK

###画个树看看样本有没有离群值！！！
if(T){
   
   
  # Plot a line to show the cut
  sampleTree <- hclust(dist(datExpr0),method = "average")
  par(cex=.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
  )
}

##然后保存数据
dim(datExpr0) ##17，8222
datExpr <- datExpr0
nGenes <- ncol(datExpr)
nSamples <- nrow(datExpr)
save(nGenes,nSamples,datExpr,datTraits,file="input.Rdata")

结果如图：
在这里插入图片描述

这结果挺好的，符合实验设计，都聚类好了。但是！但是如果不对原始数值进行log2处理，会出现一个outlier老鼠，如图：
在这里插入图片描述

看到没，那个新生的NEO1号鼠有点调皮！他超脱于所有老鼠之上，统领其他所有老鼠！这要么是问题少年，要么是智商超群，比如科大少年班的那种！而我大概深究了一下，控制这个离群值的大概是500个基因，哈哈哈。。。

2. MDS分析

专业名词叫多维标度法（Multidimensional Scaling）MDS，是一种维数缩减方法，把高维的数据点映射到一个低维的流形上；同时也是一种可视化方法，实践中通常利用2D或3D的MDS 结果观察（投影后）点的分布和聚集来研究数据的性质。说人话叫做：通过基因表达log2组间差异，看看这17个老鼠是否能够分到一个类别，比如成年鼠归到一类，EAE归到一类，新生鼠归到一类。和上面的有点类似，又不完全一样。

参考学习了stat的视频，也参考了他提供的代码。如下：

rm(list = ls())
load(file = 'input.Rdata')
log2.data.matrix <- as.data.frame(t(datExpr))
## now create an empty distance matrix
log2.distance.matrix <- matrix(0,
                               nrow=ncol(log2.data