BIOINF 545 single-nucleus multiome lab

Introduction

In this lab we’ll learn to process single-nucleus multiome (ATAC + gene expression) data generated on the popular 10x Genomics platform.

Many researchers use 10x Genomics’ cellranger software to process this data. Cellranger attempts to be an easy, all-in-one solution for basic processing of 10x Genomics single-cell and single-nucleus data, taking raw sequencing reads as input and outputting BAM files, cell type clusters, and some downstream analysis results. But certain steps of the analysis, such as quality control, are difficult to reliably perform in an automated fashion and this can lead to suboptimal results. Therefore, we’ll take a hands-on approach.

The dataset we’ll be using is an example dataset published by 10x, composed of human peripheral blood mononuclear cells (PBMCs), a subset of immune cells.

Singularity

Throughout this lab you’ll make use of Singularity containers. If you have heard of Docker containers, Singularity containers are a similar idea. They allow one to package up software and necessary dependencies into a single unit. As long as Singularity is installed on a computer, one can simply move a Singularity container onto the machine, and use the software inside the container without having to install the software or it’s dependencies on the host machine itself. This allows for seamless transfer of dependencies between different machines and eases analysis reproducibility by ensuring that software versions remain constant.

To briefly demonstrate this, try to run the MACS2 peak calling software:

macs2 --version

Unless you’ve installed macs2 previously, this command will fail with the error macs2: command not found.

We’ll use macs2 later in this lab, and a Singularity container that includes an installation of macs2 is available on greatlakes. So let’s try to run macs2 using this container instead:

# ensure Singularity is loaded. Singularity is already installed on greatlakes.
module load singularity/3.10.4

# for convenience, set a variable defining the path to the lab materials
export LAB_ROOT=/nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab

# start up a new shell inside the Singularity container:
singularity shell ${LAB_ROOT}/singularity/porchard-default-general-20220107.img

# try to run macs2
macs2 --version

Now macs2 is available to run. Once you’re finished in the container, you can exit the container using exit or CTRL+D.

Preprocessing of 10x multiome data

Much of the initial processing of 10x multiome data can be done independently for the RNA and ATAC components. In this section, we’ll examine preprocessing of the ATAC component, and we’ll describe preprocessing of the RNA component after that.

Preprocessing of the ATAC component

Many of the processing steps involved for snATAC are similar to the processing done on bulk ATAC-seq data. But some steps are performed on a per-10x barcode basis rather than a per-library basis, and a few steps are unique to bulk or single-nucleus data:

In addition, as a QC step it’s often useful to process the snATAC-seq data as if it were bulk ATAC-seq data, in order to verify that the single-nucleus library as a whole displays similar properties as bulk ATAC-seq data, such as TSS enrichment.

We’ve packaged all of these snATAC-seq processing steps into a single pipeline (https://github.com/porchard/snATACseq-NextFlow). This pipeline is built using a popular bioinformatics workflow manager, NextFlow. NextFlow is essentially a framework for the specification of a processing pipeline, and it can simplify pipeline reuse across datasets and pipeline portability across computing environments. Snakemake is another popular workflow manager you may have encountered before. We’ll employ our NextFlow pipeline here to process the snATAC-seq component of a multiome dataset (and to demonstrate the usefulness of a good workflow manager!).

To explore the individual steps of the pipeline and the general anatomy of a NextFlow pipeline, please have a look at the main pipeline file, which can be browsed at: https://github.com/porchard/snATACseq-NextFlow/blob/master/main.nf

You’ll notice larger sections in this file labeled ‘process’. Each process defines a single step in the pipeline. For example, here’s the process that defines a job to map sequencing reads with bwa, with some comments added:

process bwa {

    memory '50 GB'        // memory required for the job
    cpus 12               // number of CPUs to use for the job
    errorStrategy 'retry' // what to do if the job fails
    maxRetries 1          // in case of failure, don't try to re-run the job more than once
    time '48h'            // time requested to run the job
    container 'library://porchard/default/bwa:0.7.15' // the (cloud-hosted) Singularity container to use

    // this defines the input expected for the job
    // in this case, a library identifier, the name of the reference genome, the name of the sequencing readgoup, and two fastq files
    input:
    tuple val(library), val(genome), val(readgroup), path(fastq_1), path(fastq_2)

    // this defines the output expected for the job
    // in this case, the library identifer, the readgroup name, the reference genome name, and a BAM file
    output:
    tuple val(library), val(readgroup), val(genome), path("${library}-${readgroup}-${genome}.bam")

    // this defines the shell command to run
    """
    bwa mem -I 200,200,5000 -M -t 12 ${get_bwa_index(genome)} ${fastq_1} ${fastq_2} | samtools sort -m 1g -@ 11 -O bam -T sort_tmp -o ${library}-${readgroup}-${genome}.bam -
    """

}

At the very bottom of the file, you’ll notice a section called workflow. Briefly, this section just defines how the processes are strung together; for example, this is where you would tell NextFlow to take the fastq files, pass them to the ‘trim’ process (which trims sequencing adapters), then pass that output as input to the ‘bwa’ process (which maps the trimmed reads).

Now let’s see how one configures and runs a NextFlow pipeline. Given our pipeline definition file and some information about our data, NextFlow will run the desired processing steps for us, including submitting all of the processing jobs to greatlakes’ job queue if requested. In order to do this, NextFlow requires some details about the computing environment. A suitable configuration file is available at /nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab/pipelines/snATACseq-NextFlow/nextflow.config:

// require a minimum NextFlow version
nextflowVersion = '>=20.10.0'

// output details about the resource usage of each job
trace.enabled = true
report.enabled = true

// if one wishes to use SLURM, uncomment this section
// process {
//     executor='slurm'
//     clusterOptions='--account=bioinf545w23_class'
//     containerOptions='--bind "/nfs:/nfs" --bind "/gpfs:/gpfs"'
//     shell = ['/bin/bash', '-ueo', 'pipefail']
// }

// if one wishes to execute jobs within your active terminal session, uncomment this section
// for the sake of this demo we will do this, to avoid jobs sitting in the SLURM queue and because
// we'll be running the pipeline from within an interactive SLURM job anyway
process {
    executor='local' // each job should be run on the local machine
    containerOptions='--bind "/nfs:/nfs" --bind "/gpfs:/gpfs"'
    shell = ['/bin/bash', '-ueo', 'pipefail']
}

// if running jobs locally rather than through the SLURM queue, 
// don't use more than one CPU and 11GB memory
executor {
    name = 'local'
    cpus = 1
    memory = '11 GB'
}

// Our pipeline uses singularity containers at each step
// We'll set a singularity cache location to avoid unnecessarily re-downloading containers that are already available on greatlakes
singularity.enabled = true
singularity.autoMounts = true
singularity.cacheDir = '/nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab/pipelines/snATACseq-NextFlow/singularity-cache'

//
// Below are the paths to necessary reference files (chromosome size files, BWA index for mapping, etc)
//

params.blacklist = ["hg19": ["/nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab/data/reference/hg19.blacklist.1.bed.gz", "/nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab/data/reference/hg19.blacklist.2.bed.gz"]]
params.chrom_sizes = ["hg19": "/nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab/data/reference/hg19.chrom_sizes"]
params.bwa_index = ["hg19": "/nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab/data/reference/bwa-index/hg19/hg19"]
params.tss = ["hg19": "/nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab/data/reference/hg19.tss.refseq.bed.gz"]
params.plot_signal_at_genes = ['GAPDH', 'MYH1', 'MYH2', 'MYH7'] // if these genes exist in the TSS file, ATAC signal near their TSS will be visualized
// this is the barcode whitelist for the snATAC component of 10x Genomics multiome experiments
params['barcode-whitelist'] = "/nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab/pipelines/snATACseq-NextFlow/737K-arc-v1.txt"

We also need to tell the pipeline what reference genome to use to process the data, and where the raw sequencing data (fastq files) is located. A config file with this information is at /nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab/pipelines/snATACseq-NextFlow/library-config.json:

{
    "libraries": {
        "pbmc_granulocyte_sorted_10k": {
            "genome": [
                "hg19"
            ],
            "readgroups": {
                "L001": {
                    "1": "/nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab/data/fastq/atac/pbmc_granulocyte_sorted_10k_S16_L001_R1_001.fastq.gz",
                    "2": "/nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab/data/fastq/atac/pbmc_granulocyte_sorted_10k_S16_L001_R3_001.fastq.gz",
                    "index": "/nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab/data/fastq/atac/pbmc_granulocyte_sorted_10k_S16_L001_R2_001.fastq.gz"
                }
            }
        }
    }
}

This file says that there is only one library (named pbmc_granulocyte_sorted_10k), with a single readgroup (named L001) that we’d like to process. It says to use hg19 as the reference genome, and it points towards three fastq files. The first two fastq files are the two paired-end sequencing reads. The third fastq file tells us the 10x barcode that corresponding to each read pair in the other fastq files.

Note that, even though the pipeline uses many software packages that are not installed by default on greatlakes, we don’t need to install those software packages ourselves. This is because the NextFlow pipeline leverages Singularity containers that our lab previously built and placed in the cloud.

Now we can run the pipeline. For the sake of time we’ll just be processing a subset of ~45k snATAC-seq reads, which will probably take about 2 minutes to go through the pipeline:

# load a recent version of java (necessary for NextFlow)
module load openjdk/18.0.1.1

# request interactive SLURM session. We'll use this session for the rest of the lab
salloc --account=bioinf545w24_class --mem-per-cpu=20G --time=02:00:00

# launch the pipeline
# --results defines where the output of the pipeline should be placed
${LAB_ROOT}/nextflow run \
    -resume \
    -c ${LAB_ROOT}/pipelines/snATACseq-NextFlow/nextflow.config \
    -params-file ${LAB_ROOT}/pipelines/snATACseq-NextFlow/library-config.json \
    --results /nfs/turbo/dcmb-class/bioinf545/sec001/${USER}/multiome-lab/nextflow-pipeline-results \
    ${LAB_ROOT}/pipelines/snATACseq-NextFlow/main.nf

As the pipeline runs, it will display a list of all the processes with the current progress. Once it finishes, have a look at the --results directory (which should be /nfs/turbo/dcmb-class/bioinf545/sec001/{your_unique_name}/multiome-lab/nextflow-pipeline-results). You should see the following directories:

ataqv  bigwig  bwa-corrected-barcodes  corrected-barcodes  macs2  mark_duplicates  merge  plot-barcodes-matching-whitelist  prune  transformed-barcodes  trim

Preprocessing of the RNA component

The key steps of RNA preprocessing are:

1. Map sequencing reads
2. Correct sequencing errors in 10x barcodes (and UMIs)
3. Collapse PCR duplicates using UMIs
4. Compute gene count matrix (i.e., a matrix showing the number of deduplicated read counts each gene has for each 10x barcode)
5. Compute QC metrics for each 10x barcode

For the sake of time we won’t actually run an RNA processing pipeline, but we do have a NextFlow pipeline for this available at: https://github.com/porchard/snRNAseq-NextFlow. This pipeline leverages starsolo (part of the STAR RNA-seq aligner software) to perform steps 1-4 above. It computes per-barcode QC metrics (# of UMIs, fraction of UMIs from the mitochondrial chromosome, etc) using a small custom python script.

Selection of pass-QC barcodes

After the above processing of the ATAC and RNA components of a multiome library, one must perform QC of the data.

We’ll practice this using previously-generated NextFlow pipeline results for our example library. Pipelines results for the ATAC and RNA components of our example multiome dataset can be found at /nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab/data/full-run/atacseq and /nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab/data/full-run/rnaseq, respectively.

In a typical 10x multiome library, less than 10k nuclei are profiled. But there are > 700k unique 10x barcodes in the 10x multiome kit, and usually hundreds of thousands of these barcodes will appear at least once in the dataset. As these numbers suggest, the great majority of barcodes that appear in a library do not represent quality nuclei. These junk barcodes are capturing loose nucleic acids that have escaped from ruptured cells and nuclei, and must be filtered out prior to clustering and downstream analysis. Therefore, during QC, we’ll attempt to determine which barcodes represent quality nuclei.

So, let’s load up the QC metrics produced in the RNA and ATAC pipelines, and use these to figure out which barcodes represent quality nuclei:

# we'll use a singularity container with an R installation and some installed R libraries for this
singularity shell --bind /nfs/turbo/dcmb-class:/nfs/turbo/dcmb-class \
                    ${LAB_ROOT}/singularity/porchard-default-r-general-20220112.img

# start up R
R

Now, within R:

library(dplyr)
library(tidyr)
library(glue)
library(ggplot2)

# load the per-barcode quality control metrics
LAB_ROOT <- '/nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab'

# ATAC QC metrics
atac_qc <- read.table(
    file.path(LAB_ROOT, 'data/full-run/atacseq/ataqv/single-nucleus/pbmc_granulocyte_sorted_10k-hg19.txt'),
     sep='\t', as.is=T, head=F, col.names=c('atac_barcode', 'metric', 'value'))

# RNA QC metrics
rna_qc <- read.table(
    file.path(LAB_ROOT, 'data/full-run/rnaseq/qc/pbmc_granulocyte_sorted_10k-hg19.qc.txt'),
     sep='\t', as.is=T, head=T) %>% dplyr::filter(barcode != '-')

# ATAC and RNA barcodes have a known 1-to-1 correspondence, but the barcodes aren't identical.
# The 10x barcode whitelists can be used to link the ATAC barcode to the corresponding RNA barcode.
# The first barcode in the RNA list corresponds to the first barcode in the ATAC list, etc.
RNA_WHITELIST <- file.path(LAB_ROOT, 'data/barcode-whitelists/rna-barcode-whitelist.txt')
ATAC_WHITELIST <- file.path(LAB_ROOT, 'data/barcode-whitelists/atac-barcode-whitelist.txt')

rna_barcodes <- read.table(RNA_WHITELIST, head=F)[,1]
atac_barcodes <- read.table(ATAC_WHITELIST, head=F)[,1]

atac_to_rna_barcode <- setNames(rna_barcodes, atac_barcodes)
atac_qc$rna_barcode <- atac_to_rna_barcode[atac_qc$atac_barcode]

# subset to the most informative ATAC QC metrics and reshape the table
atac_qc <- atac_qc[atac_qc$metric %in% c('tss_enrichment', 'hqaa'),] %>% tidyr::spread(key=metric, value=value)

# rename columns for clarity
atac_qc <- dplyr::rename(atac_qc, 
        atac_pass_filter_reads=hqaa, 
        atac_tss_enrichment=tss_enrichment)
rna_qc <- dplyr::rename_with(rna_qc, 
        function(x){glue('rna_{x}')})

# join ATAC and RNA columns
joint_qc <- rna_qc %>% dplyr::full_join(atac_qc)

## Joining, by = "rna_barcode"

# we now have this table of joint ATAC + RNA QC metrics
head(joint_qc)

Quality nuclei should have:

1. High sequencing coverage in both the RNA and ATAC components
2. High TSS enrichment (clearly above 1) in the ATAC component
3. Preferably low mitochondrial contamination

Let’s plot these metrics, starting with the RNA and ATAC sequencing coverage:

# plot RNA vs ATAC sequencing coverage
p <- ggplot(joint_qc) + 
    geom_point(aes(x=rna_umis, y=atac_pass_filter_reads), alpha=0.05) + 
    scale_x_continuous(trans = "log10") + 
    scale_y_continuous(trans = "log10") + 
    geom_hline(yintercept=1e4, color='red', linetype='dashed') +
    geom_vline(xintercept=1e3, color='red', linetype='dashed') +
    labs(x = "RNA UMIs", y = "ATAC pass filter reads") +
    theme_bw()

plot(p)

## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 79075 rows containing missing values (geom_point).

As expected, there is a population of barcodes with high RNA coverage (here, greater than 1000 UMIs as indicated by the vertical red line) and high ATAC coverage (here, greater than 10k pass-filter reads; horizontal red line).

Let’s ensure that this high sequencing coverage population has good TSS enrichment as well:

# plot ATAC pass filter reads vs ATAC TSS enrichment
p <- ggplot(joint_qc) + 
    geom_point(aes(x=atac_pass_filter_reads, y=atac_tss_enrichment), alpha=0.05) + 
    scale_x_continuous(trans = "log10") + 
    scale_y_continuous(trans = "log10") + 
    geom_hline(yintercept=5, color='red', linetype='dashed') +
    geom_vline(xintercept=1e4, color='red', linetype='dashed') +
    labs(x = "ATAC pass filter reads", y = "ATAC TSS enrichment") +
    theme_bw()

plot(p)

## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 523455 rows containing missing values (geom_point).

The TSS enrichment of the high coverage nuclei in the upper right is greater than 5, which is pretty good.

We expect this high coverage population to have relatively low mitochondrial fraction as well:

# plot RNA UMIs vs fraction of RNA reads from the mitochondrial chromosome
p <- ggplot(joint_qc) + 
    geom_point(aes(x=rna_umis, y=rna_fraction_mitochondrial), alpha=0.05) + 
    scale_x_continuous(trans = "log10") + 
    geom_hline(yintercept=0.2, color='red', linetype='dashed') +
    geom_vline(xintercept=1e3, color='red', linetype='dashed') +
    labs(x = "RNA UMIs", y = "RNA fraction mitochondrial") +
    theme_bw()

plot(p)

## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Removed 36129 rows containing missing values (geom_point).

The population of barcodes with high RNA coverage here (lower right quadrant) actually has higher mitochondrial contamination than we usually hope to see, possibly indicating the library prep was a little messy. Nevertheless, the mitochondrial contamination is not completely out of control (and as we’ll see in the ‘clustering’ section this data appears to yield reasonable cell type clusters), so this data will be fine for our purposes.

Let’s create a list of pass QC barcodes we can use in downstream analysis, thresholding on the variables explored in the above plots:

pass_qc_barcodes <- joint_qc %>% 
    dplyr::filter(rna_umis>=1e3 & 
                  atac_pass_filter_reads>=1e4 & 
                  rna_fraction_mitochondrial<=0.2 & 
                  atac_tss_enrichment>=5) %>% 
    pull(rna_barcode)

# how many barcodes pass QC?
length(pass_qc_barcodes)

## [1] 10298

Clustering

Once we know which barcodes represent quality nuclei, we can cluster those nuclei in order to characterize the cell types present in our data and assign each nucleus to one of these cell types.

There are a variety of clustering algorithms, some of which are built for a specific individual modality (only RNA or only ATAC) and some of which include functionality for clustering on more than one modality.

Clustering RNA data is usually easier than clustering ATAC data. For this multiome data, let’s try to cluster only on RNA, ignoring the ATAC modality (once RNA clusters are known, ATAC clusters are known too, since each RNA barcode maps to a single ATAC barcode). We’ll use the Seurat R package for this.

Identification of cell type clusters with Seurat

The RNA gene count matrix is available in the previously-generated NextFlow RNA pipeline results. This matrix is stored in market matrix format, which is useful for storing sparse matrices (in snRNA data, most genes in any given pass-QC nucleus will have 0 counts, so the data is very sparse indeed). It’s not important for this exercise, but for future reference market matrix format is defined here; the matrix actually consists of three individual files, which we’ll load into R for clustering:

library(Seurat)

## Attaching SeuratObject

# RNA gene counts are stored in a market matrix format file
# this function reads them into an R sparse matrix:
load_mm <- function(matrix_file, features_file, barcodes_file) {
    tmp <- as(Matrix::readMM(matrix_file), 'dgCMatrix')
    features <- read.table(features_file, as.is=T, sep='\t', head=F)
    features <- paste0(features$V1, ' (', features$V2, ')')
    barcodes <- read.table(barcodes_file, as.is=T, head=F)[,1]
    dimnames(tmp) <- list(features, barcodes)
    return(tmp)
}

COUNT_MATRIX_DIR <- '/nfs/turbo/dcmb-class/bioinf545/shared/multiome-lab/data/full-run/rnaseq/starsolo/pbmc_granulocyte_sorted_10k-hg19/pbmc_granulocyte_sorted_10k-hg19.Solo.out/GeneFull_ExonOverIntron/raw'

RNA_MTX <- file.path(COUNT_MATRIX_DIR, 'matrix.mtx')
RNA_FEATURES <- file.path(COUNT_MATRIX_DIR, 'features.tsv')
RNA_BARCODES <- file.path(COUNT_MATRIX_DIR, 'barcodes.tsv')

mtx <- load_mm(RNA_MTX, RNA_FEATURES, RNA_BARCODES)

# subset to our pass QC barcodes
mtx <- mtx[,pass_qc_barcodes]

# have a look at a small piece of the matrix
# rows are genes, columns are 10x barcodes
mtx[1:5,1:2]

## 5 x 2 sparse Matrix of class "dgCMatrix"
##                               CTCATTGTCCTCCTAA CAGGACACAACCCTAA
## ENSG00000162441.7 (LZIC)                     .                .
## ENSG00000202415.1 (RN7SKP269)                .                .
## ENSG00000156869.8 (FRRS1)                    .                .
## ENSG00000201491.1 (RNU4-75P)                 .                .
## ENSG00000202254.1 (Y_RNA)                    .                .

In Seurat, the basic clustering algorithm can be described in a few steps:

1. Normalize the data (for example, this helps adjust for the fact that some nuclei have many more reads than others)
2. Perform dimensionality reduction on the data. This is often done using PCA, and we’ll use PCA here.
3. Figure out which nuclei are nearby each other in PC space, and group them into clusters (this might be done using, for example, the Louvain clustering algorithm)
4. Visualize the nuclei and clusters by using UMAP or tSNE on the PCs.

We’ll do each of these in turn. First, we need to create a Seurat object from our count matrix:

# create the Seurat object
# ignore genes that appear in very few (less than 5) nuclei
rna <- CreateSeuratObject(counts = mtx, min.cells=5)

## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')

Now, we need to normalize our data to take into account different sequencing depths across nuclei:

rna <- suppressWarnings(SCTransform(rna, verbose = FALSE))

Now we can perform PCA on this normalized matrix. After doing this, we need to decide how many PCs to use in downstream analysis. Here, we’ll examine the amount of variance explained by the top 50 PCs by making an elbow plot, and then select the number of PCs to use based on where this curve levels off:

rna <- RunPCA(rna, npcs=50, verbose=F)
ElbowPlot(rna, ndims=50)

The elbow seems to be leveling off by around 25 PCs. So we’ll use 25 PCs for the rest of this analysis:

PCS <- 25

Now, we can find sets of nuclei which are similar in PC space, and cluster them together:

rna <- FindNeighbors(rna, dims = 1:PCS, k.param = 20)

## Computing nearest neighbor graph

## Computing SNN

rna <- FindClusters(rna, resolution=0.5)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 10298
## Number of edges: 368051
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9119
## Number of communities: 16
## Elapsed time: 0 seconds

We’ll additionally produce a UMAP embedding that we can use to visualize the clusters, and plot this UMAP embedding, coloring each nucleus by the cell type cluster to which it was assigned:

rna <- RunUMAP(rna, reduction='pca', dims=1:PCS)

## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

## 19:43:33 UMAP embedding parameters a = 0.9922 b = 1.112

## 19:43:33 Read 10298 rows and found 25 numeric columns

## 19:43:33 Using Annoy for neighbor search, n_neighbors = 30

## 19:43:33 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 19:43:34 Writing NN index file to temp file /tmp/RtmpIg4ZF8/file17b77d37141ebb
## 19:43:34 Searching Annoy index using 1 thread, search_k = 3000
## 19:43:37 Annoy recall = 100%
## 19:43:38 Commencing smooth kNN distance calibration using 1 thread
## 19:43:39 Initializing from normalized Laplacian + noise
## 19:43:40 Commencing optimization for 200 epochs, with 425646 positive edges
## 19:43:44 Optimization finished

DimPlot(rna, reduction = "umap")

Determination of cell type identities

Once we have clusters, we must determine which cell type is represented by each cluster.

This is commonly done using marker genes. There are software packages available that try to automate this process (e.g., Garnett), often using either marker genes or comparing your dataset to previously-annotated datasets using correlation or machine-learning type approaches. In practice, if one has a fairly limited number of cell types from a tissue that is well-explored, it’s easy to do cell type determination by hand, by (1) visualizing per-cluster expression of known marker genes, and (2) examining the list of genes that are most specific to each cluster and poking through the literature to see if these genes are known markers of any cell types or play a special role in any cell types that might be expected to appear in your dataset.

Let’s start by visualizing some known blood cell marker genes – MS4A1 (a marker of B cells), CD14 (CD14+ Monocytes), and CD8A (CD8+ T cells).

MARKERS <- c("MS4A1", "CD14", "CD8A")

# this next line is only necessary because our feature names in the input
# data are in the format gene_id (gene_name), so we need to map the gene names in the
# MARKERS variable to the corresponding feature names
PLOT_FEATURES <- unlist(lapply(glue('\\({MARKERS}\\)'), function(x){grep(x, rownames(mtx), value=T, ignore.case=T)}))

One useful function for this is the FeaturePlot function, which overlays (normalized) gene expression for our desired genes onto the UMAP we generated:

FeaturePlot(rna, features = PLOT_FEATURES)

Another useful function is the VlnPlot function, which generates a violin plot displaying the normalized expression of our desired genes:

VlnPlot(rna, features = PLOT_FEATURES)

In the above two figures, we see for example that MS4A1, a B cell marker, is expressed in clusters 5 and 8, suggesting that these clusters are composed of B cells.

In the case that you can’t identify the cell type of each cluster based on your prior marker genes, you can try to identify genes that are specific to each cluster and use these cluster-specific genes to try to identify the corresponding cell types using literature searches. Many clustering packages will include functionality that helps you identify the genes that have high impact in the clustering or are specific to certain cell types.

In Seurat, you can use the FindMarkers() and FindAllMarkers() functions for this. Let’s use the FindMarkers() function to try to figure out what the differences bewteen the two B cell clusters are (clusters 5 and 8):

markers_distinguishing_cluster_5_from_cluster_8 <- FindMarkers(rna, ident.1=5, ident.2=8)
head(markers_distinguishing_cluster_5_from_cluster_8, 20)

We see here that IL4R and TCL1A are expressed in the majority of nuclei in cluster 8 (pct.2) but not in many nuclei in cluster 5 (pct.1). Searching the literature, one finds that these two genes tend to be more highly expressed in naive B cells than in other B cell subtypes, suggesting that cluster 8 could be naive B cells. The literature also indicates that another major B cell subtype (memory B cells) is distinguished from naive B cells by high expression of CD27. If we plot the expression of these markers in clusters 5 and 8:

MARKERS <- c("IL4R", 'TCL1A', 'CD27')
PLOT_FEATURES <- unlist(lapply(glue('\\({MARKERS}\\)'), function(x){grep(x, rownames(mtx), value=T, ignore.case=T)}))
VlnPlot(rna, features = PLOT_FEATURES, idents=c(5, 8))

These three genes clearly distinguish one B cell cluster from ther other. Cluster 5, which expresses CD27, might be memory B cells while cluster 8, which expresses IL4R and TCL1A, represent naive B cells.

Examination of per-cluster QC metrics

Note that sometimes clustering may be impacted by technical variables, such as sample batch. It is therefore good practice to look for clustering by sample, batch, QC metrics, or any other variables that may indicate confounding of clustering by technical variables.

In this case, we can add our QC metrics to the Seurat object, and again use the VlnPlot function to visualize our QC metrics on a per-cluster basis to verify that none of the clusters are obvious outliers in terms of, for example, sequencing coverage. We don’t expect QC metrics to be perfectly concordant across clusters, in part because biological differences between cell types may result in QC differences (for example, a cell type with high RNA content per cell might show greater UMI counts in the RNA data)

# add QC metrics to Seurat object
rownames(joint_qc) <- joint_qc$rna_barcode
pass_qc_barcode_metrics <- joint_qc[pass_qc_barcodes,]
pass_qc_barcode_metrics$rna_umis_log10 <- log10(pass_qc_barcode_metrics$rna_umis)
rna <- AddMetaData(rna, pass_qc_barcode_metrics)
head(rna@meta.data)

# plot
VlnPlot(rna, features = c('rna_umis_log10', 'rna_fraction_mitochondrial'))

There are some differences across clusters; for example, nuclei in clusters 9, 10, 14, and 15 generally have higher RNA UMI counts than nuclei in other clusters. It’s difficult to know from just these plots why this is, and whether it’s a problem. Perhaps those clusters represent cell types known to have higher RNA content, or perhaps those clusters represent doublets or some other technical artifact. In an actual analysis, one might wish to dig into this further, paying particular attention to how those clusters behave in downstream analysis. If any of those clusters behave strangely in downstream analysis, show difficult-to-interpret marker gene patterns, or stand out in other suspicious ways, you might decide to flag that cluster as artifactual and remove those nuclei from further analysis.

Lastly, let’s save our barcode –> cluster assignments for downstream analysis:

# output clusters
clusters <- as.data.frame(Idents(rna))
colnames(clusters) <- c('cluster')
clusters$barcode <- rownames(clusters)
write.table(clusters[,c('barcode', 'cluster')], file = 'clusters.txt', append = F, quote = F, sep = '\t', row.names = F, col.names = F)

quit(save='no')

Per-cluster processing of ATAC data

Once we have the clusters, the ATAC data can be processed on a per-cluster basis for downstream analysis.

To do this, let’s leave the current Singularity container and hop into another one:

exit
singularity shell --bind /nfs/turbo/dcmb-class:/nfs/turbo/dcmb-class \
            ${LAB_ROOT}/singularity/porchard-default-general-20220107.img

For the sake of time, we’ll select a single cluster to process. Let’s take the naive B cell cluster, cluster 8:

# get RNA barcodes from cluster 8
grep -w 8 clusters.txt | cut -f1 > rna-barcodes-in-cluster-8.txt

# we need to convert these RNA barcodes to ATAC barcodes, since we'll use them to subset the full ATAC BAM file:
# note this command assumes that the RNA barcodes are all distinct from the ATAC barcodes, which is the case for these particular 10x whitelists
paste ${LAB_ROOT}/data/barcode-whitelists/rna-barcode-whitelist.txt \
      ${LAB_ROOT}/data/barcode-whitelists/atac-barcode-whitelist.txt \
      > rna-barcode-to-atac-barcode.txt

grep -f rna-barcodes-in-cluster-8.txt rna-barcode-to-atac-barcode.txt | cut -f2 > atac-barcodes-in-cluster-8.txt

# subset the full ATAC bam file to reads from these nuclei.
# for the sake of time, we'll only use reads from chromosome 21.
# barcodes are stored in the 'CB' tag in the BAM file.
samtools view -h -b -D CB:atac-barcodes-in-cluster-8.txt \
            ${LAB_ROOT}/data/full-run/atacseq/prune/pbmc_granulocyte_sorted_10k-hg19.pruned.bam \
            chr21 > cluster_8.bam

Now we can do standard ATAC analyses on this pseudobulk BAM file. Let’s call peaks using MACS2:

bedtools bamtobed -i cluster_8.bam > cluster_8.bed

macs2 callpeak -t cluster_8.bed \
                --outdir . \
                -f BED \
                -n cluster_8 \
                --SPMR \
                -g hs \
                --nomodel \
                --shift -100 \
                --seed 762873 \
                --extsize 200 \
                -B \
                --broad \
                --keep-dup all

Per-cluster visualization of ATAC data

MACS2 also outputs a bedGraph file, which displays continuous ATAC signal across the genome. We frequently convert these to bigWig files, which are convenient for visualizing ATAC-seq signal at specific genomic loci via the UCSC genome browser. Let’s generate a bigWig file we could use for this:

bedGraphToBigWig cluster_8_treat_pileup.bdg \
                ${LAB_ROOT}/data/reference/hg19.chrom_sizes \
                cluster_8.bw

For visualization in the genome browser, the bigwig must be placed on a world-accessible server so that it can be accessed via URL. In this case we’ll use two previously-generated bigwig files already sitting on such a server. If you’d like to visualize your own data and do not have access to such a hosting server, you may be able to get free web hosting services through the university.

The bigwig files that we’ll be visualizing are available at: https://theparkerlab.med.umich.edu/gb/custom-tracks/2020-sn-muscle-cell-types/0-hg19.bw and https://theparkerlab.med.umich.edu/gb/custom-tracks/2020-sn-muscle-cell-types/1-hg19.bw

To visualize them:

1. Navigate to the UCSC genome browser home page (https://genome.ucsc.edu/)
2. In the navigation bar at the top of the page, hover your mouse over ‘My Data’, then click ‘Custom Tracks’
3. This data is human data mapped to reference genome hg19, so select genome ‘Human’ and assembly ‘Feb. 2009 (GRCh37/hg19)’ from the pulldown menus near the top of the page
4. Now we need to give the genome browser ‘track lines’, which define the URLs at which each bigwig file can be found and gives the browser a little information about the file and how we’d like to view it. We’ll use these track lines:

track type=bigWig db=hg19 name='Type I muscle fibers' visibility=full color=68,58,131 alwaysZero=on maxHeightPixels=50:50:50 windowingFunction=mean smoothingWindow=3 autoScale=off viewLimits=0:3 bigDataUrl=https://theparkerlab.med.umich.edu/gb/custom-tracks/2020-sn-muscle-cell-types/1-hg19.bw

track type=bigWig db=hg19 name='Type II muscle fibers' visibility=full color=68,1,84 alwaysZero=on maxHeightPixels=50:50:50 windowingFunction=mean smoothingWindow=3 autoScale=off viewLimits=0:3 bigDataUrl=https://theparkerlab.med.umich.edu/gb/custom-tracks/2020-sn-muscle-cell-types/0-hg19.bw

Paste this line into the “Paste URLs or data” box.

5. Press the ‘submit’ button, and on the next page press ‘go to first annotation’.

You should see two tracks at the top of the browser, labeled ‘Type I muscle fibers’ and ‘Type II muscle fibers’, representing our two bigWig tracks. You can navigate around the genome as desired to visualize the ATAC signal in these two cell types. In this case, it might be interesting to visualize loci that distinguish these cell types, such as the MYH7 locus (which is a marker for type I fibers). If you type ‘MYH7’ into the top search bar, you can see that there is much more open chromatin at this locus in Type I muscle fibers than in Type II muscle fibers: