查看: 2540|回复: 8

[R语言] 使用CHETAH包进行单细胞类型注释分析

[复制链接]

版主

Rank: 10Rank: 10Rank: 10

主题
25
注册时间
2016.5.11
在线时间
496 小时

活跃会员突出贡献


发表于 2020.6.6 14:51:02 | 显示全部楼层 |阅读模式
本帖最后由 bioinfomics 于 2020.6.6 14:50 编辑
本文首发于微信公众号“bioinfomics”:https://mp.weixin.qq.com/s/tOjE1BEcOhYdDX0zARWCnw
CHETAH包简介CHETAH(CHaracterization of cEll Types Aided by Hierarchical classification,通过层级分类辅助鉴定细胞类型)是用于单细胞RNA-seq测序(scRNA-seq)数据的细胞类型识别的R包。

CHETAH包通过以层级分类方式将输入数据与参考数据集相关联来分配细胞类型。如果输入数据不能完全分类为参考数据中的一种细胞类型,则可能被分配给中间类型。CHETAH可以与scRNA-seq参考数据集一起使用,也可以与常规RNA-seq或micro-array的参考数据集一起使用。因此,要运行CHETAH,我们只需要提供两个输入文件,它们都是SingleCellExperiment对象:
  • 输入数据(input data)
  • 参考数据集(reference dataset),带有已知的细胞类型

CHETAH包通过对参考数据集进行分层聚类(hierarchically clustering)来构建分类树,并以此分类树为指导,在分类树的每个节点中,输入细胞要么分配给右边的分支,要么分配给左边的分支。CHETAH会为每个分配结果计算一个置信度分数,当分配的置信度得分低于阈值(default = 0.1)时,则该细胞的分类将在该节点处停止。
这会产生两种类型的分类结果:
  • 最终类型(final types):细胞被分类为分类树的叶节点之一(即参考的细胞类型)。
  • 中间类型(intermediate types):将置信度得分低于某个节点中的阈值的细胞分配给分类树的中间节点。当一个细胞与某个节点的右分支和左分支中的细胞类型具有近似相同的相似性时,就会产生这种情况。




CHETAH包的安装
CHETAH will be a part of Bioconductor starting at release 2.9 (30th of April), and will be available by:
[AppleScript] 纯文本查看 复制代码
## Install BiocManager is neccesary
if (!require("BiocManager")) {
    install.packages("BiocManager")
}
BiocManager::install('CHETAH')

# Load the package
library(CHETAH)

The development version can be downloaded from the development version of Bioconductor (in R v3.6).
[AppleScript] 纯文本查看 复制代码
## Install BiocManager is neccesary
if (!require("BiocManager")) {
    install.packages("BiocManager")
}
BiocManager::install('CHETAH', version = "devel")

# Load the package
library(CHETAH)


At a glance
To run chetah on an input count matrix input_counts with t-SNE coordinates in input_tsne, and a reference count matrix ref_counts with celltypes vector ref_ct, run:
[AppleScript] 纯文本查看 复制代码
# Load the package
library(CHETAH)

## Make SingleCellExperiments
reference <- SingleCellExperiment(assays = list(counts = ref_counts),
                                  colData = DataFrame(celltypes = ref_ct))

input <- SingleCellExperiment(assays = list(counts = input_counts),
                              reducedDims = SimpleList(TSNE = input_tsne))

## Run CHETAH
input <- CHETAHclassifier(input = input, ref_cells = reference)

## Plot the classification
PlotCHETAH(input)

## Extract celltypes:
celltypes <- input$celltype_CHETAH


准备输入数据
Required data
如果输入数据存储为SingleCellExperiments对象,请继续执行下一步。否则,在开始运行CHETAH之前,我们需要准备以下数据:
  • 输入细胞的scRNA-seq表达计数矩阵:数据框或矩阵格式,列为细胞,行为基因
  • 标准化后的参考细胞的scRNA-seq表达计数矩阵
  • 参考细胞的细胞类型
  • (可选)用于可视化的输入数据的2D降维表示形式:如t-SNE,PCA。

作为有关如何准备数据的示例,我们将使用Tirosh等人的黑色素瘤数据作为输入数据,以Puram等人的头颈肿瘤数据作为参考数据集。
[AppleScript] 纯文本查看 复制代码
## To prepare the data from the package's internal data, run:
# 参考数据的细胞类型
celltypes_hn <- headneck_ref$celltypes
# 参考数据的表达矩阵
counts_hn <- assay(headneck_ref)
# 输入数据的表达矩阵
counts_melanoma <- assay(input_mel)
# 输入数据的降维信息
tsne_melanoma <- reducedDim(input_mel)

## The input data: a Matrix
class(counts_melanoma)
#> [1] "dgCMatrix"
#> attr(,"package")
#> [1] "Matrix"

counts_melanoma[1:5, 1:5]
#> 5 x 5 sparse Matrix of class "dgCMatrix"
#>              mel_cell1 mel_cell2 mel_cell3 mel_cell4 mel_cell5
#> ELMO2           .         .         .         4.5633    .     
#> PNMA1           .         4.3553    .         .         .     
#> MMP2            .         .         .         .         .     
#> TMEM216         .         .         .         .         5.5624
#> TRAF3IP2-AS1    2.1299    4.0542    2.4209    1.6531    1.3144

## The reduced dimensions of the input cells: 2 column matrix
tsne_melanoma[1:5, ]
#>               tSNE_1    tSNE_2
#> mel_cell1  4.5034553 13.596680
#> mel_cell2 -4.0025667 -7.075722
#> mel_cell3  0.4734054  9.277648
#> mel_cell4  3.2201815 11.445236
#> mel_cell5 -0.3354758  5.092415

all.equal(rownames(tsne_melanoma), colnames(counts_melanoma))
#> [1] TRUE

## The reference data: a Matrix
class(counts_hn)
#> [1] "matrix" "array"

counts_hn[1:5, 1:5]
#>              hn_cell1 hn_cell2 hn_cell3 hn_cell4 hn_cell5
#> ELMO2         0.00000        0  0.00000  1.55430   4.2926
#> PNMA1         0.00000        0  0.00000  4.55360   0.0000
#> MMP2          0.00000        0  7.02880  4.50910   6.3006
#> TMEM216       0.00000        0  0.00000  0.00000   0.0000
#> TRAF3IP2-AS1  0.14796        0  0.65352  0.28924   3.6365

## The cell types of the reference: a named character vector
str(celltypes_hn)
#>  Named chr [1:180] "Fibroblast" "Fibroblast" "Fibroblast" "Fibroblast" ...
#>  - attr(*, "names")= chr [1:180] "hn_cell1" "hn_cell2" "hn_cell3" "hn_cell4" ...

## The names of the cell types correspond with the colnames of the reference counts:
all.equal(names(celltypes_hn), colnames(counts_melanoma)) 
#> [1] "Lengths (180, 150) differ (string compare on first 150)"
#> [2] "150 string mismatches"


SingleCellExperiments
CHETAH包以SingleCellExperiment对象格式的数据作为输入,这是一种将不同类型的数据存储在一起的简便方法。
A SingleCellExperiment holds three things:
  • counts: assays(as a list of Matrices)
  • meta-data: colData(as DataFrames)
  • reduced dimensions (e.g. t-SNE, PCA): ReducedDims(as a SimpleList of 2-column data.frames or matrices)

CHETAH needs:
  • a reference SingleCellExperiment with:
1) an assay
2) a colData column with the corresponding cell types (default “celltypes”)
  • an input SingleCellExperiment with:
1) an assay
2) a reducedDim (e.g. t-SNE)

For the example data, we would make the two objects by running:
[AppleScript] 纯文本查看 复制代码
## For the reference we define a "counts" assay and "celltypes" metadata
# 构建参考数据SingleCellExperiment对象
headneck_ref <- SingleCellExperiment(assays = list(counts = counts_hn),
                                     colData = DataFrame(celltypes = celltypes_hn))

headneck_ref
class: SingleCellExperiment 
dim: 7943 180 
metadata(0):
assays(1): counts
rownames(7943): ELMO2 PNMA1 ... SLC39A6 CTSC
rowData names(0):
colnames(180): hn_cell1 hn_cell2 ... hn_cell179 hn_cell180
colData names(1): celltypes
reducedDimNames(0):
spikeNames(0):

## For the input we define a "counts" assay and "TSNE" reduced dimensions
# 构建输入数据SingleCellExperiment对象
input_mel <- SingleCellExperiment(assays = list(counts = counts_melanoma),
                                  reducedDims = SimpleList(TSNE = tsne_melanoma))

input_mel
class: SingleCellExperiment 
dim: 7943 150 
metadata(0):
assays(1): counts
rownames(7943): ELMO2 PNMA1 ... SLC39A6 CTSC
rowData names(0):
colnames(150): mel_cell1 mel_cell2 ... mel_cell149 mel_cell150
colData names(0):
reducedDimNames(1): TSNE
spikeNames(0):


运行CHETAH
Now that the data is prepared, running chetah is easy:
[AppleScript] 纯文本查看 复制代码
# 使用CHETAHclassifier函数运行CHETAH进行细胞类型分类注释
input_mel <- CHETAHclassifier(input = input_mel,
                              ref_cells = headneck_ref)
#> Preparing data....
#> Running analysis...

input_mel
class: SingleCellExperiment 
dim: 7943 150 
metadata(0):
assays(1): counts
rownames(7943): ELMO2 PNMA1 ... SLC39A6 CTSC
rowData names(0):
colnames(150): mel_cell1 mel_cell2 ... mel_cell149 mel_cell150
colData names(1): celltype_CHETAH
reducedDimNames(1): TSNE
spikeNames(0):

head(input_mel$celltype_CHETAH)
    mel_cell1     mel_cell2     mel_cell3     mel_cell4     mel_cell5 
"CD8 T cell" "Endothelial"       "Node7"  "CD8 T cell" "reg. T cell" 
    mel_cell6 
      "Node1" 


The output
CHETAH returns the input object, but added:
  • input$celltype_CHETAH:a named character vector that can directly be used in any other workflow/method.
  • “hidden” int_colData and int_metadata, not meant for direct interaction, but which can all be viewed and interacted with using: PlotCHETAH and CHETAHshiny

Standard plots
我们可以使用PlotCHETAH函数查看CHETAH的分类结果,此函数可绘制分类树和t-SNE(或其他提供的降维)图。在这些图中,最终类型(final types)或中间类型(intermediate types)都带有颜色,非彩色类型(non-colored types)以灰度表示。
To plot the final types:
[AppleScript] 纯文本查看 复制代码
PlotCHETAH(input = input_mel)



请注意,每种类型的“NodeX”对应于编号为X的节点,而“Unassigned”类型对应于节点0Conversely, to color the intermediate types:
[AppleScript] 纯文本查看 复制代码
PlotCHETAH(input = input_mel, interm = TRUE)


If you would like to use the classification, and thus the colors, in another package (e.g. Seurat), you can extract the colors using:
[AppleScript] 纯文本查看 复制代码
colors <- PlotCHETAH(input = input_mel, return_col = TRUE)
colors
       Unassigned             Node1             Node2             Node3 
         "gray90"          "gray86"          "gray82"          "gray78" 
            Node4             Node5             Node6             Node7 
         "gray74"          "gray70"          "gray66"          "gray62" 
            Node8        Fibroblast            B cell        Macrophage 
         "gray58"            "blue"            "gold"           "cyan3" 
      Endothelial         Dendritic              Mast        CD4 T cell 
           "navy"     "forestgreen"          "orange" "darkolivegreen3" 
      reg. T cell        CD8 T cell 
          "brown"           "green" 


CHETAHshiny
CHEETAH的分类结果及其他输出(如配置文件和置信度得分)可以在一个shiny应用程序中显示,该应用程序可以轻松、交互式地进行分类分析。
Here you can view:
  • the confidence of all assignments
  • the classification in an interactive window
  • the genes used by CHETAH, an it’s expression in the input data
  • a lot more

[AppleScript] 纯文本查看 复制代码
CHETAHshiny(input = input_mel)


Changing classification
Confidence score
CHETAH计算输入细胞到任一节点分支的每次分配的置信度得分,这些置信度得分:
  • 值介于0到2之间
  • 通常在0到1之间
  • 0表示细胞分配的置信度最低,1表示高的置信度。

CHETAH的默认置信度阈值为0.1。
这意味着只要将一个细胞分配给一个分支,并且该分配的置信度小于0.1,分类就会在该节点停止。

可以调整置信度阈值,以便将更多或更少的细胞分类为最终类型:
  • 如果将置信度阈值设置为0,所有的输入细胞将会被分类为最终类型。请注意,此分类可能很嘈杂,并且可能包含错误的分类。
  • 使用0.2、0.3、0.4等阈值会将逐渐减少被分类为最终类型的细胞的数量,其余的细胞在分类树中所有节点上的置信度越来越高。

For example, to only classify cells with very high confidence:
[AppleScript] 纯文本查看 复制代码
input_mel <- Classify(input = input_mel, 0.8)
PlotCHETAH(input = input_mel, tree = FALSE)


Conversely, to classify all cells:
[AppleScript] 纯文本查看 复制代码
input_mel <- Classify(input_mel, 0)
PlotCHETAH(input = input_mel, tree = FALSE)


Renaming types
CHETAH包提供了RenameBelowNode函数,可以对分类树中的细胞类型进行重命名。
对于示例数据,假设我们对T细胞的不同子类型(在Node6和Node7下)都不感兴趣,我们可以通过运行以下命令将所有这些细胞命都名为“T细胞”:
[AppleScript] 纯文本查看 复制代码
input_mel <- RenameBelowNode(input_mel, whichnode = 6, replacement = "T cell")
PlotCHETAH(input = input_mel, tree = FALSE)


To reset the classification to its default, just run Classify again:
[AppleScript] 纯文本查看 复制代码
input_mel <- Classify(input_mel) ## the same as Classify(input_mel, 0.1)
PlotCHETAH(input = input_mel, tree = FALSE)


Creating a reference
Step 0: Obtain a reference.
  • 下载或使用自己的scRNA-seq数据集,这些数据集具有已知可用的细胞类型标签。
  • 可以在此处下载包含肿瘤数据的所有主要细胞类型的肿瘤微环境参考数据集
  • 构建SingleCellExperiment对象。
Step 1: good reference characteristics
CHETAH可以使用任何scRNA-seq数据集作为参考,但是使用不同的参考数据集会极大地影响最终的细胞分类结果。
以下一般规则可以用于选择和创建参考数据集:
  • 通过使用来自相同生物学类型或至少由处于相同细胞状态的细胞组成的参考和输入数据集,可以获得更好的分类结果。例如。对于PBMC的输入数据集,可以使用骨髓参考数据集,但是由于这些细胞大多数naive或前体细胞,这可能会对分类产生负面影响。在这种情况下,另一个PBMC的数据集将发挥最佳作用。
  • 参考数据集细胞注释结果的好坏会直接影响最终的细胞分类结果。参考数据集的细胞类型标签越准确,最终的分类结果会越好。
  • 参考数据的数据越稀疏,则需要更多的参考细胞来创建可靠的参考。对于高覆盖率的Smart-Seq2数据,每个细胞类型仅需要10-20个细胞。而对于稀疏的10X Genomics数据,通常需要100多个细胞可获得最佳结果。

为了减少使用大参考数据集的计算时间,我们尝试将每种细胞类型下采样为100-200个细胞。对于具有细胞类型元数据“celltypes”的SingleCellExperiment对象的“ref”参考数据集,可以通过以下方式完成:
[AppleScript] 纯文本查看 复制代码
cell_selection <- unlist(lapply(unique(ref$celltypes), function(type) {
    type_cells <- which(ref$celltypes == type)
    if (length(type_cells) > 200) {
        type_cells[sample(length(type_cells), 200)]
    } else type_cells
}))

ref_new <- ref[ ,cell_selection]


Step 1: normalization
对于输入数据,CHETAH不需要进行了标准化的输入数据,但是参考数据必须事先进行标准化处理。此示例中使用的参考数据已经做了标准化处理。但出于示例的目的,我们还是对此数据集进行标准化处理:
[AppleScript] 纯文本查看 复制代码
assay(headneck_ref, "counts") <- apply(assay(headneck_ref, "counts"), 2, 
                                       function(column) log2((column/sum(column) * 100000) + 1))

head(headneck_ref@assays$data$counts)[,1:5]
             hn_cell1 hn_cell2 hn_cell3 hn_cell4 hn_cell5
ELMO2        0.000000        0 0.000000 4.249018 5.939126
PNMA1        0.000000        0 0.000000 5.748898 0.000000
MMP2         0.000000        0 6.716108 5.734995 6.485250
TMEM216      0.000000        0 0.000000 0.000000 0.000000
TRAF3IP2-AS1 2.206329        0 3.417146 2.121777 5.704061
ERCC5        7.551864        0 6.326401 0.000000 0.000000


Step 2: discaring of house-keeping genes
当然,在drop-out率较高的情况下,CHETAH可能会受到高表达基因(具有高度可变性)的影响。根据我们的经验,主要是核糖体蛋白基因会产生这种影响。因此,我们在此处删除那些“核糖体”基因。
[AppleScript] 纯文本查看 复制代码
ribo <- read.table("~/ribosomal.txt", header = FALSE, sep = '\t')
headneck_ref <- headneck_ref[!rownames(headneck_ref) %in% ribo[,1], ]


Step 3: Reference QC
CHETAH的性能在很大程度上取决于参考数据集的质量。
参考数据集的质量会受到以下因素的影响:
  • scRNA-seq数据的质量
  • 细胞类型标签的准确性

要检查CHETAH在参考数据集中区分不同细胞类型的能力,我们可以运行CorrelateReference函数,也可以使用ClassifyReference函数进行查看。
1)CorrelateReference
CorrelateReference函数会对两种细胞类型的每种组合,寻找两者之间具有最高差异倍数变化的基因,并使用这些基因对它们进行相关性分析。如果参考数据集是好的,则所有细胞类型的关联性都会变差甚至更好,并且会反相关。
[AppleScript] 纯文本查看 复制代码
CorrelateReference(ref_cells = headneck_ref)
#> Running... in case of 1000s of cells, this may take a couple of minutes



在上图中,大多数细胞的类型都是可以区分的:许多细胞类型不相关或反相关。但是,某些细胞类型非常相似,如调节性和CD4阳性的T细胞,或CD4和CD8阳性的T细胞,这些细胞可能很难在输入数据中区分开。

2)ClassifyReference
我们还可以使用ClassifyReference函数检查CHETAH是否可以区分参考数据中的细胞类型。此函数使用参考数据对本身进行分类,如果CHETAH与参考数据集分类的结果较好,则分类中应该几乎没有混淆的细胞类型,即所有A型细胞应被分类为A型。
[AppleScript] 纯文本查看 复制代码
ClassifyReference(ref_cells = headneck_ref)
#> Preparing data....
#> Running analysis...


在此图中,每一行是原始细胞类型标签,每一列是在CHETAH分类后分配的细胞标签。正方形的颜色和大小指示行类型细胞中的哪一部分被分类为列类型。如,第4行第2列显示5-10%的CD4 T细胞被归类为调节性T细胞。
在此参考数据集中,两种细胞类型之间的混合比例不会超过10%。另外,低百分比的细胞被分类为中间类型。大多数混淆发生在T细胞的不同亚型之间。在这种情况下,我们应注意,这些细胞类型标签具有最高的互换机会。

Optimizing the classification
CHETAH包经过优化,可以在大多数的分析中获得良好的结果,但是也可能会发生分类不完善的情况。如果CHETAH无法提供所需的输出(分类的细胞太少,视觉上的随机分类等),我们可以执行以下步骤(按此顺序)进行完善:
  • 检查我们的参考数据集是否正确创建(请参见上文,这是最重要的步骤!)
  • 如果分类的细胞太少,则将置信度阈值降低到0.05或0.01(注意假阳性!始终检查结果是否有意义。)
  • 使用input[!(grepl("^RP", rownames(input))), ]是不完善的方法,但是非常快捷。
  • 尝试使用其他数量的基因进行分类(n_genes参数)。默认为200,但有时100或(在稀疏数据中)500可以产生更好的结果。
  • 寻找其他更好的参考数据集
  • 尝试其他可用的对scRNA-seq数据进行细胞类型分类注释的方法。

参考来源:http://www.bioconductor.org/pack ... ntroduction.htmlEND

▼更多精彩推荐,请关注我们▼


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
加油啦。
回复

使用道具 举报

草履虫

Rank: 2

主题
0
注册时间
2020.6.29
在线时间
0 小时

发表于 2020.7.6 18:11:43 | 显示全部楼层
这个很实用欸
回复 支持 反对

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2019.8.6
在线时间
5 小时

发表于 2020.7.25 09:26:11 | 显示全部楼层
有点难懂
回复

使用道具 举报

帝王蝶

Rank: 4

主题
12
注册时间
2020.4.30
在线时间
18 小时

突出贡献


发表于 2020.7.31 22:27:34 | 显示全部楼层
今天也是元气满满的一天啊!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2020.2.12
在线时间
17 小时

活跃会员


发表于 2020.8.4 16:02:22 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

草履虫

Rank: 2

主题
0
注册时间
2020.8.19
在线时间
0 小时

发表于 2020.8.19 16:00:24 | 显示全部楼层
R语言好处多,就是比较难
回复 支持 反对

使用道具 举报

帝王蝶

Rank: 4

主题
5
注册时间
2020.3.21
在线时间
14 小时

发表于 2020.9.6 16:14:19 | 显示全部楼层
看起来比较难
终于周五了
回复 支持 反对

使用道具 举报

帝王蝶

Rank: 4

主题
12
注册时间
2020.4.30
在线时间
18 小时

突出贡献


发表于 2020.9.6 16:38:02 | 显示全部楼层
今天也是元气满满的一天啊!
回复

使用道具 举报

帝王蝶

Rank: 4

主题
12
注册时间
2020.4.30
在线时间
18 小时

突出贡献


发表于 2020.9.12 10:29:49 | 显示全部楼层
今天也是元气满满的一天啊!
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表