seurat
seurat copied to clipboard
In SingleExIPlot(type = type, data = data[, x, drop = FALSE], idents = idents, : All cells have the same value of blacklist_ratio.
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 defaultCreate 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