R语言实现单细胞测序的表达矩阵复原

一粒沙

发布于 2019-09-12 15:28:33

3K0

发布于 2019-09-12 15:28:33

单细胞测序成为当下热门领域，很多新的观点以及新的机制通过单细胞测序得到确认以及放大。但是单细胞测序同时存在很大的噪音，因此如何校正单细胞测序所带来的噪音也成为研究的重点。故学者们开发了很多相应的算法去解决这个问题例如：MAGIC，scImpute等，今天给大家介绍另外一个在R语言中实现的算法SAVER，于2018年发表在nature method：

此包的安装很简单那就是：

Install.packages(“SAVER”)

接下来就是其使用，因为其包含的函数很少，我们主要看下主函数saver的相关系数：

主要的参数：

Do.fast 指是否展示计算过程中的参数显示，默认是TRUE。

Npred 指调用的进程数。

Pred.cells,pred.genes 指需要恢复的特定样本或者特定基因。

Pred.genes.only 指是否对仅特定的基因进行评估。

Estimates.only 指只是显示评估后的表达矩阵。

接下来我们就直接进行实例讲解：

数据的载入，此包的数据输入主要是各基因的reads数：

首先要下载所需要的示例数据，下载地址：https://storage.googleapis.com/linnarsson-lab-www-blobs/blobs/cortex/expression_mRNA_17-Aug-2014.txt

接下来就是数据的载入：

library(SAVER)
data.path <-"expression_mRNA_17-Aug-2014.txt"
# Need to remove first 10 rows of metadata
raw.data <- read.table(data.path, header= FALSE, skip = 11, row.names = 1, check.names = FALSE)
cortex <- as.matrix(raw.data[, -1])
cellnames <- read.table(data.path, skip= 7, nrows = 1, row.names = 1, stringsAsFactors = FALSE)
colnames(cortex) <- cellnames[-1]

然后是模型的构建：

我们这次选择的运行环境是ubuntu系统下的3.6.1版本的R语言：

cortex.saver <- saver(cortex, ncores =12)

下图是计算所调用的进程数：

最后我们看下结果的信息：

str(cortex.saver)

上图所包含的主要参数信息如下：

当然还可以对指定的基因集进行评估：

genes <- c("Thy1","Mbp", "Stim2", "Psmc6", "Rps19")
genes.ind <- which(rownames(cortex) %in%genes)
 
# Generate predictions for those genes andreturn entire dataset
cortex.saver.genes <- saver(cortex,pred.genes = genes.ind,
                            estimates.only =TRUE)

cortex.saver.genes.only <- saver(cortex,pred.genes = genes.ind,
                                pred.genes.only = TRUE, estimates.only = TRUE)

如果基因数很多，但是呢计算机内存不够的话，此包还提供了一个功能那就是分离后再合并计算结果：

                                                       
saver1 <- saver(cortex, pred.genes =1:5, pred.genes.only = TRUE,
                do.fast = FALSE)
saver2 <- saver(cortex, pred.genes =6:10, pred.genes.only = TRUE,
                do.fast = FALSE)
saver3 <- saver(cortex, pred.genes =11:15, pred.genes.only = TRUE,
                do.fast = FALSE)
 
 
saver.all <- combine.saver(list(saver1,saver2, saver3))