tanaylab / mcATAC

Metacell analysis for ATAC data
https://tanaylab.github.io/mcATAC/
Other
1 stars 0 forks source link

Are metacell numbers 0- or 1-based #34

Closed yonatans2 closed 2 years ago

yonatans2 commented 2 years ago

I think they should be 1-based. But it should at least be consistent...

> wd <- "/net/mraid14/export/tgdata/users/yonshap/proj/mmcortex/"
> setwd(wd)
> library(metacell)
> scdb_init('scdb')
initializing scdb to scdb
> devtools::load_all("/home/feshap/src/mcATAC")
ℹ Loading mcATAC
ℹ Parallelization enabled. Using 77 threads.
> mca <- mcATAC::import_from_h5ad(file = '/net/mraid14/export/tgdata/users/yonshap/proj/mmcortex/data/test_anndata_atac_metacell.h5ad', genome = 'mm10')
• Reading /net/mraid14/export/tgdata/users/yonshap/proj/mmcortex/data/test_anndata_atac_metacell.h5ad
! h5ad file doesn't have the mc_size_eps_q at the uns section. Using the default: 0.1
→ No id was given, setting id to Jennifer_Carter
• Setting egc cell size to 11312.1833139651 (the 0.1 quantile of metacell sizes)
✔ Successfully loaded an `McATAC` object with 454 metacells and 98905 ATAC peaks
> mc_rna <- scdb_mc('pl_cort')
> mca <- add_mc_rna(atac_mc = mca, mc_rna = mc_rna)
Warning messages:
1: mc_rna contains 1 metacells not present in the McATAC object: "454"
2: McATAC object contains 1 metacells not present in mc_rna: "0"
>
aviezerl commented 2 years ago

The names are taken as strings from the column names in the input "anndata" file. How did you create it?

yonatans2 commented 2 years ago
wd <- "/net/mraid14/export/tgdata/users/yonshap/proj/mmcortex/"
setwd(wd)
library(anndata)
pmc <- readRDS('./data/pl_cort_peak_mc_smoothed_mg.rds')
rownames(pmc) <- gsub("-", "_", rownames(pmc))
peaks <- misha.ext::mat_to_intervs(pmc) %>% select(chrom, start, end)
ad <- anndata::AnnData(X = t(mcp), var = peaks)
anndata::write_h5ad(anndata = ad, filename = './data/test_anndata_atac_metacell.h5ad')
aviezerl commented 2 years ago

So set the colnames of pmc before creating it:

colnames(pmc) <- 1:ncol(pmc)
ad <- anndata::AnnData(X = t(mcp), var = peaks)
anndata::write_h5ad(anndata = ad, filename = './data/test_anndata_atac_metacell.h5ad')
yonatans2 commented 2 years ago

I did as much, but after import.........................................................................