Hi @timoast , I was trying to produce the region matrix from my Seurat object, but I faced this error:
hub_df_cluster_3_module_4_TSS_gr
GRanges object with 100 ranges and 6 metadata columns:
seqnames ranges strand | gene_name module kME id GeneID region
<Rle> <IRanges> <Rle> | <character> <character> <numeric> <character> <character> <character>
[1] chr2 169632646-169634646 + | Tshz2 Cluster_tre4 0.667834 ENSMUSG00000047907 228911 chr2-169632646-16963..
[2] chr11 44617134-44619134 + | Ebf1 Cluster_tre4 0.609192 ENSMUSG00000057098 13591 chr11-44617134-44619..
[3] chr5 136882107-136884107 - | Col26a1 Cluster_tre4 0.608089 ENSMUSG00000004415 140709 chr5-136882107-13688..
[4] chr14 67232292-67234292 + | Ebf2 Cluster_tre4 0.598051 ENSMUSG00000022053 13592 chr14-67232292-67234..
[5] chr3 141464564-141466564 + | Unc5c Cluster_tre4 0.530794 ENSMUSG00000059921 22253 chr3-141464564-14146..
... ... ... ... . ... ... ... ... ... ...
[96] chr15 93594891-93596891 - | Prickle1 Cluster_tre4 0.288270 ENSMUSG00000036158 106042 chr15-93594891-93596..
[97] chr7 51878145-51880145 + | Gas2 Cluster_tre4 0.287304 ENSMUSG00000030498 14453 chr7-51878145-51880145
[98] chr13 88820472-88822472 + | Edil3 Cluster_tre4 0.286328 ENSMUSG00000034488 13612 chr13-88820472-88822..
[99] chr19 14596983-14598983 - | Tle4 Cluster_tre4 0.286068 ENSMUSG00000024642 21888 chr19-14596983-14598..
[100] chr1 21960942-21962942 - | Kcnq5 Cluster_tre4 0.283888 ENSMUSG00000028033 226922 chr1-21960942-21962942
-------
seqinfo: 20 sequences from an unspecified genome; no seqlengths
obj
An object of class Seurat
169476 features across 9995 samples within 3 assays
Active assay: ATAC (134191 features, 0 variable features)
2 other assays present: integrated_rna, RNA
8 dimensional reductions calculated: pca_rna, umap_rna, integrated_lsi, umap_atac, umap_wnn, UMAP, pca, umap
Idents(obj) <- "treatment_timepoint_cluster"
prova = RegionMatrix(obj, regions = hub_df_cluster_3_module_4_TSS_gr, assay="ATAC", group.by = "treatment_timepoint_cluster", idents = c("3_KO_d8", "3_WT_d8", "3_KO_d6", "3_WT_d6"), key = "picchi", upstream = 1000, downstream = 1000, verbose = TRUE)
[W::hts_idx_load3] The index file is older than the data file: /media/olga/14F0454F24001A5C/single_cell/Fragment/KO_d6_r1_fragments.tsv.gz.tbi
Errore in res$start - start(x = plus.strand[j]) :
non-numeric argument passed to a binary operator
I used the develop branch of Signac:
if (!requireNamespace("remotes", quietly = TRUE))
install.packages("remotes")
remotes::install_github("stuart-lab/signac", ref = "develop")
EDIT::
If I use the stable version install.packages("Signac") version 1.9.0, when I was trying to create the region matrix it seems to work fine until I see that some regions in the region matrix there were 0 counts, although I can see that there are fragments in the coverageplot, as I saw in this thread:
#1076
1) extract one region matrix stored in "picchi" key
c=prova@assays$ATAC@positionEnrichment$picchi$`3_KO_d8`
rownames(c) <- hub_df_cluster_3_module_4_TSS$gene_name
> dim(c)
[1] 100 2001
2) calculate which regions have 0 counts
c_df= as.data.frame(c)
c_df$tot = rowSums(c_df)
c_zero= as.data.frame(c_df[c_df$tot == 0, ])
c_zero$name= rownames(c_zero)
dim(c_zero)
[1] 55 2003
3) calculate which regions have more than a tot of zero counts to be plotted in the heatmap
c_mat = c_df[c_df$tot != 0, ]
c_mat= as.matrix(c_mat[,-2002])
dim(c_mat)
[1] 45 2001
4) plot a coverage of some regions with zero counts
rownames(c_zero)
[1] "Cdh11" "Pbx1" "Grid2" "Tmtc1" "Emilin1" "Lima1" "Klf12" "Slc2a13"
[9] "Mfap4" "Astn2" "Mrc2" "Eda" "Adgrl3" "Trps1" "Stat3" "Pdgfrb"
[17] "Ptprd" "Lhfpl2" "Col25a1" "Pcolce" "Col13a1" "Tmem131l" "Lpp" "Junb"
[25] "Ghr" "Tbx1" "Sncaip" "Fbn1" "Dock1" "Gab2" "Ppm1l" "Zeb1"
[33] "Lix1" "Nrg1" "Egr1" "Arhgap28" "Basp1" "Hmcn1" "Sobp" "Colec12"
[41] "Kalrn" "Fosb" "Foxn3" "Boc" "2610307P16Rik" "Arhgap26" "Ror1" "4930438E09Rik"
[49] "Cntln" "Six1" "Prickle1" "Gas2" "Edil3" "Tle4" "Kcnq5"
DefaultAssay(prova) <-"ATAC"
cI= CoveragePlot(
object = try,
region = hub_df_cluster_3_module_4_TSS[hub_df_cluster_3_module_4_TSS$gene_name == c_zero$name[2],]$region,
extend.upstream = 0,
extend.downstream = 0,
group.by = "treatment_timepoint_cluster",
idents = c("3_KO_d8", "3_WT_d8"),
features = hub_df_cluster_3_module_4_TSS[hub_df_cluster_3_module_4_TSS$gene_name == c_zero$name[2],]$gene_name,
peaks = T,
annotation = T)
cI
can u help me?
RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4
