学习单细胞WGCNA（一）

学生信的或多或少都听说过WGCNA，其实就是一种相关性分析，你设置一种表型，然后就能找到几群相关表达的基因，具体内容就不啰嗦了，网上介绍的资料很多。简书上也有不少WGCNA的代码，但是按照别人的代码跑下来总会出现一些未定义的变量。这边为了减少大家的学习负担，我把我的代码共享出来，今天主要介绍单细胞表达矩阵的处理，也不一定正确，希望大大多多批评指正。

本人生信小白，参加了生信星球学习小组，0基础学习R语言。今天也是第一次在简书上写文章。

进入正题

安装R包

library(WGCNA)
library(Seurat)
library(tidyverse)
library(reshape2)
library(stringr)

GEO随便找一个单细胞的数据集

我这边网上找了血管的AS数据集，提取了巨噬细胞分析，首先看一下细胞分群和每群的细胞数量

> table(pbmc$celltype)
    Inf_MP   ISGhi_MP     Res_MP TREM2hi_MP 
       702        762        418       1289

由于单细胞矩阵的稀疏性（存在很多0），这样会导致相关性分析结果不好-相关性很低。因此我们需要合并一下细胞，把几个临近相似的细胞合并为一个样本。

构建表达矩阵

> datadf <- as.matrix(pbmc@assays$RNA@data )
> idd1 <- pbmc@meta.data
> Inter.id1<-cbind(rownames(idd1),idd1$celltype)
> rownames(Inter.id1)<-rownames(idd1)
> colnames(Inter.id1)<-c("CellID","celltype")
> Inter.id1<-as.data.frame(Inter.id1)
> head(Inter.id1)
                                           CellID celltype
ZsLdlr0_AACGTTGAGACGCTTT ZsLdlr0_AACGTTGAGACGCTTT   Res_MP
ZsLdlr0_AACGTTGTCAAACGGG ZsLdlr0_AACGTTGTCAAACGGG   Res_MP
ZsLdlr0_AAGACCTCAGCTGCTG ZsLdlr0_AAGACCTCAGCTGCTG   Res_MP
ZsLdlr0_ACAGCTACAACACCCG ZsLdlr0_ACAGCTACAACACCCG   Res_MP
ZsLdlr0_ACATGGTAGGAGTACC ZsLdlr0_ACATGGTAGGAGTACC ISGhi_MP
ZsLdlr0_ACATGGTAGTTGTCGT ZsLdlr0_ACATGGTAGTTGTCGT   Res_MP
> Inter1<-datadf[,Inter.id1$CellID]
> Inter2<-as.matrix(Inter1)
> Inter2[1:4,1:4]
       ZsLdlr0_AACGTTGAGACGCTTT ZsLdlr0_AACGTTGTCAAACGGG ZsLdlr0_AAGACCTCAGCTGCTG ZsLdlr0_ACAGCTACAACACCCG
Sox17                         0                        0                        0                 0.000000
Mrpl15                        0                        0                        0                 0.000000
Lypla1                        0                        0                        0                 0.000000
Tcea1                         0                        0                        0                 1.044136

开始合并细胞

> pseudocell.size = 10 ## 10个细胞合并为一个
> new_ids_list1 = list()
> length(levels(as.factor(Inter.id1$celltype)))
[1] 4
> #pseudo cell-合并细胞
> for (i in 1:length(levels(as.factor(Inter.id1$celltype)))) {
+   cluster_id = levels(as.factor(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_list1长这个样子

image.png

此时每10个细胞共享一个名字

接下来无脑运行

new_ids <- as.data.frame(new_ids)
head(new_ids)
new_ids_length <- table(new_ids)
new_ids_length

new_colnames <- rownames(new_ids)  ###add
#rm(all.data1)
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)]
result <- t(new.data)
dim(result)
idd1$pseudocell <- new_ids$new_ids[match(rownames(pbmc@meta.data),rownames(new_ids))]
head(idd1)
head(new_ids)
> dim(idd1)
[1] 3171   33
> 
> new_idd1 <- idd1[!duplicated(idd1$pseudocell),]
> dim(new_idd1)
[1] 319  33

可以看到原来3171个细胞，现在还剩319个，WGCNA每组样本量不能少于15个，大家到时候注意一下。

过滤基因

###过滤基因-只保留高变基因####
pbmc <- FindVariableFeatures(pbmc,nfeatures = 2000)
result[1:4, 1:4]
colnames(result)
Cluster1 <- result[intersect(Seurat::VariableFeatures(pbmc),rownames(result)),]
#最终我们的表达谱是这样的：
Cluster1[1:4,1:4]
#       Inf_MP_Cell0 Inf_MP_Cell1 Inf_MP_Cell10 Inf_MP_Cell11
#Retnla    0.0000000    0.0000000    0.00000000     0.0000000
#Ccl8      0.0848587    0.3074741    0.27450422     0.5210496
#Mfge8     0.2808082    0.3869427    0.09270133     0.0000000
#S100a9    0.0000000    0.0000000    0.00000000     0.0000000
dim(Cluster1)
#[1] 1949  319

好了，表达矩阵就处理完了，下次讲WGCNA步骤了，我也要继续去学习了。