seurat icon indicating copy to clipboard operation
seurat copied to clipboard

In SingleExIPlot(type = type, data = data[, x, drop = FALSE], idents = idents, : All cells have the same value of blacklist_ratio.

Open GaF123 opened this issue 6 months ago • 0 comments

R version 4.3.3 (2024-02-29) -- "Angel Food Cake" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-conda-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R.

load("/vast/palmer/pi/fan/fg349/analysis/mouse_ATAC/all_cells/my_environment.RData") cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr1-186356415-186437897"
  • ) Error in CoveragePlot(object = combined, group.by = "dataset", region = "chr1-186356415-186437897") : could not find function "CoveragePlot"

#library(spatstat.core) library(Signac) library(Seurat) Loading required package: SeuratObject Loading required package: sp

Attaching package: ‘SeuratObject’

The following objects are masked from ‘package:base’:

intersect, t

library(GenomicRanges) Loading required package: stats4 Loading required package: BiocGenerics

Attaching package: ‘BiocGenerics’

The following object is masked from ‘package:SeuratObject’:

intersect

The following objects are masked from ‘package:stats’:

IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

anyDuplicated, aperm, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated,
eval, evalq, Filter, Find, get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply, setdiff,
sort, table, tapply, union, unique, unsplit, which.max, which.min

Loading required package: S4Vectors

Attaching package: ‘S4Vectors’

The following object is masked from ‘package:utils’:

findMatches

The following objects are masked from ‘package:base’:

expand.grid, I, unname

Loading required package: IRanges

Attaching package: ‘IRanges’

The following object is masked from ‘package:sp’:

%over%

Loading required package: GenomeInfoDb

library(future) library(harmony) Error in library(harmony) : there is no package called ‘harmony’ library(ggplot2) library(EnsDb.Hsapiens.v86) Error in library(EnsDb.Hsapiens.v86) : there is no package called ‘EnsDb.Hsapiens.v86’ library(EnsDb.Mmusculus.v79) Loading required package: ensembldb Loading required package: GenomicFeatures Loading required package: AnnotationDbi Loading required package: Biobase Welcome to Bioconductor

Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see
'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: AnnotationFilter

Attaching package: ‘AnnotationFilter’

The following object is masked from ‘package:future’:

value

Attaching package: 'ensembldb'

The following object is masked from 'package:stats':

filter

library(openxlsx) Error in library(openxlsx) : there is no package called 'openxlsx' cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr1-186356415-186437897"
  • )

ggsave("cov_plot_tgfb2.pdf", plot = cov_plot, device = "pdf", width = 8, height = 6) cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr1-186356415-186437897"
  • )

ggsave("cov_plot_tgfb2.pdf", plot = cov_plot, device = "pdf", width = 8, height = 6) setwd("/vast/palmer/pi/fan/fg349/analysis/mouse_ATAC/all_cells") cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr1-186356415-186437897"
  • )

ggsave("cov_plot_tgfb2.pdf", plot = cov_plot, device = "pdf", width = 8, height = 6) cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr4-47353222-47414931"
  • )

ggsave("cov_plot_tgfbr1.pdf", plot = cov_plot, device = "pdf", width = 8, height = 6) cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr9-115917997-116004428"
  • )

ggsave("cov_plot_tgfbr2.pdf", plot = cov_plot, device = "pdf", width = 8, height = 6) cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr11-106545044-106606220"
  • )

ggsave("cov_plot_pecam1.pdf", plot = cov_plot, device = "pdf", width = 8, height = 6) cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr19-34218490-34232990"
  • )

ggsave("cov_plot_acta2.pdf", plot = cov_plot, device = "pdf", width = 8, height = 6) cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr1-71624683-71692334"
  • )

ggsave("cov_plot_FN1.pdf", plot = cov_plot, device = "pdf", width = 8, height = 6) cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr2-13578763-13587637"
  • )

ggsave("cov_plot_VIM.pdf", plot = cov_plot, device = "pdf", width = 8, height = 6) cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr10-24471340-24474576"
  • )

ggsave("cov_plot_ctgf.pdf", plot = cov_plot, device = "pdf", width = 8, height = 6) cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr5-138962726-138983125"
  • )

ggsave("cov_plot_pdgfa.pdf", plot = cov_plot, device = "pdf", width = 8, height = 6) cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr15-79880075-79899009"
  • )

ggsave("cov_plot_pdgfb.pdf", plot = cov_plot, device = "pdf", width = 8, height = 6) cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr5-30218112-30224973"
  • )

ggsave("cov_plot_il6.pdf", plot = cov_plot, device = "pdf", width = 8, height = 6) cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr9-20927281-20940113"
  • )

ggsave("cov_plot_icam1.pdf", plot = cov_plot, device = "pdf", width = 8, height = 6) cov_plot <- CoveragePlot(

  • object = combined,
  • group.by = 'dataset',
  • region = "chr11-81926397-81928279"
  • )

ggsave("cov_plot_ccl2.pdf", plot = cov_plot, device = "pdf", width = 8, height = 6) combined An object of class Seurat 166563 features across 35953 samples within 1 assay Active assay: ATAC (166563 features, 166563 variable features) 2 layers present: counts, data 2 dimensional reductions calculated: lsi, umap #mice:

extract gene annotations from EnsDb

annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Mmusculus.v79) |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01m 20s There were 21 warnings (use warnings() to see them)

change to UCSC style since the data was mapped to hg19

seqlevels(annotations) <- paste0('chr', seqlevels(annotations)) genome(annotations) <- "mm10"

add the gene information to the object

Annotation(combined) <- annotations #Computing QC Metrics combined <- NucleosomeSignal(object = combined) Found 4713 cell barcodes Done Processing 91 million linesFound 8364 cell barcodes Done Processing 123 million linesFound 7776 cell barcodes Done Processing 116 million linesFound 15100 cell barcodes Done Processing 179 million lines> View(combined) View(combined) View(combined) [email protected][["blacklist_region_fragments"]] [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [60] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [119] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [178] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [237] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [296] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [355] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [414] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [473] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [532] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [591] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [650] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [709] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [768] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [827] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [886] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [945] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [ reached getOption("max.print") -- omitted 34953 entries ] combined$nucleosome_group <- ifelse(combined$nucleosome_signal > 4, 'NS > 4', 'NS < 4') FragmentHistogram(object = combined, group.by = 'nucleosome_group', region = 'chr1-1-10000000') Warning messages: 1: Removed 78 rows containing non-finite outside the scale range (stat_bin()). 2: Removed 2 rows containing missing values or values outside the scale range (geom_bar()). combined$nucleosome_group <- ifelse(combined$nucleosome_signal > 4) Error in ifelse(combined$nucleosome_signal > 4) : argument "no" is missing, with no default

Create the plot and store it in a variable

frag_hist_plot <- FragmentHistogram(

  • object = combined,
  • group.by = 'nucleosome_group',
  • region = 'chr1-1-10000000'
  • )

Save the plot as a PDF

ggsave("frag_hist_plot.pdf", plot = frag_hist_plot, device = "pdf", width = 8, height = 6) Warning messages: 1: Removed 78 rows containing non-finite outside the scale range (stat_bin()). 2: Removed 2 rows containing missing values or values outside the scale range (geom_bar()). table(combined$nucleosome_group)

NS < 4 35953

summary(combined$nucleosome_signal) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00125 0.20547 0.26922 0.27224 0.33745 2.57838 combined An object of class Seurat 166563 features across 35953 samples within 1 assay Active assay: ATAC (166563 features, 166563 variable features) 2 layers present: counts, data 2 dimensional reductions calculated: lsi, umap combined <- TSSEnrichment(combined, fast = FALSE) Extracting TSS positions Finding + strand cut sites Finding - strand cut sites Computing mean insertion frequency in flanking regions Normalizing TSS score combined$high.tss <- ifelse(combined$TSS.enrichment > 2, 'High', 'Low') TSSPlot(combined, group.by = 'high.tss') + NoLegend()

Create the plot

p <- TSSPlot(combined, group.by = 'high.tss') + NoLegend()

Save the plot to a PDF file

ggsave("TSSPlot.pdf", plot = p, dpi = 600, height = 8, width = 8) combined <- TSSEnrichment(combined, fast = FALSE) Extracting TSS positions Finding + strand cut sites Finding - strand cut sites Computing mean insertion frequency in flanking regions Normalizing TSS score

combined$high.tss <- ifelse(combined$TSS.enrichment > 2, 'High', 'Low') TSSPlot(combined, group.by = 'high.tss') + NoLegend()

Create the plot

p <- TSSPlot(combined, group.by = 'high.tss') + NoLegend()

Save the plot to a PDF file

ggsave("TSSPlot.pdf", plot = p, dpi = 600, height = 8, width = 8) View(C1039G) combined$pct_reads_in_peaks <- combined$peak_region_fragments / combined$passed_filters * 100 combined$blacklist_ratio <- combined$blacklist_region_fragments / combined$peak_region_fragments p <- VlnPlot(

  • object = combined,
  • features = c('pct_reads_in_peaks', 'peak_region_fragments',
  •            'TSS.enrichment', 'blacklist_ratio', 'nucleosome_signal'),
    
  • pt.size = 0.1,
  • ncol = 5
  • ) **Warning message: In SingleExIPlot(type = type, data = data[, x, drop = FALSE], idents = idents, : All cells have the same value of blacklist_ratio.

This issue even happens when I use the mouse brain example dataset.

And I checked the data:

table(combined$blacklist_ratio)

0 

35953

summary(combined$blacklist_ratio) Min. 1st Qu. Median Mean 3rd Qu. Max. 0 0 0 0 0 0

GaF123 avatar Aug 14 '24 23:08 GaF123