Lab 4: QC for snATAC-seq data

sc/snATAC-seq Analysis Pipeline

Notes

The estimated time for this lab is around 1h.

Aims

• Load 10X Genomics Multiome ATAC data into R.
• Create a Signac ChromatinAssay object.
• Add ATAC data to an existing Seurat object.
• Calculate ATAC quality control metrics and filter.
• Run the SCT, ATAC, and WNN analysis workflow.

5.1 Connect to RStudio Server

Open a browser and copy-paste the following address:

URL

http://10.35.229.71:8787/

Or click here: http://10.35.229.71:8787/

An RStudio log in page will appear; to log in, use your user ID for both ID and password.

5.2 Load Libraries

5.2.1 (RStudio env)

Load required R packages

Click

# Load
library(Seurat)
library(Signac)
library(ensembldb)
library(BSgenome.Ggallus.Ensembl.GRCg7bCustom)
library(ggplot2)
library(dplyr)
library(readr)
library(RColorBrewer)
library(hdf5r)
library(colorspace)

5.3 Define Variables and Load Data

5.3.1 (RStudio env)

Define file paths and load the data

Click

#######
# Variables
#######
h5_file <- "/data/processed/cellranger/E4_Dev/filtered_feature_bc_matrix.h5"
gtf_file <- "/data/shared/source/Gallus_gallus.bGalGal1.mat.broiler.GRCg7b.110.gtf"
frag_file <- "/data/processed/cellranger/E4_Dev/atac_fragments.tsv.gz"
seqinfo_file <- "/data/shared/source/GRCg7bSeqInfo.csv"
inputdata.10x.filtered <- Read10X_h5(h5_file)
Multiome <- readRDS("/data/shared/source/E4_seurat_obj.rds") #object with RNAseq data
# Extract peaks count data from the filtered dataset
atac_counts <- inputdata.10x.filtered$Peaks

5.4 Prepare Genome Annotation

5.4.1 (RStudio env)

Build gene annotations and chromosome (sequence) information for the GRCg7b genome, which the ChromatinAssay needs to map ATAC peaks

Click

# Get annotations to GRCg7b
edb <- EnsDb(ensDbFromGtf(gtf = gtf_file))
grange <- StringToGRanges(rownames(atac_counts), sep = c(":", "-"))
annotations <- GetGRangesFromEnsDb(ensdb = edb, standard.chromosomes = F)
genome(annotations) <- "GRCg7b"
frag.file <- frag_file
seqInfo <- read.csv(seqinfo_file)
seqInfo <- Seqinfo(seqInfo$seqnames, seqlengths=seqInfo$length,
                   isCircular=seqInfo$isCircular, genome="GRCg7b")

5.5 Create ATAC assay

5.5.1 (RStudio env)

Create a Signac ChromatinAssay

Click

# Create a ATAC assay
chrom_assay <- CreateChromatinAssay(counts = atac_counts,
  sep = c(":", "-"),
  genome = seqInfo,
  fragments = frag.file,
  min.cells = 5,
  annotation = annotations)

5.6 Add ATAC Assay to Seurat Object

5.6.1 (RStudio env)

Subset and add the ATAC assay

Click

# Subsets the ChromatinAssay object to include only the cells found in combined
chrom_assay <- subset(chrom_assay, cells = colnames(Multiome))
# Adds the ChromatinAssay
Multiome[["ATAC"]] <- chrom_assay
DefaultAssay(Multiome) <- "ATAC"

5.7 Calculate ATAC Quality Control Metrics

5.7.1 (RStudio env)

Calculate nucleosome signal and TSS enrichment

Click

# Calculates the nucleosome signal for each cell
Multiome <- NucleosomeSignal(Multiome)
# Calculates the transcription start site (TSS) enrichment for each cell
Multiome <- TSSEnrichment(Multiome, process_n = 2000)

5.8 Filtering

5.8.1 (RStudio env)

Plot QC metrics and subset

Click

# Filtering:
VlnPlot(Multiome, features = c("nCount_ATAC", "TSS.enrichment", "nucleosome_signal"), ncol = 3, pt.size = 1)
# Subset
Multiome <- subset(x = Multiome, subset = nCount_ATAC > 1000 & nCount_ATAC < 40000 &
    nucleosome_signal < 1 & TSS.enrichment > 1.5)

5.9 SCT Analysis

5.9.1 (RStudio env)

Run SCTransform and PCA

Click

# SCT analysis
Multiome <- SCTransform(Multiome, verbose = FALSE) |> RunPCA(reduction.name = "pca.sct")

5.10 ATAC Analysis

5.10.1 (RStudio env)

Run TF-IDF, top features, SVD and UMAP

Click

# ATAC analysis
DefaultAssay(Multiome) <- "ATAC"
Multiome <- RunTFIDF(Multiome)
Multiome <- FindTopFeatures(Multiome)
Multiome <- RunSVD(Multiome)

5.11 WNN Analysis

For the weighted nearest neighbor (WNN) analysis, follow the official Seurat tutorial: https://satijalab.org/seurat/articles/weighted_nearest_neighbor_analysis#wnn-analysis-of-10x-multiome-rna-atac