【原】单细胞水平的基因共表达网络网络分析之scGENA

健明 2022-08-09 发布于广东

展开全文

说起基因共表达网络网络分析大家第一时间想到的应该是weighted gene coexpression network analysis (WGCNA) ，其实还有：

coexpression of RNA-seq data (Coseq),
coexpression modules identification tool (CEMiTool),
petal,
CoXpress,
coexpressed biological processes (CoP)

当然了，并不需要大家都掌握，WGCNA就足够了。我们《生信技能树》这些年有很多关于WGCNA的实战细节建议分享，见：

但是这些算法都是为传统的bulk转录组测序或者转录组表达量芯片数据设计的，但是单细胞表达量矩阵的特点是绝大部分基因在绝大部分细胞里面的表达量都是0，比如以大家熟知的pbmc3k数据集为例。

读取pbmc3k的seurat对象结果，大家先安装这个数据集对应的包，并且对它进行降维聚类分群，参考前面的例子：人人都能学会的单细胞聚类分群注释，而且每个亚群找高表达量基因，都存储为Rdata文件。

library(SeuratData) #加载seurat数据集  
getOption('timeout')
options(timeout=10000)
#InstallData("pbmc3k")  
data("pbmc3k")  
ct = pbmc3k.final@assays$RNA@counts
dim(ct) 
# [1] 13714  2638
sum(as.numeric(ct==0)) / (ncol(ct) * nrow(ct))
# [1] 0.9381182

可以看到这个单细胞数据集，2638个细胞，只有13714个基因，并不是人类的两万个基因。而且可以看到， 93.8%的基因的表达量都是0 ，虽然对很多高表达量基因来说0值不会太多，但是这样的表达量分布情况也让我们很难直接对两个基因做相关性分析，单细胞水平的基因共表达网络网络分析也就无从谈起。

这里我们推荐一个刚刚发表的文献，标题是：《scGENA: A Single-Cell Gene Coexpression Network Analysis Framework for Clustering Cell Types and Revealing Biological Mechanisms》，这个流程相关的代码在 GitHub (https://github.com/zpliulab/scGENA) ，并不是一个R包哦。可以看到主要是5个步骤：

Phase 1 set up and preprocesses the scRNA-seq dataset to filter low-dimensional and noisy single-cell expression genes.
Phase 2 performs a differentially expressed (DE) genes analysis to determine which genes are expressed significantly dif- ferent in different conditions. These genes can reveal biological information about the processes that are influenced by the conditions of interest.
Phase 3 applies the SAVER imputation method to estimate and replace dropout values in each gene cross cell’s actual missing expression level, reducing technical differences while preserving biological vari- ability across cells [23].
Phase 4 constructs a coexpression network analysis to shed light on the transcriptional regulatory mechanisms underpinning numerous biological processes [24].
Phase 5 performs further analyses, including a functional enrichment analysis, differential coexpression network analysis, and overlapping genes identification across different cell-types to better interpret the biological insights.

可以看到，其实就是最开始调用seurat里面的质量控制步骤去除低质量的细胞和不重要的基因，然后根据细胞的生物学来源分组后做差异分析，仍然是seurat里面的函数即可，依赖于MAST (model-based analysis of single-cell transcriptomics) 包，接着就是对有统计学意义的差异基因进行插值（因为单细胞技术限制很多基因被错误的0值化了，插值相当于纠正那些0值），这里作者选择了 SAVER method in this pipeline because it imputed original zero values to actual values ，表达量矩阵到这里，就基本上很少0值了，这3个步骤的效果如下所示：

解决单细胞0值问题

然后，“精彩的”地方来了，其实作者仍然是使用了WGCNA算法做单细胞水平的基因共表达网络网络分析，

流程图如下所示：

流程图

这个流程相关的代码在 GitHub (https://github.com/zpliulab/scGENA) ，并不是一个R包哦。我打开看了看：

  103 1_Data_preparation.R
  184 2_Data_preprocessing.R
  147 3_DEgenes.R
   71 4_imputation.R
  269 5_co_ex_net.R
  117 further_analysis.R

代码量并不大，而且很容易看懂哦，好像是一个博士毕业课题。

感兴趣的小伙伴可以拿作者代码里面的示例数据试试看，就是GSE81608，是Xin et al., 2016文章的，基于 Fluidigm C1 仪器的SMART-seq单细胞哦，并不是常规的10x单细胞。因为10x单细胞的0值更恐怖，所以这个流程能否在10x单细胞表达量矩阵应用或者说效果怎么样，都不好说哦。

写在文末

我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的，有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 jmzeng1314@163.com

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。