Error: "No matching chromosomes found in fragment file" while calculating gene activity

[rylanci]

Hello. I am attempting to integrate two merged objects, one sci-atac the other snRNA, and having an issue while calculating gene activity. Despite the error it appears the chromosomes in the frag files are consistent with the annotations. Not sure if I need to merge fragment files for the atac component as in the 0.2.5 signac merging vignette. It seems that method was phased out. Any advice would be appreciated. Thanks!

Here is the code and traceback.

gene.activities <- GeneActivity(atac.sobj, features = VariableFeatures(atac.sobj))

1. GeneActivity(atac.sobj, features = VariableFeatures(atac.sobj))
2. FeatureMatrix(fragments = frags, features = transcripts, cells = cells, 
 .     verbose = verbose, ...)
3. sapply(X = obj.use, FUN = function(x) {
 .     SingleFeatureMatrix(fragment = fragments[[x]], features = features, 
 .         cells = cells, sep = sep, verbose = verbose, process_n = process_n)
 . })
4. lapply(X = X, FUN = FUN, ...)
5. FUN(X[[i]], ...)
6. SingleFeatureMatrix(fragment = fragments[[x]], features = features, 
 .     cells = cells, sep = sep, verbose = verbose, process_n = process_n)
7. stop("No matching chromosomes found in fragment file.")

Here are samples of the annotation and fragments.

Annotation:
GRanges object with 6 ranges and 5 metadata columns:
seqnames ranges strand | tx_id
|
ENSMUSE00001236884 chr3 3508030-3508332 + | ENSMUST00000108393
ENSMUSE00000676606 chr3 3634150-3634347 + | ENSMUST00000108394
ENSMUSE00001345708 chr3 3638059-3638230 + | ENSMUST00000108393
ENSMUSE00001345708 chr3 3638059-3638230 + | ENSMUST00000108394
ENSMUSE00000149313 chr3 3641223-3641317 + | ENSMUST00000108393
ENSMUSE00000149313 chr3 3641223-3641317 + | ENSMUST00000108394
gene_name gene_id gene_biotype type

ENSMUSE00001236884 Hnf4g ENSMUSG00000017688 protein_coding exon
ENSMUSE00000676606 Hnf4g ENSMUSG00000017688 protein_coding exon
ENSMUSE00001345708 Hnf4g ENSMUSG00000017688 protein_coding exon
ENSMUSE00001345708 Hnf4g ENSMUSG00000017688 protein_coding exon
ENSMUSE00000149313 Hnf4g ENSMUSG00000017688 protein_coding exon
ENSMUSE00000149313 Hnf4g ENSMUSG00000017688 protein_coding exon
seqinfo: 22 sequences from mm10 genome

Fragment object (one of six present in merged atac object):

chrom start end barcode readCount
chr1 3000826 3000898 LL_7_TACTTTGCGCAAAGGAACAGAC 1
chr1 3001689 3001808 LL_7_CCTATGTGTAACCGTCCCATTC 1
chr1 3001738 3001792 LL_7_GGATTGCCTTCTTAGGTAGCGT 1
chr1 3001738 3001956 LL_7_GCCAACTGAGGGTTGTCCTCTG 1
chr1 3003199 3003367 LL_7_ACTCCTATGGAAATAAGGCCAG 1
chr1 3003574 3003641 LL_7_CAACCTCATTTGACACTGTGCG 1

Our OS is CentOS. Here is the session info:

R version 4.1.0 (2021-05-18)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /home/rlancione/miniconda3/envs/renv4/lib/libopenblasp-r0.3.15.so

locale:
[1] C

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] EnsDb.Mmusculus.v79_2.99.0 ensembldb_2.16.4          
 [3] AnnotationFilter_1.16.0    GenomicFeatures_1.44.0    
 [5] AnnotationDbi_1.54.1       Biobase_2.52.0            
 [7] Matrix_1.3-4               gridExtra_2.3             
 [9] ggplot2_3.3.5              dplyr_1.0.7               
[11] GenomicRanges_1.44.0       GenomeInfoDb_1.28.1       
[13] IRanges_2.26.0             S4Vectors_0.30.0          
[15] BiocGenerics_0.38.0        Signac_1.3.0              
[17] SeuratObject_4.0.2         Seurat_4.0.3              

[timoast]

My guess is that VariableFeatures(atac.sobj) is an empty vector, or a set of peak names rather than gene names. I've improved the error checking in the GeneActivity() function to give more informative error messages (on develop branch, see instructions here: https://satijalab.org/signac/articles/install.html#development-version-1)

You do not need to merge fragment files on-disk

[rylanci]

I see. You are correct VariableFeatures() is returning peaks. I did not realize that was an issue due to peaks being the currently defined feature set in the object. It made sense to me the variable features would be peaks.

In this case how would VariableFeatures be pointed towards the genes present in our annotations vs our peaks?

[timoast]

See the documentation for GeneActivity(): https://satijalab.org/signac/reference/geneactivity

The features parameter should be a list of gene names present in the gene annotation that you want to quantify, or you can leave it blank to quantify all protein-coding genes in the gene annotation