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

百味科研芝士

发布于 2019-06-05 16:17:05

1K0

发布于 2019-06-05 16:17:05

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之外，其余的都找到了相对应的细胞类型（基于检验的，不是生物学的）。这当然可以为我们定义细胞类型做一个统计上的参考。

准备数据

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

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

Counts Matrix Number of gene tags per gene per cell.

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.

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

输入数据

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
....

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