Region Matrix won't work #1076
Description
Hi Tim, Thanks in advance...
Essentially I cannot get region matrix to work. Here is the code I used to generate my Seurat object before trying to use feature matrix:
annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Mmusculus.v79)
seqlevelsStyle(annotation) <- "UCSC"
peaks.AF1 <- read.table(
file = "AF1_atac_peaks.bed",
col.names = c("chr", "start", "end")
)
#I then filtered by peak width. This dataset contains 4 objects but I am only showing one.
peakwidths <- width(combined.peaks)
combined.peaks <- combined.peaks[peakwidths < 10000 & peakwidths > 20]
combined.peaks
# Read metadata
md.AF1 <- read.table(
file = "AF1_per_barcode_metrics.csv",
stringsAsFactors = FALSE,
sep = ",",
header = TRUE,
row.names = 1
)[-2, ]
# Filter
md.AF1 <- md.AF1[md.AF1$atac_raw_reads > 250, ]
# create fragment objects
frags.AF1 <- CreateFragmentObject(
path = "AF1_atac_fragments.tsv.gz",
cells = rownames(md.AF1)
)
# Make counts using combined peak set
AF1.counts <- FeatureMatrix(
fragments = frags.AF1,
features = combined.peaks,
cells = rownames(md.AF1)
)
AF1_assay <- CreateChromatinAssay(AF1.counts, fragments = frags.AF1, annotation = annotation, min.features = -1)
AF1 <- CreateSeuratObject(AF1_assay, assay = "ATAC", meta.data=md.AF1)
AF1[["RNA"]]<- CreateAssayObject(counts=counts$'Gene Expression'[,colnames(AF1)])
AF1
Output:
An object of class Seurat
178393 features across 2283 samples within 2 assays
Active assay: ATAC (146108 features, 0 variable features)
1 other assay present: RNA
## I then performed qc, and merged 4 different seurats together.
## I went through the weighted nearest neighbors vignett to identify clusters, make reductions, ect.
## Then I try to do the region matrix to look at the signal within a few select clusters of cells:
# First make a granges object
bed <- read.table('filepath', sep = '\t', header = FALSE)
names(bed)[1] <- 'chromosome'
names(bed)[2] <- 'start'
names(bed)[3] <- 'end'
newrange <- makeGRangesFromDataFrame(bed, start.field = 'start', end.field = 'end',
seqinfo = c('chr1', 'chr2', 'chr3', 'chr4', 'chr5', 'chr6', 'chr7', 'chr8', 'chr9', 'chr10',
'chr11', 'chr12', 'chr13', 'chr14', 'chr15', 'chr16', 'chr17', 'chr18',
'chr19', 'chr20', 'chr21', 'chrX'), seqnames.field = 'chromosome' )
newrange
Output:
GRanges object with 1112 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr1 130803960-130804946 *
[2] chr1 9858933-9859881 *
[3] chr1 77048816-77049875 *
[4] chr1 79631214-79632159 *
[5] chr1 92343525-92344341 *
... ... ... ...
[1108] chr19 10433720-10434542 *
[1109] chr19 47561819-47562748 *
[1110] chr19 40368233-40369160 *
[1111] chr19 5821966-5822899 *
[1112] chrX 159942321-159943395 *
-------
seqinfo: 22 sequences from an unspecified genome; no seqlengths
## do regionmatrix:
idents.plot <- c("x1", "x2", "x3", 'x4', "x5", "x6")
DefaultAssay(Seurat)<- 'ATAC'
obj <-RegionMatrix(
Seurat,
newrange,
key = 'myregions',
assay = NULL,
idents = idents.plot,
upstream = 1000,
downstream = 1000,
verbose = TRUE,
)
obj
output:
An object of class Seurat
200391 features across 14569 samples within 4 assays
Active assay: ATAC (146056 features, 146056 variable features)
3 other assays present: RNA, chromvar, SCT
6 dimensional reductions calculated: pca, umap.rna, lsi, umap.atac, wnn.umap, umap
## Then I try RegionHeatMap
RegionHeatmap(
obj,
key = 'myregions',
assay = 'ATAC',
idents = idents.plot,
normalize = TRUE,
upstream = 1000,
downstream = 1000,
max.cutoff = "q95",
cols = NULL,
min.counts = 1,
window = (2000)/30,
order = TRUE,
nrow = NULL
)
This returns 0 counts for any of the identities, although I know that isn't the case because there are fragments in the coverageplot. An example:
CoveragePlot(Seurat, 'chr1-130803960-130804946', idents = idents.plot)