TSSEnrichment() fails

[Seandelao]

Hello,

Great package. I am running Signac on my data, and am running into an issue at the TSSEnrichment step. When running the function, I get the following error:

Error in SetAssayData.ChromatinAssay(object = object[[assay]], slot = slot, :
Position enrichment must be provided as a matrix or sparseMatrix

I am able to follow the Vignette fine, so it's something with my data. Here is my code:

counts <- Read10X_h5(filename = "/Users/seandelao/Desktop/Signac/12wpc_epcam/filtered_peak_bc_matrix.h5")
metadata <- read.csv(file = "/Users/seandelao/Desktop/Signac/12wpc_epcam/singlecell.csv",
                     header = TRUE,
                     row.names = 1
)

chrom_assay <- CreateChromatinAssay(
  counts = counts,
  sep = c(":", "-"),
  genome = 'hg38',
  fragments = "/Users/seandelao/Desktop/Signac/12wpc_epcam/fragments.tsv.gz",
  min.cells = 10,
  min.features = 200
)

pbmc <- CreateSeuratObject(
  counts = chrom_assay,
  assay = "peaks",
  meta.data = metadata
)
pbmc
pbmc[['peaks']]
granges(pbmc)

annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75)

seqlevelsStyle(annotations) <- 'UCSC'
genome(annotations) <- "hg38"

Annotation(pbmc) <- annotations

pbmc <- NucleosomeSignal(object = pbmc)
pbmc

pbmc <- TSSEnrichment(object = pbmc, fast = FALSE)

[timoast]

Here you're using the wrong gene annotations for hg38, so it's probably looking in the wrong places for the TSS sites. EnsDb.Hsapiens.v86 is the correct version for hg38

[Seandelao]

Thank you! I changed to EnsDb.Hsapiens.v86 and am now getting this error for the TSSEnrichment step:

Extracting TSS positions
Finding + strand cut sites
Finding - strand cut sites
Error in validObject(r) :
invalid class "ngTMatrix" object: Dim slot must have length 2

[timoast]

This should now be fixed on the develop branch. See the website for installation instructions: https://satijalab.org/signac/articles/install.html#development-version-1

[Seandelao]

Thanks @timoast, I am now able to get through the TSSEnrichment function after installing the develop branch. However, I am getting a score of essentially zero (see attached plots). I've ran this data through ArchR and observe a normal looking TSS Enrichment score, so I'm not sure if it's something I'm doing wrong? Here is my current code:

''#Pre-processing workflow
#Peak/Cell matrix and Fragment file
counts <- Read10X_h5(filename = "/Users/seandelao/Desktop/Signac/12wpc_epcam/filtered_peak_bc_matrix.h5")
metadata <- read.csv(file = "/Users/seandelao/Desktop/Signac/12wpc_epcam/singlecell.csv",
header = TRUE,
row.names = 1
)

chrom_assay <- CreateChromatinAssay(
counts = counts,
sep = c(":", "-"),
genome = 'hg38',
fragments = "/Users/seandelao/Desktop/Signac/12wpc_epcam/fragments.tsv.gz",
min.cells = 10,
min.features = 200
)

pbmc <- CreateSeuratObject(
counts = chrom_assay,
assay = "peaks",
meta.data = metadata
)

pbmc
pbmc[['peaks']]
granges(pbmc)

#Get gene annotations for the peaks and add to the object
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotations) <- 'UCSC'
genome(annotations) <- "hg38"
Annotation(pbmc) <- annotations

#QC metrics
pbmc <- NucleosomeSignal(object = pbmc)
pbmc <- TSSEnrichment(object = pbmc, fast = FALSE, assay = "peaks")
pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments / pbmc$passed_filters * 100
pbmc$blacklist_ratio <- pbmc$blacklist_region_fragments / pbmc$peak_region_fragments
pbmc@meta.data$TSS.enrichment

#Plot the TSS Enrichment score
pbmc$high.tss <- ifelse(pbmc$TSS.enrichment > 2, 'High', 'Low')
TSSPlot(pbmc, group.by = 'high.tss') + NoLegend()
pbmc$nucleosome_group <- ifelse(pbmc$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
FragmentHistogram(object = pbmc, group.by = 'nucleosome_group')

VlnPlot(
object = pbmc,
features = c('pct_reads_in_peaks', 'peak_region_fragments',
'TSS.enrichment', 'blacklist_ratio', 'nucleosome_signal'),
pt.size = 0,
ncol = 5
)''
Rplot09.pdf
Rplot10.pdf

For reference, here is the score I'm getting with ArchR:

Rplot11.pdf

[timoast]

My guess is that you're using the wrong set of gene annotations. What genome did you map the reads to? How big is the fragment file and how many cells do you have? Can you show the first few lines of the fragment file?

[Seandelao]

The files are output from 10x's CellRanger, which mapped the reads from GRCH38 (based on the web summary). The fragment file is 2.58GB, have about 6.5k cells. If I do rownames(chrom_assay) after reading in the fragment file, I get:

[1] "chr1-181266-181696" "chr1-191245-191738" "chr1-629711-630163" "chr1-633794-634264"
[5] "chr1-778172-779753" "chr1-804779-804932" "chr1-826632-827932" "chr1-832099-832522"
[9] "chr1-833415-833519" "chr1-858140-858980" "chr1-869494-870388" "chr1-876749-877499"

etc.

[Seandelao]

Would I use the:
seqlevelsStyle(annotations) <- 'UCSC'

line for GRCh38? If not, is there a different 'seqlevelsStyle()' that I should set it to?

[Seandelao]

@timoast turns out I somehow got my wires crossed and mixed up sample files, it seems to be working now. Thank for your help!