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

2019-05-31 本文已影响74人周运来就是我

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

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

且慢，用非监督的方式聚出来的类（一般叫做 $Cluster_i,i$ 为分的类数），既然是细胞异质性的体现，那么它到底是什么细胞呢，可以叫做什么，是CD4细胞呢还是microglia细胞？因为 $cluster_{i}$ 是没有意义的，只有起了像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 data
toy_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 data
toy_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 expression
de_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”基因在参考数据库中的分布。

# Plot
make_ranking_violin_plot(de_table.test=de_table.toy_query, de_table.ref=de_table.toy_ref)

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

# And get group labels
make_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的输出结果，注意为了和官网数据格式保持一致，我把 $otus$ 文件夹命名为 $10X_pbmc4k$ 了，这个前后一致就可以了。

+--- 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_k7
class: 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-1
colData names(2): cell_sample group

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

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 $Watkins et al. (2009)$ . 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 $‘haemosphere’$ 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 NA
probe_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 probes
genes_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细胞的定义工作就完成了。官网也提供了小鼠细胞类型识别的教程。这里就不再复述啦。

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

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