celaref ||单细胞细胞类型定义工具

Sarah Williams (2019). celaref: Single-cell RNAseq cell cluster labelling by reference. R package version 1.0.1.

当我们拿到单细胞转录组数据的时候,不管结果怎么样,翻来覆去发现其实就是一个基因和细胞的矩阵而已!一般基于这个矩阵,会先过滤细胞均一化,降维,聚类(最重要的一步,因为它可以发现细胞的异质性),聚类结果的差异分析(声称找到了每个亚群的marker基因),最后做一下数据的可视化。

且慢,用非监督的方式聚出来的类(一般叫做

为分的类数),既然是细胞异质性的体现,那么它到底是什么细胞呢,可以叫做什么,是CD4细胞呢还是microglia细胞?因为

是没有意义的,只有起了像CD4这样的名字,这类细胞才有故事可以讲。

在定义细胞类型之前,需要确定就哪种聚类结果来做,是图聚类的结果还是k-means某一类的结果。如何来确定?看看分群的tsne图是一个不错的参考。

那么呢,之前你们家大师推荐过几个数据库single cell marker 基因数据库是基于某亚群细胞的marker基因,根据已经知道marker基因的细胞类型,这些数据库其实就是一个细胞类型和marker基因大表。根据marker基因列表与分析出来的差异基因列表以及基因丰度的关系(往往是做一个热图)来推断某个cluster属于哪一种细胞类型。

另外还有一种定义细胞类型的方法是通过每个cluster与已知细胞类型的表达谱的相似性来定义细胞类型,已有的方法为SingleR、celaref。均是R包,其中celaref是2019年的新品,也是本文主要介绍的一种方法。

celaref 简介

走遍示例流程

安装celaref包的时候,就把示例数据也下载了,首先载入数据,然后观察一下示例数据是怎么组织的,以便我们以后组织自己的数据。在接触一个新包的时候,最好的学习方法就是用示例数据走一遍然后看看有哪些功能(函数以及函数的参数可以用),多用help或者?来查看细节。一个R包封装的越好也就要求我们花时间去了解他的细节。

library(celaref)# Paths to data files.counts_filepath.query <- system.file("extdata", "sim_query_counts.tab", package = "celaref") cell_info_filepath.query <- system.file("extdata", "sim_query_cell_info.tab", package = "celaref") counts_filepath.ref <- system.file("extdata", "sim_ref_counts.tab", package = "celaref") cell_info_filepath.ref <- system.file("extdata", "sim_ref_cell_info.tab", package = "celaref")# Load datatoy_ref_se <- load_se_from_files(counts_file=counts_filepath.ref, cell_info_file=cell_info_filepath.ref) head(toy_ref_se) toy_query_se <- load_se_from_files(counts_file=counts_filepath.query, cell_info_file=cell_info_filepath.query)

可见我载入了两个数据:一个是toy_ref_se,reference;一个是toy_query_se,query.都是SummarizedExperiment对象:

对数据进行过滤,help一下这个函数,看看都有哪些过滤条件。

# Filter datatoy_ref_se <- trim_small_groups_and_low_expression_genes(toy_ref_se) toy_query_se <- trim_small_groups_and_low_expression_genes(toy_query_se)

对参考和需要鉴定的表达谱分别做差异分析:

# Setup within-experiment differential expressionde_table.toy_ref <- contrast_each_group_to_the_rest(toy_ref_se, dataset_name="ref") de_table.toy_query <- contrast_each_group_to_the_rest(toy_query_se, dataset_name="query")

分别得到差异基因结果:

绘制一组小提琴图,显示每个查询组的“top”基因在参考数据库中的分布。

# Plotmake_ranking_violin_plot(de_table.test=de_table.toy_query, de_table.ref=de_table.toy_ref)

根据参考数据集的相似性,在查询数据集中构造一些合理的标签或组/集群。

# And get group labelsmake_ref_similarity_names(de_table.toy_query, de_table.toy_ref)

可以看到除了Group2 没有相应的similarity 类似的type之外,其余的都找到了相对应的细胞类型(基于检验的,不是生物学的)。这当然可以为我们定义细胞类型做一个统计上的参考。

准备数据

那么,如何应用我们自己的数据来做呢?官网也给出了数据准备的方法:

在准备数据时需要以下数据:

A.

Counts Matrix Number of gene tags per gene per cell.

B.

Cell information A sample-sheet table of cell-level information. Only two fields are essential: cell_sample: A unique cell identifier group: The cluster/group to which the cell has been assigned.

C.

Gene information Optional. A table of extra gene-level information. ID: A unique gene identifier

输入数据

01

tab键分割的ounts matrix

gene

Cell1

cell2

cell3

cell4

cell954

GeneA

0

1

0

1

0

GeneB

0

3

0

2

2

GeneC

1

40

1

0

0

....

02

tab键分割的细胞分群信息

CellId

Sample

Cluster

cell1

Control

cluster1

cell2

Control

cluster7

cell954

KO

cluster8

From 10X pipeline output

dataset_se <- load_dataset_10Xdata('~/path/to/data/10X_mydata', dataset_genome = "GRCh38", clustering_set = "kmeans_7_clusters")

拿一个10X的例子试试吧

Prepare 10X query dataset

先要把数据下载好,发现这应用的对象还是cellranger V2 pipeline的输出结果,注意为了和官网数据格式保持一致,我把

文件夹命名为

了,这个前后一致就可以了。

+--- 10X_pbmc4k|   +--- analysis
|   |   +--- clustering
|   |   |   +--- graphclust|   |   |   |   +--- clusters.csv|   |   |   +--- kmeans_10_clusters
|   |   |   |   +--- clusters.csv
|   |   |   +--- kmeans_2_clusters
···|   |   +--- diffexp|   |   |   +--- graphclust
|   |   |   |   +--- differential_expression.csv
|   |   |   +--- kmeans_10_clusters
···|   |   +--- pca|   |   |   +--- 10_components
|   |   |   |   +--- components.csv
|   |   |   |   +--- dispersion.csv
|   |   |   |   +--- genes_selected.csv
|   |   |   |   +--- projection.csv
|   |   |   |   +--- variance.csv
|   |   +--- tsne
|   |   |   +--- 2_components|   |   |   |   +--- projection.csv|   +--- filtered_gene_bc_matrices
|   |   +--- GRCh38
|   |   |   +--- barcodes.tsv|   |   |   +--- genes.tsv
|   |   |   +--- matrix.mtx
datasets_dir <- "D:\\Users\\Administrator\\Desktop\\RStudio\\single_cell\\celaref/celaref_extra_vignette_data/datasets"dir(datasets_dir)

dataset_se.10X_pbmc4k_k7 <- load_dataset_10Xdata(
  dataset_path   = file.path(datasets_dir,'10X_pbmc4k'), 
  dataset_genome = "GRCh38", 
  clustering_set = "kmeans_7_clusters", 
  id_to_use      = "GeneSymbol")
dataset_se.10X_pbmc4k_k7 <- trim_small_groups_and_low_expression_genes(dataset_se.10X_pbmc4k_k7)

可见我选择的是kmeans_7_clusters 的聚类结果来进行细胞定义。处理过之后:

> dataset_se.10X_pbmc4k_k7class: SummarizedExperiment dim: 15407 4340 metadata(0): assays(1): ''rownames(15407): A1BG A1BG-AS1 ... ZZEF1 ZZZ3 rowData names(4): ensembl_ID GeneSymbol ID total_count colnames(4340): AAACCTGAGAAGGCCT-1 AAACCTGAGACAGACC-1 ... TTTGTCAGTTAAGACA-1 TTTGTCATCCCAAGAT-1colData names(2): cell_sample group

下一步是对数据的转化和处理,也是one-to-others的DE分析,非常的慢和消耗内存,有可能跑断的。关键是这个函数做了是什么?其实就是把我们的矩阵文件整理成上面演示的

格式。

de_table.10X_pbmc4k_k7 <- contrast_each_group_to_the_rest(dataset_se.10X_pbmc4k_k7, dataset_name="10X_pbmc4k_k7", num_cores=7) Sorry, multithreading not supported on windows. Use linux, or set num_threads = 1 to suppress this message.Running single threaded # 这个sorry 可以忽略 ,在Linux下才能指定运行进程数。`fData` has no primerid. I'll make something up. `cData` has no wellKey. I'll make something up. Assuming data assay in position 1, with name et is log-transformed.

Done! Combining coefficients and standard errors Calculating log-fold changes Calculating likelihood ratio tests Refitting on reduced model...

Done!`fData` has no primerid. I'll make something up. `cData` has no wellKey. I'll make something up. Assuming data assay in position 1, with name et is log-transformed. Completed [=======>----------------------------------------------------------------------------------------------------------------] 7% with 0 failures

windows下居然跑了一上午

准备 reference数据集

来自响应的数据库,因为我们的是血细胞,选择的是haemosphere .

A suitable PBMC reference (a ‘HaemAtlas’) has been published by

. They purified populations of PBMC cell types and measured gene expression via microarray. The data used here was downloaded in a normalised table from the

website (Graaf et al. 2016).

下载页面

下载完之后由于是微阵列数据需要做一些特殊处理。

  • Logged, normalised expression values. Any low expression or poor quality measurements should have already been removed.
  • Sample information (see contrast_each_group_to_the_rest_for_norm_ma_with_limma for details)

library(readr) this_dataset_dir <- file.path(datasets_dir, 'haemosphere_datasets','watkins') norm_expression_file <- file.path(this_dataset_dir, "watkins_expression.txt") samples_file <- file.path(this_dataset_dir, "watkins_samples.txt") norm_expression_table.full <- read.table(norm_expression_file, sep="\t", header=TRUE, quote="", comment.char="", row.names=1, check.names=FALSE) samples_table <- read_tsv(samples_file, col_types = cols()) samples_table$description <- make.names( samples_table$description) # Avoid group or extra_factor names starting with numbers, for microarrays> samples_table$description [1] "X49.years.adult" "X49.years.adult" "X49.years.adult" "X49.years.adult" "X49.years.adult" "X49.years.adult" "X49.years.adult" [8] "X49.years.adult" "X49.years.adult" "X49.years.adult" "X49.years.adult" "X49.years.adult" "X63.years.adult" "X63.years.adult"[15] "X63.years.adult" "X63.years.adult" "X63.years.adult" "X63.years.adult" "X53.years.adult" "X53.years.adult" "X53.years.adult"[22] "X53.years.adult" "X53.years.adult" "X53.years.adult" "X60.years.adult" "X60.years.adult" "X60.years.adult" "X60.years.adult"[29] "X60.years.adult" "X60.years.adult" "X48.years.adult" "X48.years.adult" "X48.years.adult" "X48.years.adult" "X48.years.adult"[36] "X48.years.adult" "X57.years.adult" "X57.years.adult" "X57.years.adult" "X57.years.adult" "X57.years.adult" "X57.years.adult"[43] "NA." "NA." "NA." "NA." "NA." "NA." "NA." [50] "NA."

从样本表中可以看出,这个数据集包含了其他组织,但是作为PBMC的参考,我们只考虑外周血样本。与其他数据加载函数一样,要从分析中删除一个样本(或单元格),从样本表中删除它就ok了。

amples_table <- samples_table[samples_table$tissue == "Peripheral Blood",] > samples_table# A tibble: 42 x 7 sampleId celltype cell_lineage surface_markers tissue description notes <chr> <chr> <chr> <lgl> <chr> <chr> <lgl> 1 1674120023_A B lymphocyte B Cell Lineage NA Peripheral Blood X49.years.adult NA 2 1674120023_B granulocyte Neutrophil Lineage NA Peripheral Blood X49.years.adult NA 3 1674120023_C natural killer cell NK Cell Lineage NA Peripheral Blood X49.years.adult NA 4 1674120023_D Th lymphocyte T Cell Lineage NA Peripheral Blood X49.years.adult NA 5 1674120023_E Tc lymphocyte T Cell Lineage NA Peripheral Blood X49.years.adult NA 6 1674120023_F monocyte Macrophage Lineage NA Peripheral Blood X49.years.adult NA 7 1674120053_A B lymphocyte B Cell Lineage NA Peripheral Blood X49.years.adult NA 8 1674120053_B monocyte Macrophage Lineage NA Peripheral Blood X49.years.adult NA 9 1674120053_C granulocyte Neutrophil Lineage NA Peripheral Blood X49.years.adult NA 10 1674120053_D Tc lymphocyte T Cell Lineage NA Peripheral Blood X49.years.adult NA # ... with 32 more rows

通常情况下,最难的部分是格式化输入。应该使用与查询数据集相同的id,将微阵列表达式值作为规范化的、日志转换的数据。任何探头或样品级的过滤也应事先进行。在这种情况下,从网站获得的数据是正常的,但仍然需要匹配id。

library("tidyverse") library("illuminaHumanv2.db") probes_with_gene_symbol_and_with_data <- intersect(keys(illuminaHumanv2SYMBOL),rownames(norm_expression_table.full))# Get mappings - non NAprobe_to_symbol <- select(illuminaHumanv2.db, keys=rownames(norm_expression_table.full), columns=c("SYMBOL"), keytype="PROBEID") probe_to_symbol <- unique(probe_to_symbol[! is.na(probe_to_symbol$SYMBOL),])# no multimapping probesgenes_per_probe <- table(probe_to_symbol$PROBEID) # How many genes a probe is annotated against?multimap_probes <- names(genes_per_probe)[genes_per_probe > 1] probe_to_symbol <- probe_to_symbol[!probe_to_symbol$PROBEID %in% multimap_probes, ] convert_expression_table_ids<- function(expression_table, the_probes_table, old_id_name, new_id_name){ the_probes_table <- the_probes_table[,c(old_id_name, new_id_name)] colnames(the_probes_table) <- c("old_id", "new_id") # Before DE, just pick the top expresed probe to represent the gene # Not perfect, but this is a ranking-based analysis. # hybridisation issues aside, would expect higher epressed probes to be more relevant to Single cell data anyway. probe_expression_levels <- rowSums(expression_table) the_probes_table$avgexpr <- probe_expression_levels[as.character(the_probes_table$old_id)] the_genes_table <- the_probes_table %>% group_by(new_id) %>% top_n(1, avgexpr) expression_table <- expression_table[the_genes_table$old_id,] rownames(expression_table) <- the_genes_table$new_id return(expression_table) }# Just the most highly expressed probe foreach gene.norm_expression_table.genes <- convert_expression_table_ids(norm_expression_table.full, probe_to_symbol, old_id_name="PROBEID", new_id_name="SYMBOL")# Go...de_table.Watkins2009PBMCs <- contrast_each_group_to_the_rest_for_norm_ma_with_limma( norm_expression_table = norm_expression_table.genes, sample_sheet_table = samples_table, dataset_name = "Watkins2009PBMCs", extra_factor_name = 'description', sample_name = "sampleId", group_name = 'celltype') > de_table.Watkins2009PBMCs# A tibble: 113,376 x 21 ID logFC AveExpr t P.Value adj.P.Val B ci CI_upper CI_lower ci_inner ci_outer log2FC fdr group sig sig_up <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <lgl> <lgl> 1 JCHA~ 6.86 8.65 27.5 7.13e-26 2.70e-23 49.0 0.0952 6.96 6.77 6.77 6.96 6.86 2.70e-23 B ly~ TRUE TRUE 2 FCRLA 6.41 8.77 23.7 1.13e-23 3.29e-21 44.0 0.103 6.51 6.30 6.30 6.51 6.41 3.29e-21 B ly~ TRUE TRUE 3 VPRE~ 6.18 8.66 40.2 1.20e-31 8.70e-29 62.0 0.0587 6.23 6.12 6.12 6.23 6.18 8.70e-29 B ly~ TRUE TRUE 4 TCL1A 6.09 8.52 42.8 1.33e-32 1.14e-29 64.1 0.0544 6.14 6.04 6.04 6.14 6.09 1.14e-29 B ly~ TRUE TRUE 5 CD79A 6.09 9.08 22.7 5.17e-23 1.40e-20 42.4 0.102 6.20 5.99 5.99 6.20 6.09 1.40e-20 B ly~ TRUE TRUE 6 CD19 5.90 8.44 38.4 6.09e-31 4.11e-28 60.5 0.0587 5.96 5.85 5.85 5.96 5.90 4.11e-28 B ly~ TRUE TRUE 7 HLA-~ 5.61 8.35 47.9 2.30e-34 2.56e-31 68.0 0.0447 5.66 5.57 5.57 5.66 5.61 2.56e-31 B ly~ TRUE TRUE 8 OSBP~ 5.48 7.83 53.2 5.60e-36 8.82e-33 71.4 0.0394 5.52 5.45 5.45 5.52 5.48 8.82e-33 B ly~ TRUE TRUE 9 CNTN~ 5.38 7.95 49.3 8.12e-35 9.60e-32 68.9 0.0416 5.42 5.34 5.34 5.42 5.38 9.60e-32 B ly~ TRUE TRUE 10 COBL~ 5.26 7.70 56.4 6.77e-37 1.28e-33 73.3 0.0356 5.30 5.23 5.23 5.30 5.26 1.28e-33 B ly~ TRUE TRUE # ... with 113,366 more rows, and 4 more variables: gene_count <int>, rank <int>, rescaled_rank <dbl>, dataset <chr>

细胞类型定义(Compare 10X query PBMCs to to reference)

数据准备好就和之前的示例是一样的了。

make_ranking_violin_plot(de_table.test=de_table.10X_pbmc4k_k7, de_table.ref=de_table.Watkins2009PBMCs)

Hmm, there’s a few clusters where different the top genes are bunched near the top for a couple of different reference cell types.

Logging the plot will be more informative at the top end for this dataset.

make_ranking_violin_plot(de_table.test=de_table.10X_pbmc4k_k7, de_table.ref=de_table.Watkins2009PBMCs, log10trans = TRUE)

可以打标签啦。

label_table.pbmc4k_k7_vs_Watkins2009PBMCs <- make_ref_similarity_names(de_table.10X_pbmc4k_k7, de_table.Watkins2009PBMCs)label_table.pbmc4k_k7_vs_Watkins2009PBMCs

可以看出,七个群有六个找到了与之相似的群名称,并给出了pval.

make_ranking_violin_plot(de_table.test=de_table.Watkins2009PBMCs, de_table.ref=de_table.10X_pbmc4k_k7, log10trans = TRUE)

行啦,人pbmc细胞的定义工作就完成了。官网也提供了小鼠细胞类型识别的教程。这里就不再复述啦。

完成了细胞类型定义,单细胞的故事才刚刚开始。所以,这不是结束,甚至不是结束的开始,而可能是开始的结束。

本文转载自简书

© 著作权归作者所有

本文分享自微信公众号 - 百味科研芝士(keyanzhishi)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-05-31

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

扫码关注云+社区

领取腾讯云代金券