Unable to link peaks to genes for multiomics data

Question

[bgbs]

Hello,

I've been following the Joint RNA and ATAC analysis: 10x multiomic tutorial on the Stuart Lab website to process output from the 10X Genomics Cellranger ARC pipeline. I've been trying to run the following steps in the 'Link peaks to genes' section for one of my samples and for 1 or two genes, but I keep getting errors and I'm not sure how to fix them in order to make coverage plots. Any insight on how to fix this would be appreciated, thanks.

# chnage from RNAseq assay to ATACseq assay
DefaultAssay(liftoff_1_MI5_V1_SO) <- "ATAC"

# first compute the GC content for each peak (worked! Got a warning though)
liftoff_1_MI5_V1_SO <- RegionStats(liftoff_1_MI5_V1_SO, genome = BSgenome.Mfascicularis.NCBI.5.0)
# Warning: Not all seqlevels present in supplied genome

liftoff_1_MI5_V1_SO <- LinkPeaks(
  object = liftoff_1_MI5_V1_SO,
  peak.assay = "ATAC",
  expression.assay = "SCT",
  genes.use = c("OSTN", "FOS")
)

Processed 35548 groups out of 35548. 100% done. Time elapsed: 3s. ETA: 0s.
Testing 1 genes and 210517 peaks
  |                                                  | 0 % ~calculating  Warning: Requested more features than present in supplied data.
            Returning 0 featuresError in density.default(x = query.feature[[featmatch]], kernel = "gaussian",  :
argument 'x' must be numeric

# > liftoff_1_MI5_V1_SO
# An object of class Seurat 
# 271812 features across 8787 samples within 3 assays 
# Active assay: ATAC (210519 features, 210519 variable features)
#  2 layers present: counts, data
#  2 other assays present: RNA, SCT
#  5 dimensional reductions calculated: pca, lsi, umap, umap.atac, wnn.umap


# check if genes present in SCT expression assay
# > all(c("OSTN", "FOS") %in% rownames(liftoff_1_MI5_V1_SO[["SCT"]]))
# [1] TRUE


# Check the peak identifiers in the ATAC assay
# > head(rownames(liftoff_1_MI5_V1_SO[["ATAC"]]))
# [1] "chr1-30919-31783" "chr1-32491-33425" "chr1-36538-37474" "chr1-54567-55405" "chr1-57748-58645" "chr1-61338-62189"

I then tried to make a coverage plot with just one gene, but I also got an error.

liftoff_1_MI5_V1_SO <- LinkPeaks(
  object = liftoff_1_MI5_V1_SO,
  peak.assay = "ATAC",
  expression.assay = "SCT",
  genes.use = c("OSTN")
)

Error in density.default(x = query.feature[[featmatch]], kernel = "gaussian",  : 
  argument 'x' must be numeric

When I set genes.use to NULL to determine genes from expression assay, I get a similar error

liftoff_1_MI5_V1_SO <- LinkPeaks(
  object = liftoff_1_MI5_V1_SO,
  peak.assay = "ATAC",
  expression.assay = "SCT",
  genes.use = NULL
)

# Testing 23041 genes and 210517 peaks
#   |                                                  | 0 % ~calculating  Warning: Requested more features than present in supplied data.
#             Returning 0 featuresError in density.default(x = query.feature[[featmatch]], kernel = "gaussian",  : 
#   argument 'x' must be numeric

[yura.grabovska]

From the object ATAC assay:

> liftoff_1_MI5_V1_SO
# An object of class Seurat 
# 271812 features across 8787 samples within 3 assays 
# Active assay: ATAC (210519 features, 210519 variable features)

From the code error:

> Processed 35548 groups out of 35548. 100% done. Time elapsed: 3s. ETA: 0s.
Testing 1 genes and 210517 peaks

Your ATAC assay lists 210519 variable features and your code lists 210517 peaks. You're missing 2 features from your variable features.

Check and fix your variable features and see if that helps.

[bgbs]

Thanks for your quick reply. I did check if the genes are in the SCT data slot and they are.

# Extract the data slot of the SCT assay
sct_data <- GetAssayData(object = liftoff_1_MI5_V1_SO, assay = "SCT", slot = "data")

# Check if the genes are in the rownames
genes_to_check <- c("OSTN", "FOS")
genes_present <- genes_to_check %in% rownames(sct_data)

# Print results
names(genes_present) <- genes_to_check
genes_present


# OSTN  FOS 
# TRUE TRUE

How can I specifically check and fix the variable features in the ATAC assay of this sample? The number of variable features in the object is 210519. I'm not sure why I'm missing 2 features/genes when I'm using the LinkPeaks function? Is it the 2 genes I specified that are missing?