RuntimeError: neither 'fragment_single' nor 'fragment_paired' is present in the '.obsm'

[jayluo2]

Dear CEMBA team,

Thank you for this great resource. I am attempting to save subclass-specific fragments files (as BED files) for one of the h5ad files available under this folder — specifically, the file “CEMBA180111_4E_rm_dlt.h5ad” — and am getting the following error (I am using SnapATAC 2.5.1):

RuntimeError                              Traceback (most recent call last)
Cell In[29], line 1
----> 1 snapatac2.ex.export_fragments(anndata_test,
      2                               groupby='Subclass',
      3                              out_dir=anndata_dir)

File ~/miniconda3/envs/snapatac2/lib/python3.9/site-packages/snapatac2/export/__init__.py:68, in export_fragments(adata, groupby, selections, ids, out_dir, prefix, suffix, compression, compression_level)
     65     elif suffix.endswith(".zst"):
     66         compression = "zstandard"
---> 68 return internal.export_fragments(
     69     adata, list(ids), list(groupby), out_dir, prefix, suffix, selections, compression, compression_level,
     70 )

RuntimeError: neither 'fragment_single' nor 'fragment_paired' is present in the '.obsm'

when running

snapatac2.ex.export_fragments(adata,
                              groupby='Subclass',
                             out_dir=anndata_dir)

where ‘adata’ is loaded with snapatac2.read(), and ’Subclass’ is a column of subclass labels retrieved from Supplementary Table 2 (downloaded from this folder). The error occurs with both backed=None and backed=‘r+’.

I verified that the keys are indeed missing by checking anndata_test.obsm, which returns AxisArrays with keys: insertion. I am wondering if there is an error in my data-loading/handling?

Best, Jay

The subclass column under .obs is as follows.

[beyondpie]

@jayluo2 Thank you for letting me know this. I will reply to you within one week. @wkl1990 Do you save subclass-level bed files ?

[jayluo2]

Thank you! I am also wondering if it would be possible to provide stranded and unshifted (no Tn5 +4/-5 correction) bedpe files (generated from BAM files) at the cell type or subclass level, like the ones provided here?

[wkl1990]

@jayluo2 Thank you for letting me know this. I will reply to you within one week. @wkl1990 Do you save subclass-level bed files ?

Yes, I exported the subclass-level bed file using SnapATAC2. I will put them in the Catlas later.

[jayluo2]

Thank you, we really appreciate it. Do you know when this will be uploaded to catlas.org so we can look out for it? Can’t wait to access this exciting dataset!

It would also be great to have class-level fragments bed files!

[beyondpie]

@jayluo2 Sorry for the late response.

I tested the data you mentioned, in both my local PC (SnapATAC2 version 2.3.1)and our cluster (SnapATAC2 version 2.4.0). Both of them are OK.

Note, when I use SnapATAC2, I notice that I don't have the function snapatac2.ex.export_fragments, instead I have the function snapatac2.ex.export_bed.

So not sure if this is the problem of SnapATAC2?

Would you mind to try SnapATAC2 with the version I have?

@wkl1990 When you have time, please also give this example a try. I think this might be SnapATAC2 version issue. And if possible, would you mind to share the codes about exporting bed files with us? I think, this will be a common request for the users. And we may need to write it in our README.

Thanks! Songpeng

[beyondpie]

@jayluo2

BTW, in order to get the bed files for all the cells in a subclass, you can try to use AnnDataSet to merge all the 234 AnnData into one AnnDataSet, then use the groupby to export the bed files.

Also, our paper is published now: https://www.nature.com/articles/s41586-023-06824-9

Best, Songpeng

[jayluo2]

Thank you! snapatac2. ex.export_bed() from SnapATAC 2.4.0 works for writing BED files (e.g. below), but there does not seem to be a column for strand information. Is there a way to get either unshifted/uncorrected Tn5 insertion sites or an additional column with ‘+’/‘-‘ (strand) for these fragments?

head -5 052_Pvalb_Gaba.bed
chr1    3649238 3649239 AGCGATAGAACCAGGTTTCATCCACCTATCCT
chr1    3649523 3649524 AGCGATAGAACCAGGTTTCATCCACCTATCCT
chr1    7188790 7188791 AGCGATAGAACCAGGTTTCATCCACCTATCCT
chr1    7188956 7188957 AGCGATAGAACCAGGTTTCATCCACCTATCCT
chr1    7583296 7583297 AGCGATAGAACCAGGTTTCATCCACCTATCCT

We are also wondering if a class-to-subclass mapping can be provided? I see the Fig. 1 subclass ordering on Supplementary Table 3 but it would be great to have class-level annotations in the cell metadata as well!

Best, Jay

[beyondpie]

@jayluo2

1. The fragment here has no "+/-" strand information after SnapATAC2 preprocessing. Why do you need those information? The fragments here are already shifted to remove the bias of Tn5 insertion.
2. If you want to get the raw sequences before correction, you can perform alignment using our raw fastq (all the raw fastq files can be downloaded from NCBI GEO: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE246791) files, and get the bam files, I think?
3. For the corresponding class information, you can check Allen's metadata here: AIT21_annotation_freeze_081523.tsv under https://github.com/beyondpie/CEMBA_wmb_snATAC/blob/main/02.integration/src/main/resource) Feel free to read the files under that directory for other information.

Songpeng

[beyondpie]

@jayluo2 @wkl1990

I've asked Kai about this. It is the SnapATAC2 version. There are some break changes in SnapATAC2 >= 2.5. And all the data I generated are under SnapATAC2 <= 2.4. So if you want to use the h5ad files I shared, you have to use SnapATAC2 <= 2.4.

I will add this note in README.

https://kzhang.org/SnapATAC2/changelog.html

Thanks! Songpeng

[jayluo2]

Thank you for your help!