【单细胞转录组实战】七、scRNAseq、M3Drop Pac

2022-09-02 本文已影响0人佳奥

这里是佳奥！我们补充了相关知识，现在进一步学习相关R包的使用方法。

文章数据：Low-coverage single-cell mRNA sequencing reveals cellular heterogeneity and activated signaling pathways in developing cerebral cortex

简要的流程在Git Hub下载的压缩包中，section02-scRNA文件夹的study_scRNAseq文件。

QQ截图20220901134733.png

打开R Project文件，我们正式开始吧。

注：从这里开始，难度会明显提高，先从运行代码开始，再通过帮助文档和参考代码逐渐理解，最后理解下列代码。这需要时间。

1 R包安装与镜像设置

rm(list = ls())
Sys.setenv(R_MAX_NUM_DLLS=999) ##Sys.setenv修改环境设置，R的namespace是有上限的，如果导入包时超过这个上次就会报错,R_MAX_NUM_DLLS可以修改这个上限
options(stringsAsFactors = F) ##options:允许用户对工作空间进行全局设置，stringsAsFactors防止R自动把字符串string的列辨认成factor

options()$repos  ## 查看使用install.packages安装时的默认镜像
options()$BioC_mirror ##查看使用bioconductor的默认镜像
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") ##指定镜像，这个是中国科技大学镜像
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")) ##指定install.packages安装镜像，这个是清华镜像
options()$repos 
options()$BioC_mirror

# http://www.bio-info-trainee.com/3727.html 周末生物信息学培训班准备工作
 
if (!requireNamespace("BiocManager", quietly = TRUE)) 
  install.packages("BiocManager") ##判断是否存在BiocManager包，不存在的话安装

if(!require('Seurat')){
  BiocManager::install('Seurat',ask = F,update = F)
}
if(!require('scran')){
  BiocManager::install(c( 'scran'),ask = F,update = F)
}
if(!require('monocle')){
  BiocManager::install(c( 'monocle'),ask = F,update = F)
}

## http://cole-trapnell-lab.github.io/monocle-release/monocle3/ 
BiocManager::install('destiny')
BiocManager::install(c( 'flexclust','mcclust'),ask = F,update = F)

BiocManager::install(c( 'scRNAseq'),ask = F,update = F)
BiocManager::install(c( 'dbscan'),ask = F,update = F)
BiocManager::install(c( 'M3Drop'),ask = F,update = F)

rm(list = ls()) # clear the environment
#load all the necessary libraries
options(warn=-1) # turn off warning message globally
suppressMessages(library(reticulate))
suppressMessages(library(devtools))
suppressMessages(library(monocle))
suppressMessages(library(flexclust))
suppressMessages(library(mcclust))

全部不报错后，我们开始第一个包的学习。

2 scRNAseq Package

需要自行下载安装一些必要的R包，主要是scRNAseq，然后是其它一些辅助包用来探索这个数据集。

if (!requireNamespace("Rtsne"))
    install.packages("Rtsne")
if (!requireNamespace("FactoMineR"))
    install.packages("FactoMineR")
if (!requireNamespace("factoextra"))
    install.packages("factoextra")
if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
if (!requireNamespace("scater"))
    BiocManager::install("scater")
if (!requireNamespace("scRNAseq"))
    BiocManager::install("scRNAseq") 
if (!requireNamespace("M3Drop"))
    BiocManager::install("M3Drop") 
if (!requireNamespace("ROCR"))
    BiocManager::install("ROCR") 

##加载R包，前提是已经成功安装了
rm(list = ls()) # clear the environment
#load all the necessary libraries
options(warn=-1) # turn off warning message globally
suppressMessages(library(scater))
suppressMessages(library(scRNAseq))
library(ggplot2)
library(tidyr)
library(cowplot)
library("FactoMineR")
library("factoextra")
library("ROCR")

scRNAseq R包中的数据集

这个包内置的是 Pollen et al. 2014 数据集，人类单细胞细胞，分成4类，分别是 pluripotent stem cells 分化而成的 neural progenitor cells (“NPC”) ，还有 “GW16” and “GW21” ，“GW21+3” 这种孕期细胞，理解这些需要一定的生物学背景知识，如果不感兴趣，可以略过。

不过本例子只使用了数据集的4种细胞类型而已，因为 scRNAseq 这个R包就提供了这些，完整的数据是 23730 features，301 samples, 地址为：

https://hemberg-lab.github.io/scRNA.seq.datasets/human/tissues/

这里面的表达矩阵是由 RSEM (Li and Dewey 2011) 软件根据 hg38 RefSeq transcriptome 得到的，总是130个文库，每个细胞测了两次，测序深度不一样。

library(scRNAseq)
## ----- Load Example Data -----
data(fluidigm) 
fluidigm<-ReprocessedFluidigmData()##作者修改了提取数据的方法，需要多运行这一步
ct <- floor(assays(fluidigm)$rsem_counts)##counts含有小数，因此需要floor()
ct[1:4,1:4] 
sample_ann <- as.data.frame(colData(fluidigm))
DT::datatable(sample_ann)

先探索表型信息

前面说到，这个数据集是130个文库，每个细胞测了两次，测序深度不一样，这130个细胞，分成4类，分别是： pluripotent stem cells 分化而成的 neural progenitor cells (“NPC”) ，还有 “GW16” and “GW21” ，“GW21+3” 这种孕期细胞。

批量，粗略的看一看各个细胞的一些统计学指标的分布情况：

```{r fig.width=10, fig.height=15}##rmarkdown画图技巧
box <- lapply(colnames(sample_ann[,1:19]),function(i) {
    dat <-  sample_ann[,i,drop=F] 
    dat$sample=rownames(dat)
    ## 画boxplot 
   ggplot(dat, aes('all cells', get(i))) +
          geom_boxplot() +
          xlab(NULL)+ylab(i)
})
plot_grid(plotlist=box, ncol=5 )
# ggsave(file="stat_all_cells.pdf")

QQ截图20220901145444.png

因为进行了简单探索，对表型数据就有了把握，接下来可以进行一定程度的过滤，因为细节太多，这里重点是批量，代码技巧非常值得学习。教程略过了这一部分。

pa <- colnames(sample_ann[,c(1:9,11:16,18,19)])
tf <- lapply(pa,function(i) {
 # i=pa[1]
  dat <-  sample_ann[,i]  
  dat <- abs(log10(dat))
  fivenum(dat)
  (up <- mean(dat)+2*sd(dat))
  (down <- mean(dat)- 2*sd(dat) ) 
  valid <- ifelse(dat > down & dat < up, 1,0 ) 
})

tf <- do.call(cbind,tf)
choosed_cells <- apply(tf,1,function(x) all(x==1))
table(sample_ann$Biological_Condition)
sample_ann=sample_ann[choosed_cells,]
table(sample_ann$Biological_Condition)
ct <- ct[,choosed_cells]

再探索基因表达情况

ct[1:4,1:4] 
counts <- ct
fivenum(apply(counts,1,function(x) sum(x>0) ))
boxplot(apply(counts,1,function(x) sum(x>0) ))
fivenum(apply(counts,2,function(x) sum(x>0) ))
hist(apply(counts,2,function(x) sum(x>0) ))
choosed_genes=apply(counts,1,function(x) sum(x>0) )>0
table(choosed_genes)
counts <- counts[choosed_genes,]

##查看一下表达不为0的基因个数
> table((apply(counts,1,function(x) sum(x>0)>0 )))

 TRUE 
17096

几乎没有基因在所有细胞都表达，中位线为3_{5，如果基因在3}5个细胞内不表达即可过滤掉（阈值）。

QQ截图20220901150747.png
接下来要利用自己的常规转录组数据分析知识。

细胞之间的所有的基因的表达量的相关性

下面的计算，都是基于log后的表达矩阵。

dat <- log2(edgeR::cpm(counts) + 1)
dat[1:4, 1:4]
dat_back <- dat

先备份这个表达矩阵，后面的分析都用得上。

exprSet <- dat_back
colnames(exprSet)
pheatmap::pheatmap(cor(exprSet))
group_list <- sample_ann$Biological_Condition
tmp <- data.frame(g = group_list)
rownames(tmp) <- colnames(exprSet)

##组内的样本的相似性应该是要高于组间的！
pheatmap::pheatmap(cor(exprSet), annotation_col = tmp)
dim(exprSet)
exprSet = exprSet[apply(exprSet, 1, function(x)
sum(x > 1) > 5), ]
dim(exprSet)
 
dim(exprSet)
exprSet <-  exprSet[names(sort(apply(exprSet, 1, mad), decreasing = T)[1:500]), ]
dim(exprSet)
M <-cor(log2(exprSet + 1))
tmp <- data.frame(g = group_list)
rownames(tmp) <-  colnames(M)
pheatmap::pheatmap(M, annotation_col = tmp)

table(sample_ann$LibraryName)

QQ截图20220901152639.png

可以看到,从细胞的相关性角度来看，到NPC跟另外的GW细胞群可以区分的很好，但是GW本身的3个小群体并没有那么好的区分度。

而且简单选取top的sd的基因来计算相关性，并没有很明显的改善。

但是可以看到每个细胞测了两次，所以它们的相关性要好于其它同类型的细胞。

对表达矩阵进行简单的层次聚类

如果计算机资源不够，这里可以先对基因进行一定程度的挑选，最简单的就是选取top的sd的基因，这里略。

dat <- dat_back
hc <- hclust(dist(t(dat))) 
plot(hc,labels = FALSE)
clus <-  cutree(hc, 4) #对hclust()函数的聚类结果进行剪枝，即选择输出指定类别数的系谱聚类结果。
group_list <-  as.factor(clus) ##转换为因子属性
table(group_list) ##统计频数
table(group_list,sample_ann$Biological_Condition)

QQ截图20220901152800.png

可以看到GW16和GW21是很难区分开来的，如果是普通的层次聚类的话。

然后看看最常规的PCA降维结果

降维算法很多，详情可以去自行搜索学习，比如：

主成分分析PCA
多维缩放(MDS)
线性判别分析(LDA)
等度量映射(Isomap)
局部线性嵌入(LLE)
t-SNE
Deep Autoencoder Networks

这里只介绍 PCA 和 t-SNE

dat <- dat_back
dat <- t(dat)
dat <- as.data.frame(dat)
plate <- sample_ann$Biological_Condition # 这里定义分组信息
dat <-  cbind(dat, plate) #cbind根据列进行合并，即叠加所有列 #矩阵添加批次信息
dat[1:4, 1:4]
table(dat$plate)

# The variable plate (index = ) is removed
# before PCA analysis
dat.pca <- PCA(dat[, -ncol(dat)], graph = FALSE)
head(dat.pca$var$coord) ## 每个主成分的基因重要性占比
head(dat.pca$ind$coord) ## 每个细胞的前5个主成分取值。
fviz_pca_ind(
      dat.pca,
      #repel =T,
      geom.ind = "point",
      # show points only (nbut not "text")
      col.ind = dat$plate,
      # color by groups
      #palette = c("#00AFBB", "#E7B800"),
      addEllipses = TRUE,
      # Concentration ellipses
      legend.title = "Groups"
)

同样的，很明显可以看到NPC跟另外的GW细胞群可以区分的很好，但是GW本身的3个小群体并没有那么好的区分度。

接着是稍微高大上的tSNE降维

因为计算量的问题，这里先选取PCA后的主成分，然进行tSNE，当然，也有其它做法，比如选取变化高的基因，显著差异基因等等。

# 选取前面PCA分析的5个主成分。
dat_matrix <- dat.pca$ind$coord
# Set a seed if you want reproducible results
set.seed(42)
library(Rtsne) 
# 如果使用原始表达矩阵进行 tSNE耗时很可怕，dat_matrix = dat_back
# 出现Remove duplicates before running TSNE 则check_duplicated = FALSE
# tsne_out <- Rtsne(dat_matrix,pca=FALSE,perplexity=30,theta=0.0, check_duplicates = FALSE) # Run TSNE
tsne_out <- Rtsne(dat_matrix,perplexity=10)
plate <- sample_ann$Biological_Condition # 这里定义分组信息
plot(tsne_out$Y,col= rainbow(4)[as.numeric(as.factor(plate))], pch=19)

对PCA或者tSNE结果进行kmeans或者dbscan算法聚类

降维是降维，聚类是聚类，需要理解其中的区别。

降维与否，不同的降维算法选择，不同参数的选择得到的结果都不一样。

聚类也是一样，不同的算法，不同的参数。

# 前面我们的层次聚类是针对全部表达矩阵，这里我们为了节省计算量，可以使用tsne_out$Y这个结果
head(tsne_out$Y)
opt_tsne=tsne_out$Y
table(kmeans(opt_tsne,centers = 4)$clust)
plot(opt_tsne,  col=kmeans(opt_tsne,centers = 4)$clust, pch=19, xlab="tSNE dim 1", ylab="tSNE dim 2")
library(dbscan)
plot(opt_tsne,  col=dbscan(opt_tsne,eps=3.1)$cluster, pch=19, xlab="tSNE dim 1", ylab="tSNE dim 2")
table(dbscan(opt_tsne,eps=3.1)$cluster)
# 比较两个聚类算法区别
table(kmeans(opt_tsne,centers = 4)$clust,dbscan(opt_tsne,eps=3.1)$cluster)

3 M3Drop Package分析scRNA数据

首先构建M3Drop需要的对象

library(M3Drop) 
Normalized_data <- M3DropCleanData(counts, 
                                   labels = sample_ann$Biological_Condition , 
                                   is.counts=TRUE, min_detected_genes=2000)

##查看对象
dim(Normalized_data$data)
length(Normalized_data$labels)
class(Normalized_data)
str(Normalized_data)

这个包设计比较简单，并没有构建S4对象，只是一个简单的list而已。

统计学算法 Michaelis-Menten

需要深入读该文章，了解其算法，这里略过，总之它对单细胞转录组的表达矩阵进行了一系列的统计检验。

fits <- M3DropDropoutModels(Normalized_data$data)

# Sum absolute residuals
data.frame(MM=fits$MMFit$SAr, Logistic=fits$LogiFit$SAr,
           DoubleExpo=fits$ExpoFit$SAr) 
# Sum squared residuals
data.frame(MM=fits$MMFit$SSr, Logistic=fits$LogiFit$SSr,
           DoubleExpo=fits$ExpoFit$SSr)

QQ截图20220901154758.png

找差异基因

DE_genes <- M3DropFeatureSelection(Normalized_data$data, 
                                   mt_method="fdr", mt_threshold=0.01)
dim(DE_genes)
head(DE_genes)

这里是针对上面的统计结果来的。

QQ截图20220901155120.png

针对差异基因画热图

par(mar=c(1,1,1,1)) 
heat_out <- M3DropExpressionHeatmap(DE_genes$Gene, Normalized_data$data, 
                                    cell_labels = Normalized_data$labels)

可视化了解一下找到的差异基因在不同的细胞类型的表达分布情况。

QQ截图20220901155236.png

聚类

这里可以重新聚类后，针对自己找到的类别来分别找marker基因，不需要使用测试数据自带的表型信息。

cell_populations <- M3DropGetHeatmapCellClusters(heat_out, k=4)
library("ROCR") 
marker_genes <- M3DropGetMarkers(Normalized_data$data, cell_populations)
table(cell_populations,Normalized_data$labels)

每个类别的marker genes

head(marker_genes[marker_genes$Group==4,],20) 
marker_genes[rownames(marker_genes)=="FOS",]

也可以针对这些 marker genes去画热图，当然，得根据AUC和P值来挑选符合要求的差异基因去绘图。

par(mar=c(1,1,1,1)) 
choosed_marker_genes=as.character(unlist(lapply(split(marker_genes,marker_genes$Group), function(x) (rownames(head(x,20))))))
heat_out <- M3DropExpressionHeatmap(choosed_marker_genes, Normalized_data$data, cell_labels =  cell_populations)

如果遇到Error in plot.new() : figure margins too large报错，则单独将heat_out这行命令复制出来运行

对感兴趣基因集进行注释

通常是GO/KEGG等数据库，通常是超几何分布，GSEA，GSVA等算法。

这里就略过。

显示运行环境

sessionInfo()

下一篇我们继续用单细胞的R包来处理数据集。

加油！我们下一篇再见！

【单细胞转录组实战】七、scRNAseq、M3Drop Pac

注：从这里开始，难度会明显提高，先从运行代码开始，再通过帮助文档和参考代码逐渐理解，最后理解下列代码。这需要时间。

1 R包安装与镜像设置

2 scRNAseq Package

scRNAseq R包中的数据集

先探索表型信息

再探索基因表达情况

细胞之间的所有的基因的表达量的相关性

对表达矩阵进行简单的层次聚类

然后看看最常规的PCA降维结果

接着是稍微高大上的tSNE降维

对PCA或者tSNE结果进行kmeans或者dbscan算法聚类

3 M3Drop Package分析scRNA数据

首先构建M3Drop需要的对象

统计学算法 Michaelis-Menten

找差异基因

针对差异基因画热图

聚类

每个类别的marker genes

对感兴趣基因集进行注释

显示运行环境

猜你喜欢

热点阅读

【单细胞转录组 实战】七、scRNAseq、M3Drop Pac

注：从这里开始，难度会明显提高，先从运行代码开始，再通过帮助文档和参考代码逐渐理解，最后理解下列代码。这需要时间。

1 R包安装与镜像设置

2 scRNAseq Package

scRNAseq R包中的数据集

先探索表型信息

再探索基因表达情况

细胞之间的所有的基因的表达量的相关性

对表达矩阵进行简单的层次聚类

然后看看最常规的PCA降维结果

接着是稍微高大上的tSNE降维

对PCA或者tSNE结果进行kmeans或者dbscan算法聚类

3 M3Drop Package分析scRNA数据

首先构建M3Drop需要的对象

统计学算法 Michaelis-Menten

找差异基因

针对差异基因画热图

聚类

每个类别的marker genes

对感兴趣基因集进行注释

显示运行环境

猜你喜欢

热点阅读

【单细胞转录组实战】七、scRNAseq、M3Drop Pac