本教程已更新,更新时间:
2019/12/27
一文介绍单细胞测序生物信息分析完整流程,这可能是最新也是最全的流程
cellranger mkfastq
可用于将单细胞测序获得的 BCL 文件拆分为可以识别的 fastq 测序数据
cellranger makefastq --run=[ ] --samplesheet=[sample.csv] --jobmode=local --localcores=20 --localmem=80
-–run :是下机数据 BCL 所在的路径;
-–samplesheet :样品信息列表--共三列(lane id ,sample name ,index name)
注意要控制好核心数和内存数
运行产出结果存在于 out
目录中
cellranger count
是 cellranger
最主要也是最重要的功能:完成细胞和基因的定量,也就是产生了我们用来做各种分析的基因表达矩阵。
cellranger count -–id=sample345 -–transcriptome=/opt/refdata-cellranger-GRCh38-1.2.0/GRCh38 -–fastqs=/home/jdoe/runs/HAWT7ADXX/outs/fastq_path -–indices=SI-3A-A1 –-cells=1000
id :产生的结果都在这个文件中,可以取几号样品(如 sample345);
fastqs :由 cellranger mkfastq 产生的 fastqs 文件夹所在的路径;fastqs :由 cellranger mkfastq 产生的 fastqs 文件夹所在的路径;
indices:sample index:SI-3A-A1;
transcriptome:参考转录组文件路径;
cells:预期回复的细胞数;
cellranger count
计算的结果只能作为初步观测的结果,如果需要进一步分析聚类细胞,还需要进行下游分析,这里使用官方推荐 R 包(Seurat 3.1)
流程参考官方(外周血分析标准流程)
install.packages('Seurat')library(Seurat)
library(dplyr)library(Seurat)# Load the PBMC datasetpbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")# Initialize the Seurat object with the raw (non-normalized data).pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)pbmc
## An object of class Seurat ## 13714 features across 2700 samples within 1 assay ## Active assay: RNA (13714 features)
这里读取的是单细胞 count 结果中的矩阵目录;
在对象生成的过程中,做了初步的过滤;
留下所有在>=3 个细胞中表达的基因min.cells
= 3;
为了除去一些质量差的细胞,留下所有检测到>=200 个基因的细胞min.genes
= 200。
# The [[ operator can add columns to object metadata. This is a great place to stash QC statspbmc[["percent.mt"]] <- PercentageFeatureSet(object = pbmc, pattern = "^MT-")
这一步 mit-开头的为线粒体基因,这里将其进行标记并统计其分布频率
# Visualize QC metrics as a violin plotVlnPlot(object = pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
对 pbmc 对象做小提琴图,分别为基因数,细胞数和线粒体占比
pbmc <- subset(x = pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
接下来,根据图片中基因数和线粒体数,分别设置过滤参数,这里基因数 200-2500,线粒体百分比为小于 5%
pbmc <- NormalizeData(object = pbmc, normalization.method = "LogNormalize", scale.factor = 10000)pbmc <- NormalizeData(object = pbmc)
pbmc <- FindVariableFeatures(object = pbmc, selection.method = "vst", nfeatures = 2000)# Identify the 10 most highly variable genestop10 <- head(x = VariableFeatures(object = pbmc), 10)# plot variable features with and without labelsplot1 <- VariableFeaturePlot(object = pbmc)plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)CombinePlots(plots = list(plot1, plot2))
all.genes <- rownames(x = pbmc)pbmc <- ScaleData(object = pbmc, features = all.genes)
这里设置对所有的基因都做了scale
,但是需要知道的是,其实后续的分析都是基于高变基因的,因此,使用默认参数就可以了,而且提升效率。
pbmc <- ScaleData(object = pbmc)
pbmc <- RunPCA(object = pbmc, features = VariableFeatures(object = pbmc))
这里有多种方法展示 pca 结果,本文采用最简单的方法
DimPlot(object = pbmc, reduction = "pca")
pbmc <- JackStraw(object = pbmc, num.replicate = 100)pbmc <- ScoreJackStraw(object = pbmc, dims = 1:20)JackStrawPlot(object = pbmc, dims = 1:15)
虚线以上的为可用维度,你也可以调整 dims 参数,画出所有 pca 查看
另外一种鉴定手段是绘制所有PC的分布点图
ElbowPlot(pbmc)
大多数软件都是通过拾取拐点处的pc作为选定数目
pbmc <- FindNeighbors(object = pbmc, dims = 1:10)pbmc <- FindClusters(object = pbmc, resolution = 0.5)
这里的 dims 为上一步计算所用的维度数,而 resolution 参数控制聚类的数目,针对3K的细胞数目,最好的范围是0.4-1.2
head(Idents(pbmc), 5)
## AAACATACAACCAC AAACATTGAGCTAC AAACATTGATCAGC AAACCGTGCTTCCG AAACCGTGTATGCG ## 1 3 1 2 6 ## Levels: 0 1 2 3 4 5 6 7 8
这里注意,这一步聚类有两种聚类方法(umap/tSNE),两种方法都可以使用,但不要混用,这样,后面的结算结果会将先前的聚类覆盖掉,只能保留一个
本文采用基于umap的聚类方法
pbmc <- RunUMAP(object = pbmc, dims = 1:10)DimPlot(object = pbmc, reduction = "umap")
完成聚类后,一定要记住保存数据,不然重新计算可要头疼了
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")
cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 1, min.pct = 0.25)head(x = cluster1.markers, n = 5)
这样是寻找单个聚类中的显著基因
cluster5.markers <- FindMarkers(object = pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)head(x = cluster5.markers, n = 5)
这样寻找所有聚类中显著基因,计算速度很慢,需要等待
另外,我们有多种方法统计基因的显著性
FeaturePlot(object = pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)DoHeatmap(object = pbmc, features = top10$gene) + NoLegend()
剩下的便是寻找基因 marker 并对细胞类型进行注释
你可能想要知道如何对多样本进行整合分析,请参考教程:整合刺激性和对照性PBMC数据集,以学习细胞类型特异性反应
你可能想了解如何灵活操作Seurat的S4对象,以便轻松的提取表达矩阵;亦或者想要做一些不一样的可视化效果,例如采用平均表达量做热图展示等,请参考Seurat3.1的灵活操作指南
你可能想对单细胞数据做进一步分析,例如功能分析,请参考10X单细胞数据针对细胞及其亚型的基因集功能分析和转录调控分析
众所周知,细胞类型的注释是最困难的一步,除非你有很强的对细胞基因的敏感度,不然很难识别细胞类型。
通常识别细胞类型的方法主要有三种
根据Marker基因,采用CellMarker或者panglaoDB数据库,进行细胞注释,可以采用超几何分布算法来进行精确性验证,方法参考ClusterProfiler:真的不只是富集分析
从文献中获取已经验证的Marker
这里简单介绍下SingleR
SingleR
:一个全自动细胞注释的 R 包,用法很简单
BiocManager::install("SingleR")browseVignettes("SingleR")
从头预测的方法请参考官方教程
因为我们刚刚从Seurat过来的,所以我们应该很想知道Seurat cluster 的细胞注释结果,因此,对Seurat的结果进行注释
我们这里采用两个人类的参考集去做细胞注释
library(Seurat)library(SingleR)library(dplyr)library(tibble)hpca.se <- HumanPrimaryCellAtlasData()bpe.se <- BlueprintEncodeData()
读入Seurat
对象转换为SingleCell
支持的对象
seurat.obj <- readRDS("../output/pbmc_tutorial.rds")seurat.obj@meta.data$cell.type <- Idents(seurat.obj)test <- as.SingleCellExperiment(seurat.obj)
采用两个参考集一起进行注释
Anno <- SingleR(test = test, ref = list(HP = hpca.se , BP = bpe.se), labels = list(hpca.se$label.main , bpe.se$label.main), method = "cluster", cluster = test$cell.type)
提取需要的细胞分类信息
Anno$cluster <- rownames(Anno)fin <- Anno %>% dplyr::tbl_df() %>% dplyr::select(cluster,labels)
你也可以将细胞注释信息重新添加到Seurat
对象中去
new.cluster.ids <- fin$labelsnames(new.cluster.ids) <- levels(seurat.obj)seurat.obj <- RenameIdents(seurat.obj, new.cluster.ids)
伪时间分析建议采用monocle3
软件
##安装依赖BiocManager::install(c('BiocGenerics', 'DelayedArray', 'DelayedMatrixStats', 'limma', 'S4Vectors', 'SingleCellExperiment', 'SummarizedExperiment', 'batchelor'))##安装monocle3devtools::install_github('cole-trapnell-lab/leidenbase')devtools::install_github('cole-trapnell-lab/monocle3')library(monocle3)
library(Seurat)library(monocle3)endo<-readRDS("../output/pbmc_tutorial.rds")data <- endo@assays$RNA@countspd <- endo@meta.data
不要直接把meta.data
放入pd
中去,太多没用的信息了,先清理下
new_pd<-select(pd,nCount_RNA,nFeature_RNA,percent.mt)new_pd$Cell.type<-Idents(endo)head(new_pd) nCount_RNA nFeature_RNA percent.mtAAACCTGGTATCAGTC-1_1 4835 1383 3.536711AACCGCGCATTCCTGC-1_1 4029 1230 2.655746AACGTTGAGCCCAATT-1_1 2576 918 3.377329AAGACCTGTGGTTTCA-1_1 2841 1044 5.455825AAGTCTGCATGAAGTA-1_1 3731 1213 3.886358ACAGCCGCATCGGACC-1_1 5915 1918 3.043111 Cell.typeAAACCTGGTATCAGTC-1_1 DC2AACCGCGCATTCCTGC-1_1 MDSCAACGTTGAGCCCAATT-1_1 MDSCAAGACCTGTGGTTTCA-1_1 Dendritic.cells.activatedAAGTCTGCATGAAGTA-1_1 MDSCACAGCCGCATCGGACC-1_1 Dendritic.cells.activatedfData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
生成monocle
的cds
对象
cds <- new_cell_data_set(data,cell_metadata = new_pd,gene_metadata = fData)
运行下游标准流程,无论如何都得跑,是必须的默认流程
cds <- preprocess_cds(cds, num_dim = 30)#umapcds <- reduce_dimension(cds,umap.n_neighbors = 20L)#color by seurat clusterplot_cells(cds,label_groups_by_cluster=FALSE,color_cells_by = "Cell.type")
#clustercds <- cluster_cells(cds,resolution = 0.5)#color by monocle clusterplot_cells(cds, color_cells_by = "partition",label_groups_by_cluster=FALSE)
cds <- learn_graph(cds)plot_cells(cds,color_cells_by = "Cell.type",label_groups_by_cluster=FALSE,label_leaves=TRUE,label_branch_points=TRUE)
##一个有用的寻找起源节点的函数get_earliest_principal_node <- function(cds, time_bin="Dendritic.cells.activated"){ cell_ids <- which(colData(cds)[, "Cell.type"] == time_bin) closest_vertex <- cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex closest_vertex <- as.matrix(closest_vertex[colnames(cds), ]) root_pr_nodes <- igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names (which.max(table(closest_vertex[cell_ids,]))))] root_pr_nodes}#order cellcds <- order_cells(cds, root_pr_nodes=get_earliest_principal_node(cds))#pseudotime analysisplot_cells(cds,color_cells_by = "pseudotime",label_cell_groups=FALSE,label_leaves=FALSE,label_branch_points=FALSE,graph_label_size=1.5)
联系客服