tanaylab / mcatac Goto Github PK

View Code? Open in Web Editor NEW

0.0 0.0 0.0 24.43 MB

Metacell analysis for ATAC data

Home Page: https://tanaylab.github.io/mcATAC/

License: Other

R 99.57% Shell 0.43%

mcatac's People

Contributors

Watchers

mcatac's Issues

generate_pheatmap_annotation returns error when fed with a cell type annotation vector

This should work (mcmd$st is a cell-type annotation):
ann_list = generate_pheatmap_annotation(clust_vec = mcmd$st, feature_type = "cell_type", feature_annotation = "color")

But gives an error:

`
Error in stop_vctrs(class = c(class, "vctrs_error_names"), ...): Names must be unique.
✖ These names are duplicated:

"color" at locations 1 and 2.
Traceback:

generate_pheatmap_annotation(clust_vec = mcmd$st, feature_type = "cell_type",
. feature_annotation = "color")
tibble(name = cts, color = chameleon::distinct_colors(length(cts))$name) %>%
. rename(:=(!!feature_annotation, name)) # at line 67-68 of file /home/feshap/src/mcATAC/R/utils.R
rename(., :=(!!feature_annotation, name))
rename.data.frame(., :=(!!feature_annotation, name))
tidyselect::eval_rename(expr(c(...)), .data)
rename_impl(data, names(data), as_quosure(expr, env), strict = strict,
. name_spec = name_spec)
with_subscript_errors(vctrs::vec_as_names(names, repair = "check_unique"))
tryCatch(instrument_base_errors(expr), vctrs_error_subscript = function(cnd) {
. cnd$subscript_action <- subscript_action(type)
. cnd$subscript_elt <- "column"
. cnd_signal(cnd)
. })
tryCatchList(expr, classes, parentenv, handlers)
tryCatchOne(expr, names, parentenv, handlers[[1L]])
doTryCatch(return(expr), name, parentenv, handler)
instrument_base_errors(expr)
withCallingHandlers(expr, simpleError = function(cnd) {
. abort(conditionMessage(cnd), parent = cnd)
. })
vctrs::vec_as_names(names, repair = "check_unique")
(function ()
. validate_unique(names = names, arg = arg))()
validate_unique(names = names, arg = arg)
stop_names_must_be_unique(names, arg)
stop_names(class = "vctrs_error_names_must_be_unique", arg = arg,
. names = names)
stop_vctrs(class = c(class, "vctrs_error_names"), ...)
abort(message, class = c(class, "vctrs_error"), ...)
signal_abort(cnd, .file)`

mcc_to_tracks

> devtools::load_all("~/src/mcATAC/")
ℹ Loading mcATAC
ℹ Parallelization enabled. Using 77 threads.
> gset_genome("mm10")
> mcc <- mcc_read("/net/mraid14/export/tgdata/users/yonshap/proj/embenh/data/mcc_lg/")
✔ Succesfully read a ScCounts object from /net/mraid14/export/tgdata/users/yonshap/proj/embenh/data/mcc_lg/
✔ Succesfully read a McCounts object from /net/mraid14/export/tgdata/users/yonshap/proj/embenh/data/mcc_lg/
> mcc_to_tracks(mc_counts=mcc, track_prefix='test_gastru_large_tracks', overwrite = T, normalize = T)
→ Normalizing each metacell by its total counts
Error in do.ply(i) :
  task 1 failed - "error in evaluating the argument 'x' in selecting a method for function 't': Cholmod error 'problem too large' at file ../Core/cholmod_dense.c, line 102"

McPeaks export to h5ad remove slots

This is the slots of the object:

@mat: a numeric matrix where rows are peaks and columns are metacells.
Can be a sparse matrix.
@peaks: a misha intervals set with the peak definitions.
@genome: genome assembly of the peaks
@egc: a numeric matrix which contains normalized metacell accessibility.
@fp: a matrix showing for each peak (row) the relative enrichment of umis
in log2 scale.
@metadata: a tibble with a column called 'metacell' and additional
metacell annotations.

But after saving it to h5ad we don't have the egc and fp. From the code:
slots <- slots[slots %!in% c("egc", "fp", "mat", "peaks", "genome", "metadata", "ignore_peaks", "ignore_pmat", "rna_egc", "tad_based")]

Now that we manipulate the matrix for normalization we probably want to save all of them to re-load them, both in R and in python. We should just use the general layers of the anndata to do this

solve unnamed peaks set without raising Error when creating McPeaks object from McCounts

if a valid peaks interval set is passed (i.e. has 'chrom','start','end'), but lacking peak_name column - call peak_names() internally to solve this issue, instead of aborting:

peak_set = fread('path/to/peak/file.txt')

mcatac_obj <- mcc_to_mcatac(mcc, peaks = peak_set, metadata = mc_md)
Error in `mcc_to_mcatac()`:
! The `peaks` must have a column called peak_name

Solution:
if(!has_name(peak_set)){ peak_set$peak_name = peak_names(peak_set) }

plot_atac_rna_markers looks for colors in cell_type column

(or vice versa)

Are metacell numbers 0- or 1-based

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"
>

generate_motif_pssm_matrix

Suggestions:

Add original peak name to output
Warn user than output intervals are not necessarily the same as input, unless the input had uniform peak width and it was specified in the relevant parameter

mcc_to_tracks promts(...) prompts Error message even when overwrite=T

small comment.
when re-creating tracks to the same directory and passing overwrite=True, it complains about an already-existing marginal track.
setting create_marginal_track=False, still results in an Error message:
"Error: File .../marginal already exists"

mcc_to_mcatac

Run this:

devtools::load_all("~/src/mcATAC")
wd <- "/net/mraid14/export/tgdata/users/yonshap/proj/mmcortex"
setwd(wd)
mcmd <- readr::read_tsv('./BonevCollab/mcmd_pl_cort.tsv')
mcc <- mcc_read('./data/mmcortex_mcc/', id = 'mmcortex')
pks_all <- readRDS('./data/all_peaks_for_mmcortex_atac_mc.rds')
atac_mc <- mcc_to_mcatac(mc_counts = mcc, peaks = pks_all, metadata = dplyr::rename(mcmd, metacell = mc, cell_type = st))

Output:

> devtools::load_all("~/src/mcATAC")
Loading mcATAC
ℹ Parallelization enabled. Using 77 threads.
> wd <- "/net/mraid14/export/tgdata/users/yonshap/proj/mmcortex"
> setwd(wd)
> mcmd <- readr::read_tsv('./BonevCollab/mcmd_pl_cort.tsv')
Rows: 454 Columns: 11
── Column specification ───────────────────────────────────────────────────────────────────────────────────[...]
> mcc <- mcc_read('./data/mmcortex_mcc/', id = 'mmcortex')
|--------------------------------------------------|
[...]
=======|
✔ Succesfully read a McCounts object from ./data/mmcortex_mcc/
> pks_all <- readRDS('./data/all_peaks_for_mmcortex_atac_mc.rds')
> atac_mc <- mcc_to_mcatac(mc_counts = mcc, peaks = pks_all, metadata = dplyr::rename(mcmd, metacell = mc, cell_type = st))
97%...100%
81%...100%
94%...100%
87%...100%
96%...100%
59%...75%...100%
100%
Error in do.ply(i) :
  task 1 failed - "error in evaluating the argument 'x' in selecting a method for function 't': error in evaluating the argument 'x' in selecting a method for function 't': contrasts can be applied only to factors with 2 or more levels"
>

plot_tracks_at_locus issue

In plot_tracks_at_locus, if you choose to plot an interval and not a gene TSS extended left and right, you lose the option to plot gene expression. There should be an option to plot a specified interval (like a TAD) and still show gene expression on the side.

many functions resulted in Error because of atac_mc@mat is 'dgCMatrix'

I'm not sure why but when I built an atac_mc object using mcc_to_mcatac many functions (such as: rna_atac_cor_knn, plot_atac_atac_cor, etc...) doesn't handle the input properly..

example:
my_mc_atac <- readRDS('/home/tomgo/raid/proj/mtec_populations/analysis_dir/multiome_mc_model/2022/atac_analysis/objs/mc_atac_obj.Rds')
plot_atac_atac_cor(atac_mc = my_mc_atac)
Error: "x" argument must be a matrix of numeric values

`atac_ignore_peaks` fails when given peak names

This breaks filter_features (and its tests):

atac_sc <- import_from_10x("pbmc_data", genome = "hg38", id = "pbmc", description = "PBMC from a healthy donor - granulocytes removed through cell sorting (10k)")
atac_sc <- filter_features(atac_sc, minimal_max_umi = 3, min_peak_length = 200, max_peak_length = 1000, max_peak_density = 250)
#> Error in `name_enhancers()` at mcATAC/R/PeakIntervals.R:86:12:
#> ! Class of `atac` is not recognized (should be either ScATAC, McATAC or PeakIntervals

write_sc_counts_from_fragments -> scc_read

> write_sc_counts_from_fragments("./raid/proj/matching/data/pbmc_granulocyte_sorted_10k_atac_fragments.tsv", out_dir="./raid/proj/matching/data/pbmc_fragment_reads", genome = "hg38", overwrite=T, cell_na→ Writing to ./raid/proj/matching/data/pbmc_fragment_reads
ℹ 'tabix' was not found or an index file doesn't exist. Using all chromosomes
ℹ Processing 504 genomic bins of maximal size 50000000
→ Processing "chr1_0_50000000"
→ Processing "chr1_50000000_100000000"
→ Processing "chr1_100000000_150000000"
→ Processing "chr1_150000000_200000000"
→ Processing "chr1_200000000_248956422"
→ Processing "chr1_GL383518v1_alt_0_182439"
→ Processing "chr1_GL383519v1_alt_0_110268"
→ Processing "chr1_GL383520v2_alt_0_366580"
→ Processing "chr1_KI270706v1_random_0_175055"
→ Processing "chr1_KI270707v1_random_0_32032"
→ Processing "chr1_KI270708v1_random_0_127682"
gzip: ./raid/proj/matching/data/pbmc_granulocyte_sorted_10k_atac_fragments.tsv: not in gzip format
→ Processing "chr1_KI270709v1_random_0_66860"
gzip: ./raid/proj/matching/data/pbmc_granulocyte_sorted_10k_atac_fragments.tsv: not in gzip format
→ Processing "chr1_KI270710v1_random_0_40176"
gzip: ./raid/proj/matching/data/pbmc_granulocyte_sorted_10k_atac_fragments.tsv: not in gzip format

...
...

✔ Finished processing "chrUn_KI270753v1_0_62944"
✔ Finished processing "chrX_0_50000000"
✔ Finished processing "chrUn_KI270756v1_0_79590"
✔ Finished processing "chrUn_KI270752v1_0_27745"
✔ Finished processing "chrX_KI270913v1_alt_0_274009"
✔ Finished processing "chrUn_KI270754v1_0_40191"
✔ Finished processing "chrX_KI270880v1_alt_0_284869"
✔ Finished processing "chrUn_KI270751v1_0_150742"
✔ Finished processing "chrUn_KI270757v1_0_71251"
✔ Finished processing "chrY_0_50000000"
✔ Finished processing "chrX_KI270881v1_alt_0_144206"
✔ Finished processing "chrX_50000000_100000000"
✔ Finished processing "chrX_100000000_150000000"
✔ Finished processing "chrY_KI270740v1_random_0_37240"
✔ Finished processing "chrY_50000000_57227415"
✔ Finished processing "chrX_150000000_156040895"
ℹ Created sparse matrices for 504 genomic bins
✔ Created an ScCounts object at ./raid/proj/matching/data/pbmc_fragment_reads
> scc <- scc_read("./raid/proj/matching/data/pbmc_fragment_reads/")
Error in do.ply(i) :
  task 1 failed - "skip=2 but the input only has 2 lines"
In addition: Warning messages:
1: Problem with `mutate()` column `start`.
ℹ `start = as.numeric(start)`.
ℹ NAs introduced by coercion
2: The following bins were not found in the genome: "chr1_GL383518v1_alt_0_182439", "chr1_GL383519v1_alt_0_110268", "chr1_GL383520v2_alt_0_366580", "chr1_KI270706v1_random_0_175055",
"chr1_KI270707v1_random_0_32032", "chr1_KI270708v1_random_0_127682", "chr1_KI270709v1_random_0_66860", "chr1_KI270710v1_random_0_40176", "chr1_KI270711v1_random_0_42210",
"chr1_KI270712v1_random_0_176043", "chr1_KI270713v1_random_0_40745", "chr1_KI270714v1_random_0_41717", "chr1_KI270759v1_alt_0_425601", "chr1_KI270760v1_alt_0_109528", "chr1_KI270761v1_alt_0_165834",
"chr1_KI270762v1_alt_0_354444", "chr1_KI270763v1_alt_0_911658", "chr1_KI270764v1_alt_0_50258", "chr1_KI270765v1_alt_0_185285", "chr1_KI270766v1_alt_0_256271", "chr1_KI270892v1_alt_0_162212",
"chr10_GL383545v1_alt_0_179254", "chr10_GL383546v1_alt_0_309802", "chr10_KI270824v1_alt_0_181496", "chr10_KI270825v1_alt_0_188315", "chr11_GL383547v1_alt_0_154407", "chr11_JH159136v1_alt_0_200998",
"chr11_JH159137v1_alt_0_191409", "chr11_KI270721v1_random_0_100316", "chr11_KI270826v1_alt_0_186169", "chr11_KI270827v1_alt_0_67707", "chr11_KI270829v1_alt_0_204059", "chr11_KI270830v1_alt_0_177092",
"chr11_KI270831v1_alt_0_296895", "chr11_KI270832v1_alt_0_210133", "chr11_KI270902v1_alt_0_106711", "chr11_KI270903v1_alt_0_214625", "chr11_KI270927v1_alt_0_218612", "chr12_GL383549v1_alt_0_120804",
"chr12_GL383550v2_alt_0_169178", "chr12_GL383551v1_alt_0_184319", "chr12_GL383552v1_alt_0_138655", "chr12_GL383553v2_alt_0_152874", "chr12_GL877875v1_alt_0_167313", "chr12_GL877876v1_alt_0_408271",
"chr12_KI270833v1_alt_0_76061", "chr12_KI270834v1_alt_0_119498", "chr12_KI270835v1_alt_0_238139", "chr12_KI270836v1_alt_0_56134", "chr12_KI270837v1_alt_0_40090", "chr12_KI270904v1_alt_0_572349",
"chr13_KI270838v1_alt_0_306913", "chr13_KI270839v1_alt_0_180306", "chr13_KI270840v1_alt_0_191684", "chr13_KI270841v1_alt_0_169134", "chr13_KI270842v1_alt_0_37287", "chr13_KI270843v1_alt_0_103832",
"chr14_GL000009v2_random_0_201709", "chr14_GL000194v1_random_0_191469", "chr14_GL000225v1_random_0_211173", "chr14_KI270722v1_random_0_194050", "chr14_KI270723v1_random_0_38115",
"chr14_KI270724v1_random_0_39555", "chr14_KI270725v1_random_0_172810", "chr14_KI270726v1_random_0_43739", "chr14_KI270844v1_alt_0_322166", "chr14_KI270845v1_alt_0_180703",
"chr14_KI270846v1_alt_0_1351393", "chr14_KI270847v1_alt_0_1511111", "chr15_GL383554v1_alt_0_296527", "chr15_GL383555v2_alt_0_388773", "chr15_KI270727v1_random_0_448248", "chr15_KI270848v1_alt_0_327382",
"chr15_KI270849v1_alt_0_244917", "chr15_KI270850v1_alt_0_430880", "chr15_KI270851v1_alt_0_263054", "chr15_KI270852v1_alt_0_478999", "chr15_KI270905v1_alt_0_5161414", "chr15_KI270906v1_alt_0_196384",
"chr16_GL383556v1_alt_0_192462", "chr16_GL383557v1_alt_0_89672", "chr16_KI270728v1_random_0_1872759", "chr16_KI270853v1_alt_0_2659700", "chr16_KI270854v1_alt_0_134193", "chr16_KI270855v1_alt_0_232857",
"chr16_KI270856v1_alt_0_63982", "chr17_GL000205v2_random_0_185591", "chr17_GL000258v2_alt_0_1821992", "chr17_GL383563v3_alt_0_375691", "chr17_GL383564v2_alt_0_133151", "chr17_GL383565v1_alt_0_223995",
"chr17_GL383566v1_alt_0_90219", "chr17_JH159146v1_alt_0_278131", "chr17_JH159147v1_alt_0_70345", "chr17_JH159148v1_alt_0_88070", "chr17_KI270729v1_random_0_280839", "chr17_KI270730v1_random_0_112551",
"chr17_KI270857v1_alt_0_2877074", "chr17_KI270858v1_alt_0_235827", "chr17_KI270859v1_alt_0_108763", …. This usually means that the genome ("hg38") is not compatible with the bins.

scc_to_peaks malfunction?

library(mcATAC)
wd <- '/net/mraid14/export/tgdata/users/yonshap/proj/mmcortex'
setwd(wd)
my_genome <- "mm10"
gset_genome(my_genome)
mca <- readRDS('./output/mcatac/mmcortex_mcatac_feat_peaks.rds')
scc <- scc_read(path = './data/frag_reads_28122022/')

 ✔ Succesfully read a ScCounts object from ./data/frag_reads_28122022/

>  scpeaks <- scc_to_peaks(sc_counts = scc, peaks = dplyr::select(mca@peaks, chrom, start, end, peak_name))
Error in data.frame(i = all_0_idxs, j = 1, x = 0) :
  arguments imply differing number of rows: 0, 1

mct_plot_region running out of memory

> devtools::load_all("~/src/mcATAC/")
ℹ Parallelization enabled. Using 77 threads.
> setwd('/net/mraid14/export/tgdata/users/yonshap/proj/mmcortex')
> my_genome <- "mm10"
> gset_genome(my_genome)
> options(gmax.data.size = 1e+9)
> mctracks <- readRDS('./data/mmcortex_mctracks.rds')
> exons <- gintervals.load('intervs.global.exon')
> gene <- 'Satb2'
> gene_exons <- dplyr::filter(exons, geneSymbol == gene)
> gene_intervals <- dplyr::summarise(gene_exons, chrom = unique(chrom), start = min(start), end = max(end))
> mct_plot_region(mct = mctracks, intervals = gene_intervals, downsample_n = 1e+5)
→ Extracting region chr1:56793980-56971340
Error: Failed to allocate shared memory: Cannot allocate memory
Memory usage of the library can be controlled via gmax.data.size option (see options, getOptions).

mcc_to_tracks

> setwd("/net/mraid14/export/tgdata/users/yonshap/proj/mmcortex")
> devtools::load_all("~/src/mcATAC")
ℹ Loading mcATAC
ℹ Parallelization enabled. Using 77 threads.
> gset_genome("mm10")
> mcc <- mcc_read('./data/frag_read_mcc/')
> set_parallel(20)
ℹ Parallelization enabled. Using 20 threads.
> mcc_to_tracks(mcc, "mmcortex_microcluster", create_marginal_track=FALSE, overwrite=T)
ℹ Creating tracks for 192 metacells
ℹ Smoothing over 201 bp window
ℹ Tracks resolution: 10 bp
→ Extracting per-metacell data
→ Creating mmcortex_microcluster.mc1 track
→ Creating mmcortex_microcluster.mc2 track
→ Creating mmcortex_microcluster.mc3 track
→ Creating mmcortex_microcluster.mc4 track
→ Creating mmcortex_microcluster.mc5 track
…
77%...→ Creating mmcortex_microcluster.mc180 track
72%...72%...77%...90%...77%...100%
→ Creating mmcortex_microcluster.mc173 track
90%...90%...100%
100%
4%...13%...4%...22%...13%...31%...22%...40%...31%...50%...40%...59%...50%...63%...59%...72%...63%...77%...72%...90%...100%
77%...90%...100%
Error in do.ply(i) :
  task 21 failed - "Problem with `filter()` input `..1`.
ℹ Input `..1` is `metacell == "21"`.
✖ cannot allocate vector of size 5.1 Gb"
In addition: Warning message:
In mclapply(argsList, FUN, mc.preschedule = preschedule, mc.set.seed = set.seed,  :
  scheduled cores 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 15, 16, 17, 19 did not deliver results, all values of the jobs will be affected

@param `atac` for get_peak_coverage_stats?

different functions under feature_selection.R don't accept mcATAC objects as input, and it's not clear what is the correct input for them.
(e.g., get_peaks_coverage_stats / plot_peak_coverage_distribution ,etc...)

documentation says it should be "an ScATAC or McATAC object", but mcATAC is raising error

tanaylab / mcatac Goto Github PK

mcatac's People

Contributors

Watchers

mcatac's Issues

Recommend Projects

Recommend Topics

Recommend Org