单细胞转录组一条龙数据分析流程之popsicleR

生信菜鸟团

发布于 2022-05-23 14:15:21

7340

发布于 2022-05-23 14:15:21

文章被收录于专栏：生信菜鸟团

之前我们在教程：为什么要用conda来安装一个R包 提到过这个popsicleR是一个单细胞转录组数据分析的R包，它蛮简单的，就是降维聚类分群和生物学注释。最基础的往往是降维聚类分群，参考前面的例子：人人都能学会的单细胞聚类分群注释 ，这样的基础认知，也可以看基础10讲：

popsicleR这个单细胞转录组数据分析的R包自己的官方文档很清晰，见：https://github.com/bicciatolab/popsicleR ，包括：

Quality controls
Filtering
Doublet detection
Data normalization
Scaling and regression
Cell clustering
Cell annotation

现在让我们跟着其官方文档熟悉一下这个popsicleR用法吧：

示例数据

其使用了Oncogene . 2021 Jan的文章《Single-cell RNA sequencing reveals heterogeneous tumor and immune cell populations in early-stage lung adenocarcinomas harboring EGFR mutations》里面的一个样品，是 The sample N1 is composed of 14,972 single cells that were sequenced on the Illumina Hiseq X platform.

这个文章本身样品数量不少，在 https://ngdc.cncb.ac.cn/gsa/browse/CRA001963 是公开可以获取的。如下所示：

文章样品队列信息

可以看到，它的第一层次降维聚类分群比较传统，而且也给出来了各个单细胞亚群的标记基因，如下所示：

蛮有意思的，它的B细胞的基因在部分髓系细胞和T细胞里面都有混杂。不过我们现在是要介绍这个popsicleR用法，并不是想探究这个肺癌单细胞文章。查看示例数据，如下所示：

outputs/
├── [ 69K]  barcodes.tsv.gz
├── [298K]  features.tsv.gz
└── [ 72M]  matrix.mtx.gz

首先使用scRNAstat粗浅的看看这个数据集

我很久以前在《生信技能树》公众号的一个教程：这也能画？，提到了一个很无聊的R包，名字是：scRNAstat ，它可以4行代码进行单细胞转录组的降维聚类分群，其实完全没有技术含量，就是把 Seurat 流程的一些步骤包装成为了4个函数：

basic_qc (查看数据质量)
basic_filter （进行一定程度的过滤）
basic_workflow （降维聚类分群）
basic_markers（检查各个亚群的标记基因）

代码如下所示；

library(scRNAstat)
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
sce = CreateSeuratObject( Read10X('outputs/'))
sce
x='test'
dir.create(x)
sce = basic_qc(sce=sce,org='human',
               dir = x)
sce = basic_filter(sce)
sce
sce = basic_workflow(sce,dir = x)
markers_figures <- basic_markers(sce,
                                 org='human',
                                 group='seurat_clusters',
                                 dir = x) 
save(sce,file = 'sce_by_scRNAstat.Rdata')

我们简单的看了看其基因和基础分群：

如果你的背景知识足够强大，已经是可以比较好的给出了第一层次降维聚类分群的生物学名字啦。

4,6,8,11,15是上皮细胞，其中17是比较特殊的上皮细胞，大家会叫它干性上皮细胞。
5，7,18是T细胞，14是mast细胞，12是t和b两种淋巴细胞的混合，可能是双细胞或者处于细胞周期，需要具体看看。
10是成纤维细胞，16是内皮细胞，都是很出名的肿瘤微环境里面的基质细胞。
其余的细胞都是髓系，巨噬细胞或者单核细胞，其中9是比较清晰的树突细胞。

需要背诵如下所示各个细胞亚群高表达量基因的列表：

# T Cells (CD3D, CD3E, CD8A), 
# B cells (CD19, CD79A, MS4A1 [CD20]), 
# Plasma cells (IGHG1, MZB1, SDC1, CD79A), 
# Monocytes and macrophages (CD68, CD163, CD14),
# NK Cells (FGFBP2, FCG3RA, CX3CR1),  
# Photoreceptor cells (RCVRN), 
# Fibroblasts (FGF7, MME), 
# Endothelial cells (PECAM1, VWF). 
# epi or tumor (EPCAM, KRT19, PROM1, ALDH1A1, CD24).
#   immune (CD45+,PTPRC), epithelial/cancer (EpCAM+,EPCAM), 
# stromal (CD10+,MME,fibo or CD31+,PECAM1,endo)

接下来看看popsicleR用法

前面我们提到过，它官方文档很清晰，见：https://github.com/bicciatolab/popsicleR ，包括：

Quality controls
Filtering
Doublet detection
Data normalization
Scaling and regression
Cell clustering
Cell annotation

代码也是超级简单，同样的需要读入如下所示示例数据，一个文件夹名字是outputs，它里面需要有3个文件：

outputs/
├── [ 69K]  barcodes.tsv.gz
├── [298K]  features.tsv.gz
└── [ 72M]  matrix.mtx.gz

全部的代码是：

library(popsicleR)
# Define a list of cell-type markers
populations.markers <-  c("EPCAM","VIM","COL1A1","PECAM1","PTPRC","CD3D","CD14","SFTPC")

# Create a Seurat object from the raw data and visualize QC metrics
sample.umi <- PrePlots(sample = 'N1',
                       input_data = 'outputs/',
                       genelist = populations.markers,
                       percentage = 0.1, gene_filter = 200, cellranger = TRUE, organism = "human",
                       out_folder = 'test_popsicleR')
# Doublet detection
sample.umi <- CalculateDoublets(UMI = sample.umi, method = "scrublet", dbs_thr ='none', dbs_remove = FALSE, out_folder = main.dir)
# Doublet removal
sample.umi <- CalculateDoublets(UMI = sample.umi, method = "scrublet", dbs_thr = 0.22, dbs_remove = TRUE, out_folder = main.dir)
# Normalization
sample.umi <- Normalize(UMI = sample.umi, variable_genes = 2000, out_folder = main.dir)
# Regression
sample.umi <- ApplyRegression(UMI = sample.umi, organism = "human", variables = "none", explore_PC = FALSE, out_folder = main.dir)
# Regression
sample.umi <- ApplyRegression(UMI = sample.umi, organism = "human", variables = "none", explore_PC = TRUE, out_folder = main.dir)
# Cell type annotation
sample.umi <- MakeAnnotation(UMI = sample.umi, organism = "human", marker.list = "none", cluster_res= 0.8, out_folder = main.dir)

这里面的 scrublet包的CalculateDoublets 函数需要去检查数据集里面的双细胞情况，并且去除，所以非常耗费数据，因为单个样品拿到了14,972 个细胞确实数量有点多，不过其实我们的降维聚类分群也可以看到如果有双细胞也可以在命名之后再去除。参考；单个样品测序了近2万个单细胞怎么办

另外，这个popsicleR确实也值得尝试，毕竟产生了二十多张各种各样的图表，而且它代码有接近2000行，读一下源代码，学一学也挺好的。