3 Data Loading

3.1 Setup

Load all required libraries for multiome (RNA + ATAC) data processing,quality control, and visualisation.

library(Seurat)
library(hdf5r)
library(fs)
library(scCustomize)
library(clustree)
library(shiny)
library(SeuratDisk)
library(Signac)
library(Matrix)
library(dplyr)
library(tidyr)

path_wd <- getwd()
setwd(path_wd)

3.2 Genome annotation

Query AnnotationHub for Ensembl 98 human gene annotations (GRCh38/hg38). These annotations are used to link ATAC peaks to nearby genes and to compute TSS enrichment scores. Sequence levels are set to UCSC style (chr1, chr2…) to match the cellranger output.

library(AnnotationHub)
ah <- AnnotationHub()
ensdbs <- query(ah, c("EnsDb.Hsapiens"))
ensdb_id <- ensdbs$ah_id[grep(paste0(" 98 EnsDb"), ensdbs$title)]
ensdb <- ensdbs[[ensdb_id]]
seqlevelsStyle(ensdb) <- "UCSC"
annotations <- GetGRangesFromEnsDb(ensdb = ensdb)
genome(annotations) <- "hg38"

3.3 Per-sample data loading

Load the filtered feature-barcode matrix from the 10x cellranger multiome output for each sample. A Seurat object is created from the RNA Gene Expression matrix. Donor identity is added from the SNP-based demultiplexing results (scSNPdemux). The ATAC Peaks matrix is added as a ChromatinAssay with the fragment file and hg38 genome annotations attached. ATAC features are restricted to standard chromosomes to remove unplaced contigs.

3.3.1 Sample B

raw.lib <- Read10X_h5("~/data/tasks/liam.kealy/processing/cellranger_SC_ATAC_covid_B/outs/filtered_feature_bc_matrix.h5")

cnts <- CreateSeuratObject(counts = raw.lib$`Gene Expression`,
                           assay = "RNA",
                           project = "SampleB",
                           names.delim = "-", names.field = 2)

# annotate the demultiplexing
demux <- read.table("~/data/tasks/liam.kealy/processing/scSNPdemux/cellranger_SC_ATAC_covid_B.demux/results/donor_ids.tsv",
                    header = T)
rownames(demux) <- demux$cell
cnts <- AddMetaData(object = cnts, metadata = demux)

# add ATAC assay
cnts[['ATAC']] <- CreateChromatinAssay(counts = raw.lib$`Peaks`,
                                       annotation = annotations,
                                       fragments = "~/data/tasks/liam.kealy/processing/cellranger_SC_ATAC_covid_B/outs/atac_fragments.tsv.gz",
                                       sep = c(":", "-"),
                                       genome = 'hg38')

# standard chromosomes
library(BSgenome.Hsapiens.UCSC.hg38)
standard_chroms <- standardChromosomes(BSgenome.Hsapiens.UCSC.hg38)
idx_standard_chroms <- which(as.character(seqnames(granges(cnts[['ATAC']]))) %in% standard_chroms)
cnts[["ATAC"]] <- subset(cnts[["ATAC"]],
                         features = rownames(cnts[["ATAC"]])[idx_standard_chroms])
seqlevels(cnts[['ATAC']]@ranges) <- intersect(seqlevels(granges(cnts[['ATAC']])),
                                               unique(seqnames(granges(cnts[['ATAC']]))))
cntsB <- cnts

3.3.2 Sample C

raw.lib <- Read10X_h5("~/data/tasks/liam.kealy/processing/cellranger_SC_ATAC_covid_C/outs/filtered_feature_bc_matrix.h5")

cnts <- CreateSeuratObject(counts = raw.lib$`Gene Expression`,
                           assay = "RNA",
                           project = "SampleC",
                           names.delim = "-", names.field = 2)

# annotate the demultiplexing
demux <- read.table("~/data/tasks/liam.kealy/processing/scSNPdemux/cellranger_SC_ATAC_covid_C.demux/results/donor_ids.tsv",
                    header = T)
rownames(demux) <- demux$cell
cnts <- AddMetaData(object = cnts, metadata = demux)

# add ATAC assay
cnts[['ATAC']] <- CreateChromatinAssay(counts = raw.lib$`Peaks`,
                                       annotation = annotations,
                                       fragments = "~/data/tasks/liam.kealy/processing/cellranger_SC_ATAC_covid_C/outs/atac_fragments.tsv.gz",
                                       sep = c(":", "-"),
                                       genome = 'hg38')

# standard chromosomes
idx_standard_chroms <- which(as.character(seqnames(granges(cnts[['ATAC']]))) %in% standard_chroms)
cnts[["ATAC"]] <- subset(cnts[["ATAC"]],
                         features = rownames(cnts[["ATAC"]])[idx_standard_chroms])
seqlevels(cnts[['ATAC']]@ranges) <- intersect(seqlevels(granges(cnts[['ATAC']])),
                                               unique(seqnames(granges(cnts[['ATAC']]))))
cntsC <- cnts

3.3.3 Sample D

raw.lib <- Read10X_h5("~/data/tasks/liam.kealy/processing/cellranger_SC_ATAC_covid_D/outs/filtered_feature_bc_matrix.h5")

cnts <- CreateSeuratObject(counts = raw.lib$`Gene Expression`,
                           assay = "RNA",
                           project = "SampleD",
                           names.delim = "-", names.field = 2)

# annotate the demultiplexing
demux <- read.table("~/data/tasks/liam.kealy/processing/scSNPdemux/cellranger_SC_ATAC_covid_D.demux/results/donor_ids.tsv",
                    header = T)
rownames(demux) <- demux$cell
cnts <- AddMetaData(object = cnts, metadata = demux)

# add ATAC assay
cnts[['ATAC']] <- CreateChromatinAssay(counts = raw.lib$`Peaks`,
                                       annotation = annotations,
                                       fragments = "~/data/tasks/liam.kealy/processing/cellranger_SC_ATAC_covid_D/outs/atac_fragments.tsv.gz",
                                       sep = c(":", "-"),
                                       genome = 'hg38')

# standard chromosomes
idx_standard_chroms <- which(as.character(seqnames(granges(cnts[['ATAC']]))) %in% standard_chroms)
cnts[["ATAC"]] <- subset(cnts[["ATAC"]],
                         features = rownames(cnts[["ATAC"]])[idx_standard_chroms])
seqlevels(cnts[['ATAC']]@ranges) <- intersect(seqlevels(granges(cnts[['ATAC']])),
                                               unique(seqnames(granges(cnts[['ATAC']]))))
cntsD <- cnts