基因组学

MAKER的并行化初探

MAKER并行分为两种，一种基于MPI，运行方式为mpiexec -n 线程数 maker, 一种是在同一个项目中运行多次maker。前者需要在安装MAKER时进行设置，后者相当于你手动按照染色体数目进行拆分，然后分开运行MAKER。本片主要介绍基于MPI的并行策略。

下图是MAKER是基于MPI的并行化流程。分为三个层次，contig水平，更小的片段水平，不同程序的并行。

contig水平就是对每条contig进行分别注释。而由于每个contig的长度不一，当较短的contig运行结束后，我们还需要等待较长的contig结束。因此，我们还需要将contig继续拆分，让每条contig分成更小的片段，成为基本的分析单元。对于每个分析单元，还可以考虑将不同程序进行并行。

Fig1

在MAKER实际运行时，如果用gotop进行查看你会发现系统中会出现很多maker进程(以mpiexec -n 10为例)

Fig2

但是你发现只有少量的和注释有关的程序(如blastx)。如果用ps aux | grep maker查找和maker有关的进程，你会发现大量的maintain.pl进程，并且这些进程后面跟着一大串你不认识的符号

1
/usr/bin/perl /opt/biosoft/maker/bin/../lib/File/maintain.pl 46535 30 %04%09%08%31%3...

为了理解MAKER的并行，你需要通过pstree -ap [mpiexec进程ID]更细致的了解MAKER的进程树

Fig3

从中你可以发现，因为我们用mpiexec启动了10个maker进程，所以hydra_pmi_proxy有10个maker子进程。对于这10个maker子进程，每个进程都至少会有一个maintain.pl。通过阅读maintain.pl的源代码，我们可以得知该程序后接参数中的46535是PID(进程号), 30是sleep time, 最后一个是URI编码字符，是利用Storable::freeze持久化的数据结构，可以通过如下的代码进行解码（命名为decode.pl)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#!/usr/bin/env perl
use warnings;
use strict;

use URI::Escape;
use Storable;
use vars qw($LOCK);

my $serial = shift;

$serial = uri_unescape($serial);
$LOCK = Storable::thaw($serial);

print "$LOCK->{lock_file} \n"

可以批量对maintain.pl里的信息进行解码

1
ps aux | grep maintain.pl | grep -v grep | awk '{print $15}' | xargs -i perl decode.pl {}

因此，那么多的maker并不是在话说，而是在负责任务调度，差不多一个任务会有3个maker进程保驾护航。

回到之前的进程数，这里有一些maker进程下会跟着一些正在运行的命令，blastx, exonerate

Fig4

而有些只有maker和maintain.pl进程。

Fig5

如果等待一会，再次运行pstree -ap [mpiexec进程ID]，你会发现这两种状态会发生转换，这说明任务启动也需要一段时间。

Linux

Singularity和MPI应用

MPI(Message Passin Interface)广泛应用于高性能服务器中，可用于单系统多节点或者多个计算平台间通讯，目前主流的两个开源软件分别是OpenMPI和MPICH。Singularity同时支持这两个开源工具，本篇教程介绍如何在Singularity容器开发和运行MPI程序。

在Singularity容器是执行MPI程序最常见的方法就是使用宿主机中的MPI。由于我们会同时使用宿主机的MPI和容器中的MPI，因此这被称为Host MPI或者混合模型。考虑到大部分高性能计算平台并不都支持宿主机和容器之间的文件系统共享，因此我们不讨论如何将存储盘挂载到容器，从而使用容器中宿主机MPI。

混合途径的基本想法就是，当你要执行SIngularity容器的MPI代码时，你会使用mpiexec类似命令去调用singularity命令。容器外部的MPI进程会和容器内的MPI进行写作，容器后的MPI代码会实例化任务。

Open MPI/Singularity的工作流程如下

  1. 用户在shell中启动MPI(mpirun, mpiexec)
  2. Open MPI接着调用进程管理守护进程(process management daemon, ORTED)
  3. ORTED进程启动所需的Singularity容器
  4. Singularity构建容器和命名空间环境
  5. Singularity启动容器内的MPI应用
  6. MPI进程启动，加载Open MPI库
  7. Open MPI库通过进程管理接口(Process Management, PMI)连接回ORTED进程

此时在容器运行MPI就和直接宿主机运行MPI程序一样。优势在于它能够整合到Slurm等资源管理工具，并且和之前运行MPI应用一样。但是你需要保证容器中的MPI和宿主机的MPI版本兼容，也就是在构建MPI容器时保证和系统的MPI一致。

我们举个简单的例子。假设我们已经有了一个叫做mpitest.cMPI应用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>

int main (int argc, char **argv) {
        int rc;
        int size;
        int myrank;

        rc = MPI_Init (&argc, &argv);
        if (rc != MPI_SUCCESS) {
                fprintf (stderr, "MPI_Init() failed");
                return EXIT_FAILURE;
        }

        rc = MPI_Comm_size (MPI_COMM_WORLD, &size);
        if (rc != MPI_SUCCESS) {
                fprintf (stderr, "MPI_Comm_size() failed");
                goto exit_with_error;
        }

        rc = MPI_Comm_rank (MPI_COMM_WORLD, &myrank);
        if (rc != MPI_SUCCESS) {
                fprintf (stderr, "MPI_Comm_rank() failed");
                goto exit_with_error;
        }

        fprintf (stdout, "Hello, I am rank %d/%d", myrank, size);

        MPI_Finalize();

        return EXIT_SUCCESS;

 exit_with_error:
        MPI_Finalize();
        return EXIT_FAILURE;
}

MPI只是一个库接口，里面的函数能够被其他编程语言所调用，所以支持Fortran, C, Python, R等。

接下来根据宿主机的MPI版本，编写一个定义文件。假如我们宿主机的MPI是MPICH，那么定义文件(命名为mpitest)如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
Bootstrap: docker
From: ubuntu:latest

%files
    mpitest.c /opt

%environment
    export MPICH_DIR=/opt/mpich-3.3
    export SINGULARITY_MPICH_DIR=$MPICH_DIR
    export SINGULARITYENV_APPEND_PATH=$MPICH_DIR/bin
    export SINGULAIRTYENV_APPEND_LD_LIBRARY_PATH=$MPICH_DIR/lib

%post
    echo "Installing required packages..."
    apt-get update && apt-get install -y wget git bash gcc gfortran g++ make

    # Information about the version of MPICH to use
    export MPICH_VERSION=3.3
    export MPICH_URL="http://www.mpich.org/static/downloads/$MPICH_VERSION/mpich-$MPICH_VERSION.tar.gz"
    export MPICH_DIR=/opt/mpich

    echo "Installing MPICH..."
    mkdir -p /tmp/mpich
    mkdir -p /opt
    # Download
    cd /tmp/mpich && wget -O mpich-$MPICH_VERSION.tar.gz $MPICH_URL && tar xzf mpich-$MPICH_VERSION.tar.gz
    # Compile and install
    cd /tmp/mpich/mpich-$MPICH_VERSION && ./configure --prefix=$MPICH_DIR && make install
    # Set env variables so we can compile our application
    export PATH=$MPICH_DIR/bin:$PATH
    export LD_LIBRARY_PATH=$MPICH_DIR/lib:$LD_LIBRARY_PATH
    export MANPATH=$MPICH_DIR/share/man:$MANPATH

    echo "Compiling the MPI application..."
    cd /opt && mpicc -o mpitest mpitest.c

如果是Open MPI,

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
Bootstrap: docker
From: ubuntu:latest

%files
    mpitest.c /opt

%environment
    export OMPI_DIR=/opt/ompi
    export SINGULARITY_OMPI_DIR=$OMPI_DIR
    export SINGULARITYENV_APPEND_PATH=$OMPI_DIR/bin
    export SINGULAIRTYENV_APPEND_LD_LIBRARY_PATH=$OMPI_DIR/lib

%post
    echo "Installing required packages..."
    apt-get update && apt-get install -y wget git bash gcc gfortran g++ make file

    echo "Installing Open MPI"
    export OMPI_DIR=/opt/ompi
    export OMPI_VERSION=4.0.1
    export OMPI_URL="https://download.open-mpi.org/release/open-mpi/v4.0/openmpi-$OMPI_VERSION.tar.bz2"
    mkdir -p /tmp/ompi
    mkdir -p /opt
    # Download
    cd /tmp/ompi && wget -O openmpi-$OMPI_VERSION.tar.bz2 $OMPI_URL && tar -xjf openmpi-$OMPI_VERSION.tar.bz2
    # Compile and install
    cd /tmp/ompi/openmpi-$OMPI_VERSION && ./configure --prefix=$OMPI_DIR && make install
    # Set env variables so we can compile our application
    export PATH=$OMPI_DIR/bin:$PATH
    export LD_LIBRARY_PATH=$OMPI_DIR/lib:$LD_LIBRARY_PATH
    export MANPATH=$OMPI_DIR/share/man:$MANPATH

    echo "Compiling the MPI application..."
    cd /opt && mpicc -o mpitest mpitest.c

接着我们从定义文件中构建出sig文件, singularity build mpitest.sig mpitest

最后，执行Singularity容器的MPI应用的标准方法就是在host中运行mpirun命令，它会启动Singularity容器，并最终运行其中的MPI程序。

1
2
3
mpirun -n 4 singularity exec mpitest.sig /opt/mpitest
# 输出信息如下
# Hello, I am rank 0/4Hello, I am rank 1/4Hello, I am rank 2/4Hello, I am rank 3/4

如果是在SLURM这类任务提交系统中，我们需要批处理脚本来执行MPI引用，举个例子

1
2
3
4
5
6
#!/bin/bash
#SBATCH --job-name singularity-mpi
#SBATCH -N $NNODES # total number of nodes
#SBATCH --time=00:05:00 # Max execution time

mpirun -n $NP singularity exec /var/nfsshare/gvallee/mpich.sif /opt/mpitest

执行方式如下

1
sbatch my_job.sh

参考资料: https://sylabs.io/guides/3.3/user-guide/mpi.html

R

使用ArchR分析单细胞ATAC-seq数据(第九章)

第9章 ArchR的拟混池重复

因为scATAC-seq数据本质上只有两种值，也就是说每个位点要么开放要么不开放，所以你会发现某些情况下无法使用单个细胞进行数据分析。此外，许多我们想要做的分析也依赖于重复才能计算统计显著性。对于单细胞数据，我们通过构建拟混池重复(pseudo-bulk replicates)来解决该问题。所谓的拟混池(pesudo-bulk)指的就是将单细胞进行合并模拟成混池测序的ATAC-seq实验得到的数据。ArchR为每个目标细胞分组构建多个拟混池样本，也就得到了拟混池重复。这个模拟过程背后的假设是将单细胞进行合并结果和实际的混池结果非常接近，以至于不需要在乎背后的差异。这些细胞分组通常都是来自于单个细胞类群，或者是一直细胞类型对应的可能聚类。我们这一章会介绍如何使用ArchR生成这些拟混池重复。

9.1 ArchR如何构建拟混池重复

ArchR使用分级优先法(tiered priority approach)构建拟混池重复。使用者定义: i)最小和最大重复数(minReps/maxReps), ii)每个重复的最少和最多细胞数(minCells/maxCells), iii)缺少足够细胞用于构建足够重复时的采样率(Sampling Ratio)。举个例子，当采样率等于0.8时，每个重复中80%的细胞来自于无放回抽样(这会导致重复间出现有放回抽样)。在这种情况下，多个重复中可能有一些相同的细胞，但是为了能在缺少足够细胞的分组里得到拟混池重复，这是必要的牺牲。

我们的拟混池重复生成过程可以用如下的决策树进行描述

pseudobulkReplicate_DecisionTree

流程有些复杂，这里概述下这个流程中的几个关键注意事项。

首先，用户确定细胞的分组，ArchR通常会称之为聚类，我们可以在第五章分析分析得到。

其次，对于每一组，ArchR尝试去创建理想的拟混池重复。所谓理想的拟混池重复，指的是每个重复只有一个样本构成，这就要求每个样本得保证足量的细胞数。这样保证了重复间的样本多样性和生物学变异，是ArchR所期望得到的最佳结果，但这一过程实际上会有5种可能的结果，根据ArchR的偏好排序如下

  1. 数据拥有足够多的样本(至少要等于maxRep定义的重复)，且每个样本的细胞数都大于 minCell, 于是每个样本都可以作为拟混池重复，每个重复的细胞都来自于同一个样本。
  2. 一些样本的细胞数超过minCell,因此能够单独形成一个重复。剩下的重复则是通过将余下的细胞进行混合，然后通过无放回抽样得到。
  3. 数据里没有一个样本的细胞数超过minCell, 但是总细胞数超过minCells * minReps。因此将所有的细胞进行混合，然后进行无放回抽样，抽样时不考虑细胞来源。
  4. 一个细胞分组中的总细胞数低于 minCells * minReps，但是大于minCells / Sample Ratio。此时单个样本的构建采取无放回抽样，重复间则需要有放回抽样，降低多个拟混池重复间的相同细胞数。
  5. 一个细胞分组中的总细胞数低于 minCells / Sample Ratio 。这意味着我们必须在单个重复和跨重复中都采取有放回抽样策略。这是最糟糕的情况，后续在使用这些拟混池重复做下游分析分析要特别小心。后续可以通过设置ArchR的minCells参数进行淘汰。

我们使用如下的数据集阐述这一过程

1
2
3
4
5
6
7
Sample  Cluster1  Cluster2  Cluster3  Cluster4  Cluster5
A       800       600       900       100       75
B       1000      50        400       150       25
C       600       900       100       200       50
D       1200      500       50        50        25
E       900       100       50        150       50
F       700       200       100       100       25

我们设置的参数为minRep=3,maxRep=5, minCells=300, maxCells=1000和sampleRatio=0.8，也就是最少有3个重复，最多是5个重复，每个重复至少有300个细胞，最多是1000个细胞，当细胞数不满住要求，抽样率设置为0.8.

9.1.1 Cluster1

对于Cluster1, 我们的6个样本(大于maxRep)的细胞数都高于minCells(300)。这是最理想的情况，对应上述第一种情况，我们将会得到5个拟混池重复，保证每个重复都来自独立的样本。

1
2
3
4
5
Rep1 = 800 cells from SampleA
Rep2 = 1000 cells from SampleB
Rep3 = 1000 cells from SampleD
Rep4 = 900 cells from SampleE
Rep5 = 700 cells from SampleF

对于这些重复，我们需要注意两个事情:(1) 因为我们的样本数足够多，能够保证每个重复都来自独立的样本，所以可以淘汰其中细胞数最少的SampleC。(2)由于maxCells设置为1000，因此最多只能有1000个细胞。

9.1.2 Cluster2

对于Cluster2, 我们有3个样本的细胞数超过minCells, 另外3个样本的细胞数都不够。这对应上述第二种情况，我们会以如下的方法构建拟混池重复。

1
2
3
4
Rep1 = 600 cells from SampleA
Rep2 = 900 cells from SampleC
Rep3 = 500 cells from SampleD
Rep4 = 350 cells [50 cells from SampleB + 100 from SampleE + 200 from SampleF]

在这个例子中，Rep4由其他几个样本的细胞混合后通过无放回抽样得到

9.1.3 Cluster3

对于Cluster3，我们只有两个样本超过minCells, 不满足minReps。但是如果我们将剩余的样本的细胞进行混合形成额外的重复，它的细胞数就超过了minCells。最终我们得到了3个拟混池重复，对应上述的情况3。我们将得到如下重复

1
2
3
Rep1 = 900 cells from SampleA
Rep2 = 400 cells from SampleB
Rep3 = 250 cells [100 cells from SampleC + 50 from SampleD + 50 from SampleE + 50 from SampleF]

和Cluster2类似，Cluster3的Rep3由其他几个样本的细胞混合后通过无放回抽样得到

9.1.4 Cluster4

对于Cluster4，总细胞数是570个，小于minCells * minReps(900). 在这个情况下，我们无法保证有足够多的细胞通过无放回抽样的方式保证每个重复都有最小的细胞数。但是，总的细胞数依旧依旧大于minCells / sampleRatio(375个细胞)，这意味着每个重复中细胞可以来自于无放回抽样，重复之间的细胞需要放回抽样。这对应着上述的情况4，我们将得到如下重复

1
2
3
Rep1 = 300 cells [250 unique cells + 25 cells overlapping Rep2 + 25 cells overlapping Rep3]
Rep2 = 300 cells [250 unique cells + 25 cells overlapping Rep1 + 25 cells overlapping Rep3]
Rep3 = 300 cells [250 unique cells + 25 cells overlapping Rep1 + 25 cells overlapping Rep2]

在这个情况中，ArchR会尽可能降低任意两个拟混池重复的相同细胞。

9.1.5 Cluster5

对于Cluster5，总共是250个细胞，同时小于minCells * minReps（900）和minCells / sampleRatio(375). 这意味着每个样本都需要有放回的抽样，重复之间也需要有放回抽样，才能得到拟混池重复。这是上述说到的第5种情况，是其中最糟糕的情况。对于这类拟混池重复，在后续的分析中需要谨慎使用。我们将得到如下的重复:

1
2
3
Rep1 = 300 cells [250 unique cells + 25 cells overlapping Rep2 + 25 cells overlapping Rep3]
Rep2 = 300 cells [250 unique cells + 25 cells overlapping Rep1 + 25 cells overlapping Rep3]
Rep3 = 300 cells [250 unique cells + 25 cells overlapping Rep1 + 25 cells overlapping Rep2]

9.2 构建拟混池重复

通过上一节了解ArchR构建拟混池重复的逻辑后，我们就可以开始实际操作了。在ArchR中，我们通过调用addGroupCoverages()函数来构建拟混池重复。它的关键参数groupBy,定义了拟混池重复需要使用的分组。这里我们用的是上一章scRNA-seq数据标记的细胞类型，也就是Cluster2

1
projHeme4 <- addGroupCoverages(ArchRProj = projHeme3, groupBy = "Clusters2")

得到这些拟混池重复后，我们就能从数据中鉴定peak了。就像之前所说的，我们不希望使用所有的细胞鉴定peak，而是单独根据每一组细胞（例如聚类）单独鉴定peak，这样才有可能分析出不同组的特异性peak。这一章得到数据就为后续鉴定peak提供了良好的开始。

R

使用ArchR分析单细胞ATAC-seq数据(第八章)

第8章: 使用scRNA-seq定义cluster类型

除了使用基因得分定义细胞类群以外，ArchR还能整合scRNA-seq数据。通过将scATAC-seq数据里的基因得分矩阵和scRNA-seq数据的基因表达量矩阵进行对比，ArchR就能将scATAC-seq的细胞比对到scRNA-seq的细胞，实现两种数据的整合。之后，我们借助scRNA-seq数据已经定义的细胞类群，或者整合后的scRNA-seq的基因表达量来注释细胞类群。在代码内部，我们调用了Seurat::FindTransferAnchors()，从而实现两种数据集之间的比较。当然，ArchR不是简单地对函数进行封装，而是在此基础上通过对数据的拆分，实现并行计算，从而能将该流程应用到更大规模的细胞中。

img

该整合过程实际上会寻找scATAC-seq和scRNA-seq两者中最相似的细胞，然后将scRNA-seq中对应的细胞表达量赋值给scATAAC-seq细胞。最后，scATAC-seq中每个细胞都会有基因表达量特征，能被用于下游分析。这一章会阐述如何利用该信息定义细胞类型，之后会介绍如何使用连接的scRNA-seq数据做更加复杂的分析，例如识别预测的顺式调控元件。我们相信由于多组学单细胞谱的商业化，这类整合分析将会越来越多。同时，在ArchR中使用公共数据里匹配细胞类型的scRNA-seq数据或者自己使用目标样本得到的scRNA-seq数据也能加强scATAC-seq分析。

8.1 scATAC-seq细胞和scRNA-seq细胞跨平台连接

为了能将我们教程中的scATAC-seq数据与其匹配的scRNA-seq数据进行整合，我们将使用 Granja* et al (2019) 里的造血细胞scRNA-seq数据。

scRNA-seq数据以 RangedSummarizedExperiment对象保存，大小为111MB。此外，ArchR还接受未经修改的Seurat对象作为整合流程的输入。我们使用download.file下载数据

1
2
3
4
5
6
if(!file.exists("scRNA-Hematopoiesis-Granja-2019.rds")){
    download.file(
        url = "https://jeffgranja.s3.amazonaws.com/ArchR/TestData/scRNA-Hematopoiesis-Granja-2019.rds",
        destfile = "scRNA-Hematopoiesis-Granja-2019.rds"
    )
}

下载之后，我们使用readRDS进行读取，并查看该对象

1
2
3
4
5
6
7
8
9
10
11
12
13
seRNA <- readRDS("scRNA-Hematopoiesis-Granja-2019.rds")
seRNA
## class: RangedSummarizedExperiment
## dim: 20287 35582
## metadata(0):
## assays(1): counts
## rownames(20287): FAM138A OR4F5 … S100B PRMT2
## rowData names(3): gene_name gene_id exonLength
## colnames(35582): CD34_32_R5:AAACCTGAGTATCGAA-1
## CD34_32_R5:AAACCTGAGTCGTTTG-1 …
## BMMC_10x_GREENLEAF_REP2:TTTGTTGCATGTGTCA-1
## BMMC_10x_GREENLEAF_REP2:TTTGTTGCATTGAAAG-1
## colData names(10): Group nUMI_pre … BioClassification Barcode

从输出信息中，我们可以发现它里面有基因表达量的count矩阵和对应的元信息。

元信息列里的BioClassification记录着scRNA-seq数据中每个细胞对应的细胞类型分类

1
2
3
4
5
colnames(colData(seRNA))
# [1] "Group"             "nUMI_pre"          "nUMI"             
# [4] "nGene"             "initialClusters"   "UMAP1"            
# [7] "UMAP2"             "Clusters"          "BioClassification"
# [10] "Barcode"

使用table()，我们可以看到scRNA-seq细胞类型每一群的细胞数

1
2
3
4
5
6
7
8
9
10
11
12
13
table(colData(seRNA)$BioClassification)
#        01_HSC 02_Early.Eryth  03_Late.Eryth  04_Early.Baso    05_CMP.LMPP 
#          1425           1653            446            111           2260 
#      06_CLP.1         07_GMP    08_GMP.Neut         09_pDC         10_cDC 
#           903           2097           1050            544            325 
#11_CD14.Mono.1 12_CD14.Mono.2   13_CD16.Mono         14_Unk       15_CLP.2 
#          1800           4222            292            520            377 
#      16_Pre.B           17_B      18_Plasma       19_CD8.N      20_CD4.N1 
#           710           1711             62           1521           2470 
#     21_CD4.N2       22_CD4.M      23_CD8.EM      24_CD8.CM          25_NK 
#          2364           3539            796           2080           2143 
#        26_Unk 
#           161

我们后续会用到两种整合方法。第一种是无约束整合，我们对scATAC-seq实验里的细胞不作任何假设，直接尝试将这些细胞和scRNA-seq实验里的任意细胞进行配对。这是一种初步可行方案，后续会根据这一步得到的结果，对整合步骤进行约束来提升跨平台配对的质量。第二种方法是约束整合，即利用对细胞类型的先验知识限制搜索范围。举个例子，如果我们知道scATAC-seq中的Cluster A,B,C对应着三种不同的T细胞，scRNA-seq中的Cluster X,Y对应着两种不同的T细胞，我们告诉ArchR只需要尝试将scATAC-seq中的Cluster A,B,C跟scRNA-seq中的Cluster X,Y进行配对。下面，我们将先以无约束整合初步地鉴定每一种聚类的类型，然后根据分析结果做更加细致的约束整合。

8.1.1 无约束整合

我们使用addGeneIntegrationMatrix()对scATAC-seq和scRNA-seq数据进行整合。正如之前所提到的，该函数的seRNA参数接受Seurat或RangedSummarizedExperiment对象作为输入。因为第一轮是探索性质的无约束整合，因此，我们不会将结果保存在Arrow文件中(addToArrow=FALSE)。整合后的矩阵将会根据matrixName进行命名，存放在ArchRProject中。该函数的其他参数对应cellColData中列名用于存放额外的信息，nameCell存放scRNA-seq中匹配的细胞ID，nameGroup存放scRNA-seq细胞中的分组ID，nameScore存放跨平台整合得分。

1
2
3
4
5
6
7
8
9
10
11
12
projHeme2 <- addGeneIntegrationMatrix(
    ArchRProj = projHeme2, 
    useMatrix = "GeneScoreMatrix",
    matrixName = "GeneIntegrationMatrix",
    reducedDims = "IterativeLSI",
    seRNA = seRNA,
    addToArrow = FALSE,
    groupRNA = "BioClassification",
    nameCell = "predictedCell_Un",
    nameGroup = "predictedGroup_Un",
    nameScore = "predictedScore_Un"
)

无约束整合的结果可能并不准确，但是为后续更加精细的约束分析奠定了基础。

8.1.2 约束整合

现在我们已经有了初步的无约束整合结果，我们就有了大致的细胞类型分布情况，接着就是优化整合结果。

因为我们教程里的数据来自于造血细胞，我们将会在理想状态下将类似的细胞整合在一起。首先，我们先确认scRNA-seq里细胞类型在我们的scATAC-seq聚类中分布情况。这一步的目标是使用无约束整合的方法找到scATAC-seq和scRNA-seq中和T细胞和NK细胞对应的聚类，后续会用到该信息进行约束整合。具体操作为，我们创建一个confusionMatrix,并关注Cluster和predictedGroup_Un的交叉部分中scRNA-seq的细胞类型。

1
2
3
cM <- as.matrix(confusionMatrix(projHeme2$Clusters, projHeme2$predictedGroup_Un))
preClust <- colnames(cM)[apply(cM, 1 , which.max)]
cbind(preClust, rownames(cM)) #Assignments

输出信息如下，展示了12个scATAC-seq聚类中对应最优可能的scRNA-seq细胞类型。

1
2
3
4
5
6
7
8
9
10
11
12
13
#      preClust              
# [1,] "03_Late.Eryth"  "C10"
# [2,] "20_CD4.N1"      "C8" 
# [3,] "16_Pre.B"       "C3" 
# [4,] "08_GMP.Neut"    "C11"
# [5,] "17_B"           "C4" 
# [6,] "11_CD14.Mono.1" "C1" 
# [7,] "01_HSC"         "C12"
# [8,] "22_CD4.M"       "C9" 
# [9,] "09_pDC"         "C5" 
# [10,] "25_NK"          "C7" 
# [11,] "12_CD14.Mono.2" "C2" 
# [12,] "06_CLP.1"       "C6"

首先，我们检查在无约束整合中用到的scRNA-seq数据里细胞类型标签。

1
2
3
4
5
6
7
unique(unique(projHeme2$predictedGroup_Un))
# [1] "08_GMP.Neut"    "25_NK"          "16_Pre.B"       "06_CLP.1"      
# [5] "07_GMP"         "11_CD14.Mono.1" "04_Early.Baso"  "22_CD4.M"      
# [9] "03_Late.Eryth"  "05_CMP.LMPP"    "17_B"           "19_CD8.N"      
#[13] "09_pDC"         "13_CD16.Mono"   "23_CD8.EM"      "12_CD14.Mono.2"
#[17] "20_CD4.N1"      "02_Early.Eryth" "21_CD4.N2"      "24_CD8.CM"     
#[21] "01_HSC"

从上面的输出中，我们发现scRNA-seq数据中与NK细胞和T细胞对应的聚类是Cluster 19 - 25。

接着我们创建一个字符串模式用来表示这些聚类，后续的约束整合会用到，其中|在正则表达式中表示”或”，我们之后使用grep根据这些字符串模式从scATAC-seq提取和scRNA-seq对应的聚类。。

1
2
3
4
#From scRNA
cTNK <- paste0(paste0(19:25), collapse="|")
cTNK
# [1] "19|20|21|22|23|24|25"

其余的聚类就称之为”Non-T cell, Non-NK cell”(例如Cluster 1 - 18)，也创建了对应的字符串模式

1
2
3
cNonTNK <- paste0(c(paste0("0", 1:9), 10:13, 15:18), collapse="|")
cNonTNK
#[1] "01|02|03|04|05|06|07|08|09|10|11|12|13|15|16|17|18"

接着再用字符串模式在preClust找到对应的scATAC-seq列名，然后使用列名从混合矩阵提取对应的列。

对于T细胞和NK细胞，scATAC-seq聚类ID就是C7, C8, C9

1
2
3
4
#Assign scATAC to these categories
clustTNK <- rownames(cM)[grep(cTNK, preClust)]
clustTNK
#[1] "C8" "C9" "C7"

对于” Non-T cells and Non-NK cells”, ID就是scATAC-seq聚类余下的部分

1
2
3
clustNonTNK <- rownames(cM)[grep(cNonTNK, preClust)]
clustNonTNK
# [1] "C10" "C3"  "C11" "C4"  "C1"  "C12" "C5"  "C2"  "C6"

接着在scRNA-seq中做相同的操作，筛选出相同的细胞类型。首先，我们鉴定scRNA-seq数据中T细胞和NK细胞

1
2
3
4
5
6
7
8
9
#RNA get cells in these categories
rnaTNK <- colnames(seRNA)[grep(cTNK, colData(seRNA)$BioClassification)]
head(rnaTNK)
#[1] "PBMC_10x_GREENLEAF_REP1:AAACCCAGTCGTCATA-1"
#[2] "PBMC_10x_GREENLEAF_REP1:AAACCCATCCGATGTA-1"
#[3] "PBMC_10x_GREENLEAF_REP1:AAACCCATCTCAACGA-1"
#[4] "PBMC_10x_GREENLEAF_REP1:AAACCCATCTCTCGAC-1"
#[5] "PBMC_10x_GREENLEAF_REP1:AAACGAACAATCGTCA-1"
#[6] "PBMC_10x_GREENLEAF_REP1:AAACGAACACGATTCA-1"

然后，鉴定scRNA-seq数据中”Non-T cell Non-NK cell cells”

1
2
3
4
5
rnaNonTNK <- colnames(seRNA)[grep(cNonTNK, colData(seRNA)$BioClassification)]
head(rnaNonTNK)
#[1] "CD34_32_R5:AAACCTGAGTATCGAA-1" "CD34_32_R5:AAACCTGAGTCGTTTG-1"
#[3] "CD34_32_R5:AAACCTGGTTCCACAA-1" "CD34_32_R5:AAACGGGAGCTTCGCG-1"
#[5] "CD34_32_R5:AAACGGGAGGGAGTAA-1" "CD34_32_R5:AAACGGGAGTTACGGG-1"

约束整合需要我们提供一个嵌套list。这是一个由多个SimpleList对象组成的SimpleList, 每一组对应一个约束情况。在案例中，我们有两个组，一个组称之为TNK,包括两个平台的T和NK细胞，另一个组为NonTNK，包括两个平台的”Non-T cell Non-NK cell”细胞。每个SimpleList对象都包含两个细胞ID的向量，一个是ATAC，一个是RNA.

1
2
3
4
5
6
7
8
9
10
groupList <- SimpleList(
    TNK = SimpleList(
        ATAC = projHeme2$cellNames[projHeme2$Clusters %in% clustTNK],
        RNA = rnaTNK
    ),
    NonTNK = SimpleList(
        ATAC = projHeme2$cellNames[projHeme2$Clusters %in% clustNonTNK],
        RNA = rnaNonTNK
    )    
)

我们将该列表传递给addGeneIntegrationMatrix()函数的groupList参数。注意，我们依旧没有将结果添加到Arrow文件中 (addToArrow = FALSE)。我们强烈建议，在保存到Arrow文件前先彻底的检查结果，看结果是否符合预期。在教程的下一节会介绍该操作。

1
2
3
4
5
6
7
8
9
10
11
12
13
projHeme2 <- addGeneIntegrationMatrix(
    ArchRProj = projHeme2, 
    useMatrix = "GeneScoreMatrix",
    matrixName = "GeneIntegrationMatrix",
    reducedDims = "IterativeLSI",
    seRNA = seRNA,
    addToArrow = FALSE, 
    groupList = groupList,
    groupRNA = "BioClassification",
    nameCell = "predictedCell_Co",
    nameGroup = "predictedGroup_Co",
    nameScore = "predictedScore_Co"
)

8.1.3 约束整合和无约束整合对比

正如之前所提到的，我们的scATAC-seq数据可以根据整合的scRNA-seq数据进行定义，并且有约束和无约束这两种方式。为了对两者进行对比，我们分别根据这两种整合结果对scATAC-seq的数据进行上色。

首先，使用ArchR内置的paletteDiscrete()函数生成调色板

1
pal <- paletteDiscrete(values = colData(seRNA)$BioClassification)

在ArchR中，调色板本质上一个命名向量，每个十六进制编码的颜色对应着一个名字。

1
2
3
4
5
6
7
8
9
10
11
12
13
pal
#        01_HSC 02_Early.Eryth  03_Late.Eryth  04_Early.Baso    05_CMP.LMPP 
#     "#D51F26"      "#502A59"      "#235D55"      "#3D6E57"      "#8D2B8B" 
#      06_CLP.1         07_GMP    08_GMP.Neut         09_pDC         10_cDC 
#     "#DE6C3E"      "#F9B712"      "#D8CE42"      "#8E9ACD"      "#B774B1" 
#11_CD14.Mono.1 12_CD14.Mono.2   13_CD16.Mono         14_Unk       15_CLP.2 
#     "#D69FC8"      "#C7C8DE"      "#8FD3D4"      "#89C86E"      "#CC9672" 
#      16_Pre.B           17_B      18_Plasma       19_CD8.N      20_CD4.N1 
#     "#CF7E96"      "#A27AA4"      "#CD4F32"      "#6B977E"      "#518AA3" 
#     21_CD4.N2       22_CD4.M      23_CD8.EM      24_CD8.CM          25_NK 
#     "#5A5297"      "#0F707D"      "#5E2E32"      "#A95A3C"      "#B28D5C" 
#        26_Unk 
#     "#3D3D3D"

我们在scATAC-seq数据根据无约束整合得到的scRNA-seq细胞类型进行可视化

1
2
3
4
5
6
7
p1 <- plotEmbedding(
    projHeme2, 
    colorBy = "cellColData", 
    name = "predictedGroup_Un", 
    pal = pal
)
p1
Plot-UMAP-RNA-Integration_1 
1
do.call(cowplot::plot_grid, c(list(ncol = 3), p2c))

同样，我们也可以根据约束整合得到scATAC-seq对应的scRNA-seq的细胞类型进行可视化

1
2
3
4
5
6
7
p2 <- plotEmbedding(
    projHeme2, 
    colorBy = "cellColData", 
    name = "predictedGroup_Co", 
    pal = pal
)
p2
Plot-UMAP-RNA-Integration_2

这两者的结果差异其实非常细微，主要是我们感兴趣的细胞类型原本就存在明显的差异。当然，仔细观察还能发现其中的不同之处，尤其是T细胞(Clusters 17-22)

我们用plotPDF()函数保存该图矢量版本。

1
plotPDF(p1,p2, name = "Plot-UMAP-RNA-Integration.pdf", ArchRProj = projHeme2, addDOC = FALSE, width = 5, height = 5)

我们现在可以用saveArchRProject()函数保存我们的projHeme2对象。

1
saveArchRProject(ArchRProj = projHeme2, outputDirectory = "Save-ProjHeme2", load = FALSE)

8.2 为每个scATAC-seq细胞增加拟scRNA-seq谱

既然我们对scATAC-seq和scRNA-seq整合的结果感到满意，我们就能用addToArrow=TRUE重新运行，在Arrow文件中添加相关联的表达量矩阵数据。根据之前所提到的，我们传入groupList约束整合，在nameCell，nameGroup和nameScore中加入列名。这些元信息列都会被添加到cellColData中。

这里，我们新建了一个projHeme3，用于后续教程。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#~5 minutes
projHeme3 <- addGeneIntegrationMatrix(
    ArchRProj = projHeme2, 
    useMatrix = "GeneScoreMatrix",
    matrixName = "GeneIntegrationMatrix",
    reducedDims = "IterativeLSI",
    seRNA = seRNA,
    addToArrow = TRUE,
    force= TRUE,
    groupList = groupList,
    groupRNA = "BioClassification",
    nameCell = "predictedCell",
    nameGroup = "predictedGroup",
    nameScore = "predictedScore"
)

现在，当我们使用getAvailableMatrices()检查哪些矩阵可用时，我们会发现GeneIntegrationMatrix已经被添加到Arrow文件中

1
2
getAvailableMatrices(projHeme3)
# [1] "GeneIntegrationMatrix" "GeneScoreMatrix" "TileMatrix"

在新的GeneIntegrationMatrix中，我们可以比较连接的基因表达量和根据基因得分推断的基因表达量

我们需要先确保在项目中加入了填充权重值(impute weights)

1
projHeme3 <- addImputeWeights(projHeme3)

现在，我们来生成一些UMAP图，里面的基因表达量值来自于GeneIntegrationMatrix

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
markerGenes  <- c(
    "CD34", #Early Progenitor
    "GATA1", #Erythroid
    "PAX5", "MS4A1", #B-Cell Trajectory
    "CD14", #Monocytes
    "CD3D", "CD8A", "TBX21", "IL7R" #TCells
  )

p1 <- plotEmbedding(
    ArchRProj = projHeme3, 
    colorBy = "GeneIntegrationMatrix", 
    name = markerGenes, 
    continuousSet = "horizonExtra",
    embedding = "UMAP",
    imputeWeights = getImputeWeights(projHeme3)
)

以及一些相同UMAP图，但是使用GeneScoreMatrix里的基因得分值

1
2
3
4
5
6
7
8
p2 <- plotEmbedding(
    ArchRProj = projHeme3, 
    colorBy = "GeneScoreMatrix", 
    continuousSet = "horizonExtra",
    name = markerGenes, 
    embedding = "UMAP",
    imputeWeights = getImputeWeights(projHeme3)
)

最后用cowplot将这些标记基因绘制在一起

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
p1c <- lapply(p1, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})

p2c <- lapply(p2, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})

do.call(cowplot::plot_grid, c(list(ncol = 3), p1c))
Plot-UMAP-Markers-RNA-W-Imputation_1 
1
do.call(cowplot::plot_grid, c(list(ncol = 3), p2c))
Plot-UMAP-Markers-RNA-W-Imputation_2

和预期的一样，两个方法推测的基因表达量存在相似性，但并不相同。

使用plotPDF()函数保存可编辑的矢量版。

1
2
3
4
plotPDF(plotList = p1, 
    name = "Plot-UMAP-Marker-Genes-RNA-W-Imputation.pdf", 
    ArchRProj = projHeme3, 
    addDOC = FALSE, width = 5, height = 5)

8.3 使用scRNA-seq信息标记scATAC-seq聚类

现在，我们确定了scATAC-seq和scRNA-seq数据间的对应关系，我们就可以使用scRNA-seq数据中细胞类型对我们的scATAC-seq聚类进行定义。

首先，我们会在scATAC-seq和整合分析得到predictedGroup之间构建一个混合矩阵

1
2
3
4
5
6
cM <- confusionMatrix(projHeme3$Clusters, projHeme3$predictedGroup)
labelOld <- rownames(cM)
labelOld
#  [1] "Cluster11" "Cluster2"  "Cluster12" "Cluster1"  "Cluster8"  "Cluster4" 
#  [7] "Cluster9"  "Cluster5"  "Cluster7"  "Cluster14" "Cluster3"  "Cluster10"
# [13] "Cluster6"  "Cluster13"

接着，对于每一个scATAC-seq聚类，我们根据predictedGroup确定最能定义聚类的细胞类型。

1
2
labelNew <- colnames(cM)[apply(cM, 1, which.max)]
labelNew

接着，我们需要对新的聚类标签进行重命名，简化分类系统。对于每一个scRNA-seq的聚类，我们会重新定义标签，以便更好解释。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
remapClust <- c(
    "01_HSC" = "Progenitor",
    "02_Early.Eryth" = "Erythroid",
    "03_Late.Eryth" = "Erythroid",
    "04_Early.Baso" = "Basophil",
    "05_CMP.LMPP" = "Progenitor",
    "06_CLP.1" = "CLP",
    "07_GMP" = "GMP",
    "08_GMP.Neut" = "GMP",
    "09_pDC" = "pDC",
    "10_cDC" = "cDC",
    "11_CD14.Mono.1" = "Mono",
    "12_CD14.Mono.2" = "Mono",
    "13_CD16.Mono" = "Mono",
    "15_CLP.2" = "CLP",
    "16_Pre.B" = "PreB",
    "17_B" = "B",
    "18_Plasma" = "Plasma",
    "19_CD8.N" = "CD8.N",
    "20_CD4.N1" = "CD4.N",
    "21_CD4.N2" = "CD4.N",
    "22_CD4.M" = "CD4.M",
    "23_CD8.EM" = "CD8.EM",
    "24_CD8.CM" = "CD8.CM",
    "25_NK" = "NK"
)
remapClust <- remapClust[names(remapClust) %in% labelNew]

接着，使用mapLables()函数进行标签转换，将旧的标签映射到新的标签上。

1
2
3
4
5
labelNew2 <- mapLabels(labelNew, oldLabels = names(remapClust), newLabels = remapClust)
labelNew2
#  [1] "GMP"        "B"          "PreB"       "CD4.N"      "Mono"      
#  [6] "Erythroid"  "Progenitor" "CD4.M"      "pDC"        "NK"        
# [11] "CLP"        "Mono"

合并labelsOld和labelsNew2，我们现在可以用mapLables()函数在cellColData里新建聚类标签。

1
projHeme3$Clusters2 <- mapLabels(projHeme3$Clusters, newLabels = labelNew2, oldLabels = labelOld)

得到新的标签后，我们绘制最新的UMAP展示聚类结果。

1
2
p1 <- plotEmbedding(projHeme3, colorBy = "cellColData", name = "Clusters2")
p1

如果被分析scATAC-seq数据对应的细胞系统也有scRNA-seq数据存在，那么这种分析范式将会特别有用。正如之前所说，这种scRNA-seq和scATAC-seq整合分析也为后续更加复杂的基因调控分析提供了出色的框架。

用plotPDF()函数保存该图矢量版本。

1
plotPDF(p1, name = "Plot-UMAP-Remap-Clusters.pdf", ArchRProj = projHeme2, addDOC = FALSE, width = 5, height = 5)

使用saveArchRProject保存我们最初的projHeme3.

1
saveArchRProject(ArchRProj = projHeme3, outputDirectory = "Save-ProjHeme3", load = FALSE)
R

使用ArchR分析单细胞ATAC-seq数据(第七章)

第7章: ArchR的基因得分和标记基因

尽管ArchR能够可靠地对细胞进行聚类，但它并不能先验(a priori)地知道每个聚类所代表的细胞类型。通常这个任务也只能靠人工注释，毕竟每个项目都不太一样。

为了实现细胞类型注释，我们会用到已知的细胞类型特异的标记基因。在scATAC-seq数据中，这些基因的表达量其实是根据染色质开放数据估计的基因得分(gene score)。所谓的基因得分，本质上就是用基因附近的调控元件去预测基因的表达量。ArchR的亮点在于，它允许用户提供复杂的自定义的距离加权开放性模型去计算这些基因得分。

7.1 在ArchR中计算基因得分

我们在文章中测试了超过50种不同的基因得分模型，并找到一种在大部分测试情况下都拥有良好表现的模型。这个模型目前是ArchR的默认设置，它包括三个主要部分

  1. 在整个基因主体内(gene body)的开放部分都对基因得分有贡献
  2. 对于潜在的远距调控元件，我们使用一个指数加权函数计算得分
  3. 我们加入了基因边界(gene boundary)用于减少无关调控元件对基因得分的贡献

那么，ArchR到底是怎么计算基因得分的呢？我建议先看图，再看文字描述。图中标注了两个TSS，分别表示两个基因的转录起始位点(TSS)，其中红色标识的基因是目标基因，我们需要根据其附近的开放情况来预测它的表达量。

基因活跃得分

对于每条染色体，ArchR都会构建一个分块矩阵(tile matrix)，块的大小由用户进行定义，在定义时计算（默认是500 bp），然后分析用户定义的另一个基因分窗(默认是基因两边100kb)和上一个分块矩阵的重叠部分。接着计算每个分块的起始位置或结束位置相对于基因主体(可额外添加上下游延伸)或基因起始位置的距离。我们发现基因表达量最好的预测者是基因区的局部开放性，包括启动子和基因主体。对于之前提到的远距调控元件，为了能正确地处理给定基因的远距离开放状态，ArchR只筛选了在基因分窗内但不跨到另一个基因区的分块。这个过滤标准既能利用远端调控元件提高基因表达量的预测准确性，同时也避免了无关开放状态的影响（例如，邻近基因启动子的开放状态）。每个分块到基因的距离会根据用户定义的开放模型里的权重进行转换，默认是e(-abs(distance)/5000) + e-1)。

我们还发现基因长度会明显影响总体的基因得分，这是因为如果基因区包括基因主体，当基因越长时，落在基因主体内的开放区也就越多，这些开放区会和最大的基于距离的权重相乘，也就会计算出更高的得分。为了调整基因长度所带来的差异，ArchR以每个基因长度的倒数(1 / gene size)作为权重，然后将该权重线性缩放到1到用户定义的最大值之间(默认是5)。保证了越短的基因对应更大的相对权重，一定程度上避免了长度的影响。

之后每个分块对应的距离权重和基因长度权重会和每个分块内的Tn5 insertion数相乘，接着将基因分窗内的所有分块的得分进行相加，相加时排除邻近基因区的分块。最后计算的总开放度就是基因得分(gene score)，接着所有基因的基因得分会按照用户定义的常数进行深度标准化(默认是10,000)。最后，计算的基因得分会被保存在相应的Arrow文件中用作下游分析。

在创建Arrow文件时，如果设置createArrowFiles的参数addGeneScoreMat = TRUE(默认为TRUE)，那么生成的Arrow文件中会包含基因得分矩阵。此外，我们也可以在任意时候使用addGeneScoreMatrix()在Arrow文件中加入基因得分矩阵。一旦计算完成，每个细胞都可以在嵌入图中根据基因得分进行展示，方便鉴定不同的细胞类型。之后的章节中会逐一介绍基因得分的应用。

当然并不是所有基因都会有基因得分，尤其是基因密度高的区域里的基因，它们更容易出现问题。因此，最好能在基因组浏览器对所有基因得分都做一遍检查，我们也会在后续章节里做相应介绍。

7.2 鉴定标记特征

除了使用已知的细胞类型相关标记基因用于聚类注释外，ArchR还能无偏地鉴定任何给定细胞分组(如聚类)的标记特征。这些特征包括单不限于peaks, 基于基因得分的基因, 基于chromVAR离差的转录因子motif。这部分工作可通过getMarkerFeatures()函数完成，它的核心参数是useMatrix和groupBy。useMatrix可以是任何的矩阵输入，如果设置为GeneScoreMatrix那么函数就会分析不同细胞类型特异出现的基因。这就能无偏的找到每个聚类的活跃基因，可以用来辅助聚类注释。

正如之前所提到的，getMarkerFeatures()能够接受任何存在Arrow文件里的矩阵作为输入，分析细胞类型的特异特征。例如useMatrix = "TileMatrix"用于鉴定细胞类群的特异基因区间，useMatrix = "PeakMatrix"用于鉴定细胞类群的特异peak。后续章节还会介绍如何使用其他特征类型作为getMarkerFeatures()的输入。

7.2.1 标记特征的鉴定过程

标记特征的鉴定过程取决于每一组细胞对应的偏好-匹配背景细胞的选择。对于所有特征，每个细胞组都会与它自己的背景细胞组进行比较，判断给定的细胞组有更显著的开放性。

markerFeature_schematic

选择合适的背景细胞组对标记特征的鉴定至关重要，取决于getMarkerFeatures)()的bias参数所构建的多维空间。对于给定分组的细胞，ArchR会在提供的多维空间中寻找每个细胞分组外的最近邻细胞。这些细胞和给定分组内细胞极其相似，因此被称为偏好-匹配细胞(bias-matched cells)。通过这一方法，即便是在细胞数目比较少的组中，我们都能为每个特征计算稳健的显著性。

ArchR以bias参数中的所有维度作为输入，以分位数对值进行归一化，让每个维度的方差分布在相同的相对标度。举一个简单的例子，如如果bias的输入是TSS和log10(Num Fragments)，未经分位数归一化的值如下图所示

ackground_preNorm

这里y轴的相对方差和于x轴的相对方差相比，就显得特别小。如果我们对这些轴进行归一化，将值缩放到0-1之间，x和y的相对方差就近乎一样。注意，这个操作还会明显地改变每个细胞的最近邻（见下图右）。

background_postNorm

ArchR会对所有维度进行归一化，在归一化后的多维空间中以欧氏距离寻找最近邻。

7.3 鉴定标记基因

为了根据基因得分来鉴定标记基因，我们需要在调用getMarkerFeatures()时设置useMatrix = "GeneScoreMatrix。此外，我们设置了groupBy="Clusters"告诉ArchR使用cellColData的”Clusters”列对细胞分组，从而得到每个聚类特异的标记基因。

1
2
3
4
5
6
7
markersGS <- getMarkerFeatures(
    ArchRProj = projHeme2, 
    useMatrix = "GeneScoreMatrix", 
    groupBy = "Clusters",
    bias = c("TSSEnrichment", "log10(nFrags)"),
    testMethod = "wilcoxon"
)

该函数返回一个SummarizedExperiment对象，里面记录着每个标记特征的相关信息。这是ArchR的一个常见输出，是下游数据分析的关键输出格式。SummarizedExperiment对象和矩阵类似，行是感兴趣的特征（例如基因），列表示样本。一个SummarizedExperiment能够记录一个或多个assay矩阵。如果你需要了解更多和SummarizedExperiment相关的内容，建议阅读相关的Bioconductor页面，这里就不深入对这方面内容进行介绍。

我们能够使用getMarkers()函数从SummarizedExperiment对象中提取出一个包含多个DataFrame对象的列表，每个DataFrame对象都对应着一个聚类，记录着该聚类里相关的标记特征。

1
2
markerList <- getMarkers(markersGS, cutOff = "FDR <= 0.01 & Log2FC >= 1.25")
markerList$C6

我们使用markerHeatmap()创建含有所有标记特征的可视化热图。通过labelMarkers()在热图中标记部分感兴趣的基因。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
markerGenes  <- c(
    "CD34", #Early Progenitor
    "GATA1", #Erythroid
    "PAX5", "MS4A1", "EBF1", "MME", #B-Cell Trajectory
    "CD14", "CEBPB", "MPO", #Monocytes
    "IRF8", 
    "CD3D", "CD8A", "TBX21", "IL7R" #TCells
  )

heatmapGS <- markerHeatmap(
  seMarker = markersGS, 
  cutOff = "FDR <= 0.01 & Log2FC >= 1.25", 
  labelMarkers = markerGenes,
  transpose = TRUE
)

我们需要用ComplexHeatmap::draw()才能绘制该热图，这是因为heatmapGS实际上是一个热图列表。

1
ComplexHeatmap::draw(heatmapGS, heatmap_legend_side = "bot", annotation_legend_side = "bot")
GeneScores Marker heatmap

使用plotPDF()函数保存可编辑的矢量版。

1
plotPDF(heatmapGS, name = "GeneScores-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme2, addDOC = FALSE)

7.4 在嵌入上可视化标记基因

在之前章节中提到，我们可以在UMAP嵌入上展示每个细胞的基因得分。这可以通过修改plotEmbedding()函数里的colorBy和name参数实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
markerGenes  <- c(
    "CD34",  #Early Progenitor
    "GATA1", #Erythroid
    "PAX5", "MS4A1", "MME", #B-Cell Trajectory
    "CD14", "MPO", #Monocytes
    "CD3D", "CD8A"#TCells
  )

p <- plotEmbedding(
    ArchRProj = projHeme2, 
    colorBy = "GeneScoreMatrix", 
    name = markerGenes, 
    embedding = "UMAP",
    quantCut = c(0.01, 0.95),
    imputeWeights = NULL
)

选择绘图列表的其中一个基因进行展示

1
p$CD14
Plot-UMAP-Marker-Genes-WO-Imputation

如果是需要绘制所有基因，那么可以使用cowplot将不同的标记基因整合到一张图中

1
2
3
4
5
6
7
8
9
10
11
12
p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
do.call(cowplot::plot_grid, c(list(ncol = 3),p2))
UMAP of Marker gene

使用plotPDF()函数保存可编辑的矢量版。

1
2
3
4
plotPDF(plotList = p, 
    name = "Plot-UMAP-Marker-Genes-WO-Imputation.pdf", 
    ArchRProj = projHeme2, 
    addDOC = FALSE, width = 5, height = 5)

7.5 使用MAGIC填充标记基因

在上一节中，你可能注意到一些基因得分图变化很大，这是因为scATAC-seq数据太过稀疏。我们使用MAGIC根据邻近细胞填充基因得分对信号进行平滑化处理。我们发现者能够极大程度地提高基因得分的可视化解读性。要执行此操作，我们需要先在我们的ArchRProject中加入填充权重。

1
projHeme2 <- addImputeWeights(projHeme2)

这些填充权重会在之后绘制UMAP嵌入图里的基因得分时传入到plotEmbedding()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
markerGenes  <- c(
    "CD34",  #Early Progenitor
    "GATA1", #Erythroid
    "PAX5", "MS4A1", "MME", #B-Cell Trajectory
    "CD14", "MPO", #Monocytes
    "CD3D", "CD8A"#TCells
  )

p <- plotEmbedding(
    ArchRProj = projHeme2, 
    colorBy = "GeneScoreMatrix", 
    name = markerGenes, 
    embedding = "UMAP",
    imputeWeights = getImputeWeights(projHeme2)
)

和之前一样，我们可以只画其中一个基因

UMAP of imputation

也可以用cowplot绘制所有的标记基因

1
2
3
4
5
6
7
8
9
10
11
12
13
#Rearrange for grid plotting
p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
do.call(cowplot::plot_grid, c(list(ncol = 3),p2))
UMAP of all imputation

使用plotPDF()函数保存可编辑的矢量版。

1
2
3
4
plotPDF(plotList = p, 
    name = "Plot-UMAP-Marker-Genes-W-Imputation.pdf", 
    ArchRProj = projHeme2, 
    addDOC = FALSE, width = 5, height = 5)

7.6 使用ArchRBrowser绘制Track

除了在UMAP上绘制每个细胞的基因得分外，我们还可以在基因组浏览器上浏览这些标记基因的局部染色体开放状态。为了执行该操作，我们使用plotBrowserTrack()函数，它会构建一系列绘图对象，每个对象对应一个标记基因。函数会根据groupBy输入的组别信息在不同track上绘制每组的开放状态。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
markerGenes  <- c(
    "CD34", #Early Progenitor
    "GATA1", #Erythroid
    "PAX5", "MS4A1", #B-Cell Trajectory
    "CD14", #Monocytes
    "CD3D", "CD8A", "TBX21", "IL7R" #TCells
  )

p <- plotBrowserTrack(
    ArchRProj = projHeme2, 
    groupBy = "Clusters", 
    geneSymbol = markerGenes, 
    upstream = 50000,
    downstream = 50000
)

通过选择列表中的给定基因来绘制最终结果

1
2
grid::grid.newpage()
grid::grid.draw(p$CD14)
Plot-Tracks-Marker-Genes_5

使用plotPDF()函数可以将多幅基因座位对应的图保存在一个PDF文件中。

1
2
3
4
plotPDF(plotList = p, 
    name = "Plot-Tracks-Marker-Genes.pdf", 
    ArchRProj = projHeme2, 
    addDOC = FALSE, width = 5, height = 5)

7.7 启动ArchRBrowser

scATAC-seq数据分析一个与生俱来的挑战就是在基因组浏览器上可视化聚类间的染色质开放水平。传统做法是，先对scATAC-seq的fragments进行分组，然后构建一个基因组覆盖度的bigwig文件，最后对track进行标准化后才能实现定量层面上的可视化。终端用户经常用于测序数据可视化的基因组浏览器包括 WashU Epigenome Browser, UCSC Genome Browser, IGV browser。这一过程需要用到多个不同软件，并且一旦细胞分组信息改变，或者增加了更多的样本，就得重新分组输出bigiwig文件，相当耗时。

为了避免重复劳动，ArchR开发了基于Shiny的交互式基因组浏览器，只需要一行ArchRBrowser(ArchRProj)命令就可启动。因为数据都存放在Arrow文件中，所以这个交互式浏览器就能动态的更改细胞分组信息，分辨率以及标准化，实现实时的track水平的可视化。ArchR Genome Browser 同样也能生成高质量的PDF格式输出文件用于发表或分享。此外，浏览器支持用户通过features参数传入GenomicRanges对象用来展示特征，或者通过loop参数传入基因组交互文件（ co-accessibility, peak-to-gene linkages, loops from chromatin conformation data）。其中loop的预期格式是GRanges，起点表示其中一个loop的中心，终点表示另一个loop中心。

使用ArchRBrowser()函数启动我们的本地交互基因组浏览器

ArchR_Browser_1

通过选择”Gene Symbol”，就可以开始浏览了。你可能需要点击”Plot Track”才能强制让你的浏览器刷新

ArchR_Browser_2

当我们绘制一个基因位点后，我们看到不同track代表的是我们数据的不同聚类。

ArchR_Browser_3

如果我们点击”Additional Parameters”侧边栏，我们可以选择部分聚类进行展示

ArchR_Browser_4

通过反选cluster1, 2, 3，我们就在绘图中移除了它们

ArchR_Browser_5

当我们返回到”Plot”时，此时Cluster1,2,3就从图片中消失了。同样你也可能需要点击”Plot Track”才能让浏览器刷新结果

ArchR_Browser_6

无论在哪个阶段，我们都可以点击”Download The Track”来输出当前的绘图结果。

R

使用ArchR分析单细胞ATAC-seq数据(第六章)

第6章: 单细胞嵌入

在ArchR中，类似于UMAP和t-SNE的嵌入方法被用于在降维空间中可视化展示单细胞数据。这些嵌入有各自的优势和缺陷。我们之所以称它们为”嵌入”是因为他们只限于对聚类进行可视化而非用于鉴定聚类(在LSI子空间中的聚类分析)。UMAP和t-SNE的主要区别在于对细胞和聚类间的距离解释，t-SNE用于保留数据的局部结构，而UMAP则是保留局部结构的同时尽可能保留全局结构。从理论上来讲，UMAP中细胞聚类间的距离存在意义，而t-SNE则没有。也就是说，你不能说t-SNE中聚类A比聚类C更接近聚类B，而UMAP在设计的时候允许这种类型的比较。不过需要注意的是，由于UMAP是新出现的方法，因此使用UMAP的文章才刚刚兴起。

需要注意的是，无论是UMAP还是t-SNE，两个的运行结果都不是确定性的，也就是相同输入会得到不完全相同的结果。尽管如此，我们使用样本重复后发现t-SNE比UMAP更加的随机。此外，使用uwot包里UMAP时，设置相同的随机数种子会得到相同的结果。选择使用UMAP还是t-SNE是有细微差别的，但在我们手中，UMAP非常适合各种应用，这是我们对scATAC seq数据的标准选择。UMAP的性能也比t-SNE快。也许最重要的是，使用UMAP可以创建一个嵌入并将新样本投影到该嵌入中，而这在t-SNE中是不可能的，因为数据的拟合和预测是同时发生的。

无论您选择哪种方法，输入参数都会对结果嵌入产生剧烈影响。因此，了解各种输入参数并调整这些参数以最好地满足您自己的数据需要是很重要的。ArchR实现了一组默认的输入参数，这些参数适用于大多数情况，但实际上没有一组参数可以直接套用，我们要根据不同的细胞数、复杂度和质量进行调整。

6.1 Uniform Manifold Approximation and Projection (UMAP)

我们使用addUMAP函数运行UMAP

1
2
3
4
5
6
7
8
projHeme2 <- addUMAP(
    ArchRProj = projHeme2, 
    reducedDims = "IterativeLSI", 
    name = "UMAP", 
    nNeighbors = 30, 
    minDist = 0.5, 
    metric = "cosine"
)

使用@操作符可以从ArchRProject中列出有哪些可用的embedding，如projHeme2@embeddings

我们使用plotEmbedding函数绘制UMAP图，设置embedding="UMAP"。通过修改colorBy和name来告诉ArchR使用给定哪个元信息矩阵的列对细胞进行上色。p1是按照样本进行上色

1
p1 <- plotEmbedding(ArchRProj = projHeme2, colorBy = "cellColData", name = "Sample", embedding = "UMAP")

除了根据”Sample”外，我们还可以根据上一张鉴定的”Cluster”进行上色

1
p2 <- plotEmbedding(ArchRProj = projHeme2, colorBy = "cellColData", name = "Clusters", embedding = "UMAP")

使用ggAlignPlots将这两个图共同展示，”type=h”表示两个图是水平排列

1
ggAlignPlots(p1, p2, type = "h")

UMAP of Seurat

用plotPDF()可以将保存图片的矢量版。

1
plotPDF(p1,p2, name = "Plot-UMAP-Sample-Clusters.pdf", ArchRProj = projHeme2, addDOC = FALSE, width = 5, height = 5)

我们还可以使用plotEmbedding()可视化之前用scran聚类的结果

1
2
3
p1 <- plotEmbedding(ArchRProj = projHeme2, colorBy = "cellColData", name = "Sample", embedding = "UMAP")
p2 <- plotEmbedding(ArchRProj = projHeme2, colorBy = "cellColData", name = "ScranClusters", embedding = "UMAP")
ggAlignPlots(p1, p2, type = "h")

UMAP of scran

同样用plotPDF()可以将保存图片的矢量版。

1
plotPDF(p1,p2, name = "Plot-UMAP-Sample-ScranClusters.pdf", ArchRProj = projHeme2, addDOC = FALSE, width = 5, height = 5)

6.2 t-Stocastic Neighbor Embedding (t-SNE)

t-SNE图可以用addTSNE()函数运行得到

1
2
3
4
5
6
projHeme2 <- addTSNE(
    ArchRProj = projHeme2, 
    reducedDims = "IterativeLSI", 
    name = "TSNE", 
    perplexity = 30
)

和之前UMAP一样，我们同样使用plotEmbedding()绘制t-SNE嵌入图。这里不需要考虑嵌入的类型，可以继续使用之前的colorBy和name参数

1
2
3
p1 <- plotEmbedding(ArchRProj = projHeme2, colorBy = "cellColData", name = "Sample", embedding = "TSNE")
p2 <- plotEmbedding(ArchRProj = projHeme2, colorBy = "cellColData", name = "Clusters", embedding = "TSNE")
ggAlignPlots(p1, p2, type = "h")

t-SNE of Seurat

用plotPDF()可以将保存图片的矢量版。

1
plotPDF(p1,p2, name = "Plot-TSNE-Sample-Clusters.pdf", ArchRProj = projHeme2, addDOC = FALSE, width = 5, height = 5)

和之前UMAP的操作类似，我们可以将scran的结果和Seurat::FindClusters()的结果进行比较

1
2
3
p1 <- plotEmbedding(ArchRProj = projHeme2, colorBy = "cellColData", name = "Sample", embedding = "TSNE")
p2 <- plotEmbedding(ArchRProj = projHeme2, colorBy = "cellColData", name = "ScranClusters", embedding = "TSNE")
ggAlignPlots(p1, p2, type = "h")

t-SNE of Scran

用plotPDF()可以将保存图片的矢量版。

1
plotPDF(p1,p2, name = "Plot-tSNE-Sample-ScranClusters.pdf", ArchRProj = projHeme2, addDOC = FALSE, width = 5, height = 5)

6.3 Harmony后降维

在第4章中，我们通过addHarmony调用Harmony对数据进行批次校正，创建了一个名为”Harmony”的reducedDims对象。我们通过UMAP或t-SNE对结果进行嵌入可视化，对迭代LSI结果和Harmony校正结果进行比较，评估Harmony的作用。

保持和之前UMAP嵌入一样的参数，只修改reducedDims="Harmony"

1
2
3
4
5
6
7
8
9
10
11
projHeme2 <- addUMAP(
    ArchRProj = projHeme2, 
    reducedDims = "Harmony", 
    name = "UMAPHarmony", 
    nNeighbors = 30, 
    minDist = 0.5, 
    metric = "cosine"
)
p3 <- plotEmbedding(ArchRProj = projHeme2, colorBy = "cellColData", name = "Sample", embedding = "UMAPHarmony")
p4 <- plotEmbedding(ArchRProj = projHeme2, colorBy = "cellColData", name = "Clusters", embedding = "UMAPHarmony")
ggAlignPlots(p3, p4, type = "h")

UMAP of Harmony

用plotPDF()可以将保存图片的矢量版。

1
plotPDF(p1,p2,p3,p4, name = "Plot-UMAP2Harmony-Sample-Clusters.pdf", ArchRProj = projHeme2, addDOC = FALSE, width = 5, height = 5)

相同的方法用t-SNE进行可视化

1
2
3
4
5
6
7
8
9
projHeme2 <- addTSNE(
    ArchRProj = projHeme2, 
    reducedDims = "Harmony", 
    name = "TSNEHarmony", 
    perplexity = 30
)
p3 <- plotEmbedding(ArchRProj = projHeme2, colorBy = "cellColData", name = "Sample", embedding = "TSNEHarmony")
p4 <- plotEmbedding(ArchRProj = projHeme2, colorBy = "cellColData", name = "Clusters", embedding = "TSNEHarmony")
ggAlignPlots(p3, p4, type = "h")

t-SNE of Harmony

用plotPDF()可以将保存图片的矢量版。

1
plotPDF(p1,p2,p3,p4, name = "Plot-TSNE2Harmony-Sample-Clusters.pdf", ArchRProj = projHeme2, addDOC = FALSE, width = 5, height = 5)
R

使用ArchR分析单细胞ATAC-seq数据(第五章)

第5章: 使用ArchR聚类

大部分单细胞聚类算法都在降维后空间中计算最近邻图，然后鉴定”社区”或者细胞聚类。这些方法效果表现都特别出色，已经是scRNA-seq的标准策略，所以ArchR直接使用了目前scRNA-seq包中最佳的聚类算法用来对scATAC-seq数据进行聚类。

5.1 使用Seurat的FindClusters()函数进行聚类

我们发现Seurat实现的图聚类方法表现最好，所以在ArchR中，addClusters()函数本质是上将额外的参数通过...传递给Seurat::FindClusters()函数，从而得到聚类结果。在分析中，我们发现Seurat::FindClusters()是一个确定性的聚类算法，也就是相同的输入总是能得到完全相同的输出。

1
2
3
4
5
6
7
8
9
10
11
projHeme2 <- addClusters(
    input = projHeme2,
    reducedDims = "IterativeLSI",
    method = "Seurat",
    name = "Clusters",
    resolution = 0.8
)
# 只展示部分信息
# Maximum modularity in 10 random starts: 0.8568
# Number of communities: 11
# Elapsed time: 1 seconds

我们可以使用$符号来获取聚类信息，输出结果是每个细胞对应的cluster

1
2
head(projHeme2$Clusters)
# [1] "C10" "C6"  "C1"  "C2"  "C2"  "C10"

我们统计下每个cluster的细胞数

1
2
3
table(projHeme2$Clusters)
#  C1  C10  C11   C2   C3   C4   C5   C6   C7   C8   C9 
# 310 1247 1436  480  323  379  852 1271  677 2550  726

为了更好了解样本在cluster的分布，我们可以使用confusionMatrix()函数通过每个样本创建一个聚类混淆矩阵(cluster confusion matrix)

1
2
cM <- confusionMatrix(paste0(projHeme2$Clusters), paste0(projHeme2$Sample))
cM

文字信息太多，这里直接用热图进行展示

1
2
3
4
5
6
7
8
library(pheatmap)
cM <- cM / Matrix::rowSums(cM)
p <- pheatmap::pheatmap(
    mat = as.matrix(cM), 
    color = paletteContinuous("whiteBlue"), 
    border_color = "black"
)
p

混淆矩阵

细胞有时在二维嵌入中的相对位置与所识别的聚类并不完全一致。更确切的说，单个聚类中的细胞可能出现在嵌入的多个不同区域中。在这些情况下，可以适当地调整聚类参数或嵌入参数，直到两者之间达成一致。

5.2 使用scran聚类

除了Seurat, ArchR还能够使用scran进行聚类分析，我们只需要修改addClusters()中的method参数即可。

1
2
3
4
5
6
7
projHeme2 <- addClusters(
    input = projHeme2,
    reducedDims = "IterativeLSI",
    method = "scran",
    name = "ScranClusters",
    k = 15
)
R

使用ArchR分析单细胞ATAC-seq数据(第四章)

第4章: ArchR降维分析

scATAC-seq数据的稀疏性让降维这一步充满了挑战性。在scATAC-seq数据中，一个特定的位置存在三种状态，分为不开放，一条链开放，两条链都开放。但即便是在高质量的scATAC-seq数据中，大部分开放区域都没有被转座酶切割，因此大部分座位都是0。此外，如果我们发现一个细胞(A)的一处peak区域有三个Tn5 insertions，而另外一个细胞(B)对应区域只有一个Tn5 insertions，由于数据的稀疏性，我们也无法断定A就比B更开放。鉴于此，许多分析框架都会把scATAC-seq数据矩阵二值化，也就是只有0和1两种情况。即便如此，因为转座酶切割的位置还是很少，所以二值矩阵大部分区域还是0。我们还需要注意的是，scATAC-seq数据中的0既有可能表示”不开放”，也可能是”没取样到”, 从生物学角度上出发，这是两种截然相反的推断。正因为数据大部分都是0，所以1比0更有信息量。scATAC-seq数据的稀疏性便是来源于这种低信息量。

如果你直接在稀疏矩阵上使用标准的降维分析，例如主成分分析(Principal Component Analysis, PCA)，然后使用前两个主成分进行绘图，你可能无法得到你想要的结果，这是由于数据的稀疏性会导致细胞间在0的位置有更高的相似性。为了解决这个问题，我们采用了分层的降维策略。

首先，我们会使用隐语义分析(Latent Semantic Indexing, LSI)进行降维。LSI方法最早用在自然语言分析中，根据字数来评估文档相似度。之所以在自然语言分析中发明该方法，是因为文档集中会有许多不同的词，每个词出现的次数也很少，最终的数据非常稀疏且充满噪音。Cusanovich et al. (Science 2015)第一次在scATAC-sq数分析中引入了LSI。在scATAC-seq场景下，不同样本就是不同的文档，不同的区域/peak就是不同的单词(注1)。计算步骤如下

  1. 对每个细胞计算根据深度标准化后的词频(Term frequency, TF)
  2. 根据逆文档频率(IDF, inverse document frequency)进行标准化，便于后续特征选择
  3. 最后对结果矩阵进行对数变换（也就是log(TF-IDF)）

接着，奇异值分解( singular value decomposition, SVD)技术能够分析出不同样本间更有价值的信息，从而降低了维度。上面这两步能够将原本成千上万维的稀疏矩阵的降维到数十个或者数百个维度。

最后，我们使用更加常见降维方法，例如UMAP( Uniform Manifold Approximation and Projection ) 或t-SNE( t-distributed stochastic neighbor embedding ) 对数据进行可视化。ArchR称这这些可视化方法为嵌入(embeddings, 注2)

  • 注1: 对LSI最直观理解就是它在自然语言分析的用途，也就是为每篇文档找到关键词，这个关键词必须要在当前文档中出现，同时最好别在其他文档里出现。在scATAC-seq的语境中，如果一个区域在所有细胞中都有insertion，那么就不特异，最好是只有几个细胞有insertion，而大部分细胞没有，这才算一个特异的insertion。
  • 注2: 我们可以认为道路是嵌入在三维空间中一维流形。我们用一维道路中的地址号码确定地址，而非三维空间的坐标。

4.1 ArchR的LSI实现

ArchR实现了多种LSI方法，并用不同测数据集对这些方法进行了测试。ArchR默认的LSI实现和Timothy Stuart在Signac引入的方法有关，也就是先将词频(TF, term frequency)根据常数(10,000)进行深度标准化，然后根据逆文档频率(IDF, inverse document frequency)进行标准化，最后对结果矩阵进行对数变换（也就是log(TF-IDF)）

LSI降维的其中一个关键输入是起始矩阵。到目前为止，scATAC-seq主要有两种策略计算矩阵，一是使用peak区域，二是对全基因组分块(genome-wide tiles)。由于在降维前，我们还没有聚类，也没有聚类特异的peak，所以不能使用peak区域作为LSI的输入。而且直接使用所有数据进行peak calling，会导致一些细胞特异性的peak被掩盖。更何况，当你在实验中增加新的样本，那么合并后的peak集就会发生改变，导致整个方法不太稳定。第二种策略通过对全基因组分块，保证了特征集(feature set)的一致性和无偏好性，缓解了第一种方法存在的问题。只不过，使用全基因组分块会得到一个非常大的细胞X特征的矩阵，因此大部分LSI实现都至少用5000 bp对基因组进行分块。计算量降低的同时也导致分辨率急剧下降，因为大部分开放区域只有几百个碱基的长度。

在Arrow文件的设计帮助下，即便是以500-bp作为基因组的分块，ArchR依旧能够快速的计算LSI。这就解决了分辨率的问题，使得在calling peak前就能进行聚类分析。使用500-bp进行分块的更大挑战是这会得到一个拥有数百万个特征的cell-by-tile矩阵。尽管ArchR能够通过读取必要矩阵的方式在R中加载部分数据，但是我们还是实现了一个”近似LSI”算法，使用一部分细胞进行最初的降维分析。近似LSI有两个主要用途，(1)提高降维速度，(2)初始降维时使用的细胞数的减少会降低数据的粒度，而这种粒度的减少有利于减少数据中的批次效应。但要注意，这可能会掩盖真实的生物学现象，因此近似LSI方法需要在人为监督下谨慎的使用。

4.2 隐语义(Latent Semantic Indexing)迭代

在scRNA-seq降维分析中(例如PCA)，一般都会鉴定变异基因。之所以这样做，是因为高变异基因通常更有生物学意义，并能够减少实验噪音。scATAC-seq数据只有0/1两种情况，所以你不能为降维分析提供变异peak。除了无法使用变异更大的peak，我们也尝试使用了更加开放的特征作为LSI的输入，然而经过多个样本的测试，我们发现处理结果存在较高的噪音，并且很难重复。为了解决该问题，我们提出了”迭代LSI法” (Satpathy*, Granja* et al. Nature Biotechnology 2019和 Granja*, Klemm* and McGinnis* et al. Nature Biotechnology 2019). 这个方法首先在最开放的分块中计算初始的LSI变换，识别未受到批次效应影响的低分辨率聚类。例如，我们可以在PBMC中鉴定出几个主要大群(T细胞、B细胞和单核细胞)。接着，ArchR计算这些聚类在每个特征中的平均开放度。ArchR识别这些聚类中变化最大的peak，用作后续LSI的特征。在第二次迭代中，变化最大的peak就类似于scRAN-seq中高变异基因。用户可以决定进行多次轮的LSI迭代。我们发现这个方法能够降低观察到的批次效应，并能使用更加合理的特征矩阵进行降维。

LSI过程

我们可以是用ArchR的addIterativeLSI()函数执行迭代LSI分析。虽然默认参数适用于绝大部分情况，但是我们还是鼓励你使用?addIterativeLSI看看有哪些可用的参数，以及它们对你的数据集的影响。一般经常修改的参数是iterations,varFeatures和resolution。注意，LSI的结果并不是确定的，也就是即便你用相同参数处理同一个数据，也会得到不一样的结果。当然只是略微不同，总体上还是相似的。因此，一旦你得到一个相对理想的降维结果，请及时保存你的ArchRProject项目或者相关的LSI信息。

处于教学目的，我们会创建一个名为”IterativeLSI”的reduceDims对象

1
2
3
4
5
6
7
8
9
10
11
12
13
projHeme2 <- addIterativeLSI(
    ArchRProj = projHeme2,
    useMatrix = "TileMatrix", 
    name = "IterativeLSI", 
    iterations = 2, 
    clusterParams = list( #See Seurat::FindClusters
        resolution = c(0.2), 
        sampleCells = 10000, 
        n.start = 10
    ), 
    varFeatures = 25000, 
    dimsToUse = 1:30
)

如果你在下游分析时看到一些批次效应，其中一个选择就是增加LSI的迭代次数，并以更低的分辨率进行聚类。此外，也可以通过降低特征变量的数目让程序关注变异更大的特征。

1
2
3
4
5
6
7
8
9
10
11
12
13
projHeme2 <- addIterativeLSI(
    ArchRProj = projHeme2,
    useMatrix = "TileMatrix", 
    name = "IterativeLSI2", 
    iterations = 4, 
    clusterParams = list( #See Seurat::FindClusters
        resolution = c(0.1, 0.2, 0.4), 
        sampleCells = 10000, 
        n.start = 10
    ), 
    varFeatures = 15000, 
    dimsToUse = 1:30
)

4.3 近似LSI

对于特别大的数据集，ArchR使用LSI投影得到近似的LSI降维结果。这个步骤和之前的迭代LSI流程相似，但是在LSI步骤存在一些不同。首先，我们会随机选择一部分”路标细胞”用于LSI降维分析。接着，使用”路标细胞”计算的逆文档频率对余下的细胞进行TF-IDF标准化。然后，这些标准化的细胞被投影到”路标细胞”定义的SVD子空间中。最终，我们就将余下的细胞根据少量的”路标细胞”实现了LSI变换和投影。因为ArchR不需要将所有细胞的read都保存在内存中，而是不断的提取每个样本的细胞，将其投影到路标细胞的LSI空间中，所以近似LSI算法非常高效。因此，即便之后用到了更大的数据集，近似LSI也能够减少内存的使用。

注意，路标细胞的数量取决于数据集中不同细胞的比例。

近似LSI

近似LSI通过修改ArchR的addIterativeLSI()函数的sampleCellsFinal和projectCellsPre参数来实现.samplesCellsFinal决定了路标细胞的数目，projectCellsPre告诉ArchR使用路标细胞对其余细胞进行投影。

4.4 使用Harmony矫正批次效应

在某些情况下，迭代的LSI方法并不能矫正严重的批次效应。ArchR使用了最初为了scRNA-seq开发的Harmony进行批次效应矫正。 ArchR提供了专门的封装函数addHarmony()将降维后的对象直接传递给Harmony::HarmonyMatrix(), 额外参数可以通过...进行传递。更详细的信息可以用?addHarmony()进行了解。用户需要根据他们的特定目的注意批次效应校正的结果。

1
2
3
4
5
6
7
8
9
10
projHeme2 <- addHarmony(
    ArchRProj = projHeme2,
    reducedDims = "IterativeLSI",
    name = "Harmony",
    groupBy = "Sample"
)
## Harmony 1/10
## Harmony 2/10
## Harmony 3/10
## Harmony converged after 3 iterations

这一步会在我们的projHeme2创建一个新的名为Harmony的reduceDims对象，可以通过projHeme2@reducedDims$Harmony访问。

R

使用ArchR分析单细胞ATAC-seq数据(第三章)

第3章: 创建ArchRProject

ArchRProject对象允许我们将多个Arrow文件整理到单个项目之中。ArchRProject比较小能够直接存放在内存中。通过操作ArchRProject，我们能够快速读取Arrow文件里的数据，并将修改后的数据写回到Arrow文件里。因此你几乎可以在这一章中找到后续分析所要用到的所有函数。此外，由于ArchRProject文件能够进行保存并在之后进行读取，那么也就保证了分析的连续性，同时也方便合作者之间的相互交流。这一章主要介绍如何创建ArchRProject对象并和它进行交互。

3.1 创建一个ArchRProject

ArchRProject至少需要设置两个参数，ArrowFiles和outputDirectory. 其中ArrowFilesArrow文件的文件路径，如果没有，清先阅读第一章将BAM/Fragments文件转成Arrow文件。outputDirectory指的是下游分析得到的文件的保存路径。ArchR会自动将之前提供的geneAnnotation和genomeAnnotation和新的ArchRProject项目进行关联。这些在我们在之前运行addArchRGenome("hg19")的时候被保存在环境变量中。

1
2
3
4
5
projHeme1 <- ArchRProject(
  ArrowFiles = ArrowFiles,
  outputDirectory = "HemeTutorial",
  copyArrows = TRUE #This is recommened so that if you modify the Arrow files you have an original copy for later usage.
)

因为这是我们hematopoiesis项目的第一轮分析，所以我们将ArchRProject命名为”proJHeme1”。在后续的分析中，我们会不断修改并更新这个ArchRProject对象，我们需要使用不同的变量名对项目进行跟踪。如果变量名始终都是同一个，那么就可能出现明明是同一个命令，有些时候能运行成功有些时候运行失败的情况，原因可能就是某些中间步骤被漏掉了，导致同一个变量名的背后的实际内容其实是不一样的。

我们可以直接检查下ArchRProject的内容

1
2
3
4
5
6
7
8
9
10
projHeme1
#输出信息
# class: ArchRProject
# outputDirectory: /path/to/HemeTutorial
# samples(3): scATAC_BMMC_R1 scATAC_CD34_BMMC_R1 scATAC_PBMC_R1
# sampleColData names(1): ArrowFiles
# cellColData names(15): Sample TSSEnrichment ... DoubletEnrichment BlacklistRatio
# numberOfCells(1): 10661
# medianTSS(1): 16.832
# medianFrags(1): 3050

从上面的输入中，我们可以发现ArchRProject在初始化的时候增加了一些额外信息

  1. outputDirectory: 输出目录
  2. samples： 样本名
  3. sampleColData: 记录每个样本的元信息
  4. cellColData: 记录每个细胞的元信息，第二章计算的”DoubleEnrichment”,”DoubletScore”就存放在此处
  5. numberOfCells: 项目总的合格细胞数，也就是不包括之前没有通过质控，或者被认为doublet的细胞
    1.medianTSS/medianFrags: 所有细胞的TSS值和Fragments数的中位数

我们可以通过下面这行代码了解下ArchRProject在R中会使用多少内存

1
2
paste0("Memory Size = ", round(object.size(projHeme1) / 10^6, 3), " MB")
# "Memory Size = 37.477 MB"

我们还可以检查当前的ArchRProject中存放了哪些矩阵数据，这些数据一旦增加后就可以在下游分析中使用。

1
2
getAvailableMatrices(projHeme1)
# "GeneScoreMatrix" "TileMatrix"

3.2 操作ArchRProject对象

既然已经在上一节创建了ArchRProject，这一节就来介绍一些比较常用的操作来从该对象中获取我们需要的信息和保存我们的数据

例1: 使用$访问cellColData

我们可以使用$访问细胞名(cellNames), 样本名(sample)和TSS富集得分(TSS enrichment Scores)

1
2
3
head(projHeme1$cellNames)
head(projHeme1$Sample)
quantile(projHeme1$TSSEnrichment)

除了这些以外，你还可以利用R的补全功能，在输入projHeme1$后看看它还能接哪些内容。

例2: 从ArchRProject中提取部分细胞

我们可以将以前学习的R语言向量/矩阵/数据框提取数据的方法无缝的迁移到ArchRProject上，但区别在于ArchR并不会直接提取数据，而是构建一个新的ArchRProject对象。

最简单的方法是根据位置和细胞名

1
2
3
4
# 根据位置
projHeme1[1:100, ]
# 根据细胞名
projHeme1[projHeme1$cellNames[1:100], ]

复杂一些就是先根据逻辑判断提取细胞名，接着根据细胞名提取数据。例如挑选scATAC_BMMC_R1的细胞，或是TSSEnrichment大于8的细胞。

1
2
3
4
5
6
7
8
# sample name
idxSample <- BiocGenerics::which(projHeme1$Sample %in% "scATAC_BMMC_R1")
cellsSample <- projHeme1$cellNames[idxSample]
projHeme1[cellsSample, ]
# TSS enrichment score
idxPass <- which(projHeme1$TSSEnrichment >= 8)
cellsPass <- projHeme1$cellNames[idxPass]
projHeme1[cellsPass, ]

例3: 在ArchRProject增加数据

我们可以通过在cellColData增加新列的方式来为我们项目中的细胞增加额外的元信息。

例如，我们可以从原来的样本名中提取更易阅读的部分添加到cellColData中。

1
2
3
bioNames <- gsub("_R2|_R1|scATAC_","",projHeme1$Sample)
head(bioNames)
#[1] "BMMC" "BMMC" "BMMC" "BMMC" "BMMC" "BMMC"

一种方式是使用$直接往cellColData添加列

1
projHeme1$bioNames <- bioNames

或者是用addCellColData()函数往cellColData中添加新的列。使用addCellColData()的好处在于，你可以只为部分增加信息。

1
2
3
4
bioNames <- bioNames[1:10]
cellNames <- projHeme1$cellNames[1:10]
projHeme1 <- addCellColData(ArchRProj = projHeme1, data = paste0(bioNames),
    cells = cellNames, name = "bioNames2")

ArchR默认会用NA填充其余部分。正因如此，如果我们对比bioNames和bioNames2时，你就可以发现NA填充了那些没有被选择的部分

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
getCellColData(projHeme1, select = c("bioNames", "bioNames2"))
DataFrame with 10661 rows and 2 columns
#                                     bioNames   bioNames2
#                                  <character> <character>
#scATAC_BMMC_R1#TTATGTCAGTGATTAG-1        BMMC        BMMC
#scATAC_BMMC_R1#AAGATAGTCACCGCGA-1        BMMC        BMMC
#scATAC_BMMC_R1#GCATTGAAGATTCCGT-1        BMMC        BMMC
#scATAC_BMMC_R1#TATGTTCAGGGTTCCC-1        BMMC        BMMC
#scATAC_BMMC_R1#TCCATCGGTCCCGTGA-1        BMMC        BMMC
#...                                       ...         ...
#scATAC_PBMC_R1#GCTGCGAAGATCCGAG-1        PBMC          NA
#scATAC_PBMC_R1#GCAGCTGGTGGCCTTG-1        PBMC          NA
#scATAC_PBMC_R1#GCAGATTGTACGCAAG-1        PBMC          NA
#scATAC_PBMC_R1#TTCGTTACATTGAACC-1        PBMC          NA
#scATAC_PBMC_R1#CGCTATCGTGAGGTCA-1        PBMC          N

例4: 从cellColData中提取列

ArchR提供了getCellColData用于从ArchRProject中提取元信息列。你可以认为这是$的功能增强版，因为$只能提取一列，而getCellColData可以提取多列信息，此外还有许多神奇的操作。

例如我们可以根据列名获取数据，比方说每个细胞的唯一fragment数

1
df <- getCellColData(projHeme1, select = "nFrags")

甚至，你还能在列名中添加一些运算操作，这样会直接输出运算后的结果

1
df <- getCellColData(projHeme1, select = c("log10(nFrags)", "nFrags - 1"))

例5: 绘制质控信息和TSS富集得分的对比图

根据上述提供的案例，我们可以很容易获取到每个细胞的标准scATAC-seq质控信息。我们发现目前最靠谱的质控标准是TSS富集的得分(评估ATAC-seq数据的信噪比)和唯一比对数(比对数如果不够，那么该细胞也没有分析的价值)

我们先用getCellColData提取我们需要的两列，其中nFrages需要进行log10运算

1
df <- getCellColData(projHeme1, select = c("log10(nFrags)", "TSSEnrichment"))

接下来，我们基于TSS富集得分和唯一比对进行绘图，函数ggPoint是ArchR对ggplot2的点图相关函数的封装。

1
2
3
4
5
6
7
8
9
10
11
12
p <- ggPoint(
    x = df[,1],
    y = df[,2],
    colorDensity = TRUE,
    continuousSet = "sambaNight",
    xlabel = "Log10 Unique Fragments",
    ylabel = "TSS Enrichment",
    xlim = c(log10(500), quantile(df[,1], probs = 0.99)),
    ylim = c(0, quantile(df[,2], probs = 0.99))
) + geom_hline(yintercept = 4, lty = "dashed") + geom_vline(xintercept = 3, lty = "dashed")

p

dot plot

该图的主要作用是找到高质量的细胞。你可能会注意到一些细胞已经在我们之前创建Arrow文件时被设置的阈值filterTSS, filterFargs过滤掉了。然而，如果我们觉得之前的质控过滤对于这个样本不太合适，我们可以根据该图调整阈值，有必要的话还可以重新产生Arrow文件。

尽管可以使用PDF的方式保存上图，但是对于这种点特别多的图，还是用PNG形式比较好。

1
2
3
4
png("TSS-vs-Frags.png")
plot(p)
dev.off()
getwd()

3.3 为ArchRProject每个样本绘制统计信息

如果一个项目中有多个不同的样本，我们很自然地就会想比较这些样本。所谓一图胜千言，ArchR提供了两种小提琴图(violin plot)和山脊图(ridge plot)用来展示不同组之间的信息。这两种类型的图可以用一个函数plotGroups进行绘制。当然这函数除了根据样本进行分组绘图外，还可以使用下游分析得到的分组信息（例如聚类）。

例1: 根据TSS富集得分为每个样本绘制山脊图。

绘制山脊图时，需要设置plotAs = "ridges"

1
2
3
4
5
6
7
8
p1 <- plotGroups(
    ArchRProj = projHeme1,
    groupBy = "Sample",
    colorBy = "cellColData",
    name = "TSSEnrichment",
    plotAs = "ridges"
   )
p1

ridges of TSSEnrichment

例2: 根据TSS富集得分为每个样本绘制小提琴图。ArchR绘制的小提琴图额外叠加了一层箱线图(a box-and-whiskers plot in the style of Tukey)，中间水平的三条线表示数据中的1/4分位数，中位数和3/4分位数。最下方是最小值，最上方是最大值（或者是1.5倍的IQR, interquartile range, 1/4分位数与3/4分位数的距离）

绘制小提琴图时，需要设置plotAs = "violin"

1
2
3
4
5
6
7
8
9
p2 <- plotGroups(
    ArchRProj = projHeme1,
    groupBy = "Sample",
    colorBy = "cellColData",
    name = "TSSEnrichment",
    plotAs = "violin",
    alpha = 0.4,
    addBoxPlot = TRUE
   )

Violin plot

例3: 根据log10(unique nuclear fragments)为每个样本绘制山脊图。

1
2
3
4
5
6
7
8
p3 <- plotGroups(
    ArchRProj = projHeme1,
    groupBy = "Sample",
    colorBy = "cellColData",
    name = "log10(nFrags)",
    plotAs = "ridges"
   )
p3

ridge plot

例4: 根据log10(unique nuclear fragments)为每个样本绘制小提琴图。

1
2
3
4
5
6
7
8
9
10
p4 <- plotGroups(
    ArchRProj = projHeme1,
    groupBy = "Sample",
    colorBy = "cellColData",
    name = "log10(nFrags)",
    plotAs = "violin",
    alpha = 0.4,
    addBoxPlot = TRUE
   )
p4

Violin Plot

我们可以使用plotPDF()将上面的图绘制在一个PDF中

1
plotPDF(p1,p2,p3,p4, name = "QC-Sample-Statistics.pdf", ArchRProj = projHeme1, addDOC = FALSE, width = 4, height = 4)

3.4 绘制样本的TSS富集谱和Fragment大小分布

鉴于数据的存储和读取方式， ArchR能够快速的从Arrow文件中计算出fragment大小分布和TSS富集谱。

Fragments大小分布: 我们使用plotFragmentSizes函数为绘制所有样本的fragment大小分布。ATAC-seq数据中的fragment大小分布可能会因样本、细胞类型和批次不同，但这和数据质量并不是严格相关。比如说我们的数据在不同组之间就存在一些差别。

1
2
p1 <- plotFragmentSizes(ArchRProj = projHeme1)
p1

Fragment size distribution

TSS富集谱: 我们使用plotTSSEnrichment()函数绘制TSS富集谱。TSS富集谱需要在中心位置有一个明显的峰，在峰的右边还会有一个核小体引起的小隆起（约147 bp）。

1
2
p2 <- plotTSSEnrichment(ArchRProj = projHeme1)
p2

TSS enrichment profiles

和之前一样，我们可以使用plotPDF()将图形以可编辑的矢量图进行保存。

1
plotPDF(p1,p2, name = "QC-Sample-FragSizes-TSSProfile.pdf", ArchRProj = projHeme1, addDOC = FALSE, width = 5, height = 5)

3.5 保存和加载ArchRProject

ArchR提供了saveArchRProject函数让我们能将内存中ArchRProject对象的保存到本地，并在之后进行读取，当然保存到本地的数据也就能分享给其他人。从根本上来讲，ArchRProject是一个指向多个Arrow文件的对象，因此使用saveArchRProject()保存ArchRProject的过程实际上分为两步

  1. 将当前的Arrow文件复制到目标输出目录下outputDirectory,确保它们能和新的ArchRProject关联
  2. 将目标ArchRProject保存到outputDirectory目录下

我们这里以projHeme1为例介绍saveArchRProject()的用法

1
saveArchRProject(ArchRProj = projHeme1, outputDirectory = "Save-ProjHeme1", load = FALSE)

运行结束后，outputDirectory对应的目录下就会保存相应的Arrow文件和projHeme1的Rds文件。

重要知识点: 这一步骤不会自动更新你当前R session里活跃的ArchRProject对象。也就说，当前R session里的projHeme1对象还是指向原来的Arrow文件，而不会指向拷贝到outputDirectory目录下的Arrow文件。

换言之，如果我们想要将当前的项目复制到新的目录，你可以设置load=TRUE。下面的代码运行结束后得到的projHemeNew就是指向拷贝后Arrow文件的ArchRProject对象。你可以使用原来的projHeme1，这就会覆盖当前R session的projHeme1。

1
projHemeNew <- saveArchRProject(ArchRProj = projHeme1, outputDirectory = "Save-ProjHeme1", load = TRUE)

3.6 使用ArchRProject过滤doublets

当我们使用addDOubletScores()增加了预测的doublet信息后，我们就能使用filterDoublets()过滤这些预测的doublets。这个函数里的一个关键参数是filterRatio, 由于不同样本的制备过程中的细胞上样量不同，那么不同样本的doublet也就不同，为了保证不同样本的过滤后结果保持一致，所以引入了filterRatio。该值越高，被认为是doublet而被过滤的细胞也就越多。

了解了filterRatio的作用后，我们再来看它的定义，它是根据未被过滤的细胞数计算的需要被过滤doublet的最大值(the maximum ratio of predicted doublets to filter based on the number of pass-filter cells)。举个例子。假如这里有5000个细胞，需要被过滤doublet的最大值就等于 filterRatio * 5000^2 / (100000), 也就是filterRatio * 5000 * 0.05。

如果还是不太了解，也不太纠结，毕竟大部分时候我们用的都是默认参数。我们现在只需要知道这个filterRatio参数很重要，得根据后续结果进行调整。这里就直接运行filterDoublets对细胞进行过滤。出于教学目的，我们将处理结果保存为一个新的ArchRProject对象，也就是projHeme2，实际分析中，你可以使用原来的名字覆盖之前的结果，降低内存消耗。

1
2
3
4
5
projHeme2 <- filterDoublets(projHeme1)
# Filtering 410 cells from ArchRProject!
#	scATAC_BMMC_R1 : 243 of 4932 (4.9%)
#	scATAC_CD34_BMMC_R1 : 107 of 3275 (3.3%)
#	scATAC_PBMC_R1 : 60 of 2454 (2.4%)

之前的projHeme1有10,661个细胞，经过上一步的过滤，projHeme2剩下了10,251个细胞，也就意味着有410个细胞(3.85%)是doublet。而且每个样本的过滤比例并不相同，这就是filterRatio的作用。

1
2
3
4
5
projHeme2
# ...
# numberOfCells(1): 10251
# medianTSS(1): 16.856
# medianFrags(1): 2991

如果你想过滤更多ArchRProject里的细胞，你可以设置更高的filterRatio。和参数调整的更多信息可以用?filterDoublets查看。

1
2
3
4
5
projHemeTmp <- filterDoublets(projHeme1, filterRatio = 1.5)
# Filtering 614 cells from ArchRProject!
#	scATAC_BMMC_R1 : 364 of 4932 (7.4%)
#	scATAC_CD34_BMMC_R1 : 160 of 3275 (4.9%)
#	scATAC_PBMC_R1 : 90 of 2454 (3.7%)

最后，既然projHemeTmp只是说明用的临时对象，我们现在可以将其从R session中删除了。

1
rm(projHemeTmp)
R

使用ArchR分析单细胞ATAC-seq数据(第二章)

第2章: 使用ArchR推断Doublet

单细胞数据分析中的一个重要问题就是”doublet”对分析的影响。”doublet”指的是单个液滴(droplet)捕获了一个条形码珠(barcode bead)和多个细胞核。这会导致原本来自于多个细胞的read结果被当成一个细胞，结果原来两个细胞被平均成一个细胞。我们在这一章中将会介绍如何使用计算的方法鉴定并过滤doublet。

2.1 ArchR是如何鉴定doublet

几乎所有平台得到的单细胞数据都或多或少的存在doublet。所谓的doublet就是一个液滴中混入了多个细胞的细胞核，导致原本的多个细胞被认为是一个细胞。在10X Genomics平台，doublets的细胞比例和一次反应中的细胞数有关，细胞越多，doublet也就越多。不过即便你按照标准流程，也依旧会存在一定比例的doublet，而只要有5%的doublet数据，就会对后续的聚类分析造成影响，尤其是后续的发育/谱系分析。因为doublet看起来像是两类细胞的中间状态，但是当你不知道这是doublet造成的假象时，你就会得到错误的结论。

为了预测哪些细胞才是真的doublet，我们通过对不同细胞进行组合得到模拟的doublet(in silico doublets)，然后将这些模拟的doublet和原来的细胞一同投影到UMAP，接着确定和它们最近的细胞。通过不断地迭代对该步骤，我们就可以找到数据中那些与模拟doublet信号最接近的细胞。这些细胞就是潜在的doublet。

Doublet分析原理

为了开发ArchR的doublet鉴定算法并验证它的可靠性，我们将10个不同类型的细胞进行混样测序。在没有doublet的情况下，我们的scATAC-seq数据最终应该有10个不同的细胞类型。但是当我们故意在10X Genomics scATAC-seq试剂中加入过量的细胞(25,000/per reaction), 我们会得到许多doublet。我们使用demuxlet根据一个细胞里是否包括两种不同细胞的基因型来判断该细胞是否是doublet。

colored by cell type

如上图所示，我们的真实结果(ground truth)被预测的doublet覆盖。ROC曲线中的AUC(area under the curve) > 0.90（见注1)

grouded truth
ArchR通过计算的方法移除这些doublets后，我们数据整体结构有了明显的变化，符合预期，也就是有10个不同的细胞类型。

Post-ArchR doublet Removal

注1: ROC曲线全称为受试者工作特征曲线 （receiver operating characteristic curve），它是根据一系列不同的二分类方式（分界值或决定阈），以真阳性率（敏感性）为纵坐标，假阳性率（1-特异性）为横坐标绘制的曲线。AUC就是衡量学习器优劣的一种性能指标，AUC越接近于1.0，表示检测方法真实性越高， 等于0.5表示真实性越低

2.2 使用ArchR推断scATAC-seq doublets

ArchR默认使用ArchR手稿中的doublet参数。我们建议用户检查移除doublet前后的数据变化，理解移除doulet对细胞的影响，根据结果对已有的参数进行调整，而不是生搬硬我们的参数设置。我们也会展示一些

我们使用ArchR的addDoubleScores()函数来移除doublet。对于教程中使用的数据，每个样本大约需要2到5分钟时间进行处理，最后每个细胞的doublet得分会添加到Arrow文件中。 你可以使用?addDoubletScores了解该函数的每个参数的意义。

1
2
3
4
5
6
doubScores <- addDoubletScores(
    input = ArrowFiles,
    k = 10, #Refers to how many cells near a "pseudo-doublet" to count.
    knnMethod = "UMAP", #Refers to the embedding to use for nearest neighbor search with doublet projection.
    LSIMethod = 1
)

我们以其中一个样本运行过程的输出日志信息进行介绍，这里的R^2用于描述样本中的细胞异质性，如果该数值非常小（例如小于0.9），说明该样本的细胞都非常相似。那么使用模拟的方法去鉴定doublet就不太合适了。这个很好理解，如果所有细胞都表达一个基因，并且表达量是1，那么你模拟的doublet也会只有一个细胞，且表达量是均值1，结果就是所有细胞都是doublet。在这种情况下，我们推荐跳过doublet预测这一步。或者你可以尝试设置knnMethod = "LSI",force = TRUE，在LSI子空间中进行投影。（相当于提高分辨率）。无论如何，你都需要手动检查结果，确保运行过程符合预期。

1
2
3
4
5
## If there is an issue, please report to github with logFile!
## 2020-04-15 09:28:44 : Batch Execution w/ safelapply!, 0 mins elapsed.
## 2020-04-15 09:28:44 : scATAC_BMMC_R1 (1 of 3) : Computing Doublet Statistics, 0.001 mins elapsed.
## scATAC_BMMC_R1 (1 of 3) : UMAP Projection R^2 = 0.9736
## scATAC_BMMC_R1 (1 of 3) : UMAP Projection R^2 = 0.9736

这里的R^2大于0.9，表示细胞存在异质性。

addDoubletScores在计算doublet得分的时候，还会在”QualityControl”中为每个样本生成三张图（一个样本一个PDF）

  1. Doublet Enrichments: 假设doublet符合均匀分布，那么每个细胞附近的模拟doublet数目都差不多。如果一个细胞附近相对于其他细胞有更多的doublet, 就认为它富集了doublet
  2. Doublet Score: 基于均匀分布假设，计算doublet富集显著性，以-log10(binomial adjusted p-value)进行展示。我们更推荐根据doublet enrichment鉴定doublet。
  3. Doublet Density: 模拟的duoblet在二维空间的投影，我们可以直观的了解我们模拟的doublet的分布情况。

BMMC

CD34 BMMC

PBMC