Weird peak plot after merging datasets

[Ping-lin14]

Hi timoast

After merging three pieces of data and extracting the data from the inspection fragment file, I found that the result was not as expected. As shown in the figure.

How can this be improved?

[timoast]

Hi @Ping-lin14, can you provide some more information. What was the expected result? What code did you run to get to this point?

[Ping-lin14]

I follow the teaching of merging objects. According to the teaching, the expected result is that there should be similar peaks. But my results show that there is a big difference in tumor 3234, at least I think.

library(Signac)
library(Seurat)
library(future)
library(ggplot2)
library(dplyr)
library(patchwork)
library(GenomeInfoDb)
library(GenomicRanges)
library(EnsDb.Hsapiens.v86)
set.seed(1234)

##Creating a common peak set

# read in peak sets
peaks_n <- read.table(
  file = "/home/pinglin/ATAC/S239N_3234N/peaks.bed",
  col.names = c("chr", "start", "end")
)

peaks_2391t <- read.table(
  file = "/home/pinglin/ATAC/s2391_T/peaks.bed",
  col.names = c("chr", "start", "end")
)

peaks_3234t <- read.table(
  file = "/home/pinglin/ATAC/S3234T/peaks.bed",
  col.names = c("chr", "start", "end")
)
# convert to genomic ranges
gr.thym_n <- makeGRangesFromDataFrame(peaks_n)
gr.thym_2391t <- makeGRangesFromDataFrame(peaks_2391t)
gr.thym_3234t <- makeGRangesFromDataFrame(peaks_3234t)

# Create a unified set of peaks to quantify in each dataset
combined.peaks <-  reduce(x = c(gr.thym_n, gr.thym_2391t, gr.thym_3234t))

# Filter out bad peaks based on length
peakwidths <- width(combined.peaks)
combined.peaks <- combined.peaks[peakwidths  < 10000 & peakwidths > 20]
combined.peaks

##Create Fragment objects
# load metadata
md.thym_n <- read.table(
  file = "/home/pinglin/ATAC/S239N_3234N/singlecell.csv",
  stringsAsFactors = FALSE,
  sep = ",",
  header = TRUE,
  row.names = 1
)[-1, ] # remove the first row

md.thym_2391t <- read.table(
  file = "/home/pinglin/ATAC/s2391_T/singlecell.csv",
  stringsAsFactors = FALSE,
  sep = ",",
  header = TRUE,
  row.names = 1
)[-1, ]

md.thym_3234t <- read.table(
  file = "/home/pinglin/ATAC/S3234T/singlecell.csv",
  stringsAsFactors = FALSE,
  sep = ",",
  header = TRUE,
  row.names = 1
)[-1, ]

md.thym_n <- md.thym_n[md.thym_n$passed_filters >1000, ]
md.thym_2391t <- md.thym_2391t[md.thym_2391t$passed_filters > 1000, ]
md.thym_3234t <- md.thym_3234t[md.thym_3234t$passed_filters > 1000, ]

# create fragment objects
frags.n <- CreateFragmentObject(
  path = "/home/pinglin/ATAC/S239N_3234N/fragments.tsv.gz",
  cells = rownames(md.thym_n)
)

frags.2391t <- CreateFragmentObject(
  path = "/home/pinglin/ATAC/s2391_T/fragments.tsv.gz",
  cells = rownames(md.thym_2391t)
)

frags.3234t <- CreateFragmentObject(
  path = "/home/pinglin/ATAC/S3234T/fragments.tsv.gz",
  cells = rownames(md.thym_3234t)
)
##Quantify peaks in each dataset
thym_n.counts <- FeatureMatrix(
  fragments = frags.n,
  features = combined.peaks,
  cells = rownames(md.thym_n)
)

thym_2391t.counts <- FeatureMatrix(
  fragments = frags.2391t,
  features = combined.peaks,
  cells = rownames(md.thym_2391t)
)

thym_3234t.counts <- FeatureMatrix(
  fragments = frags.3234t,
  features = combined.peaks,
  cells = rownames(md.thym_3234t)
)

thym_n_assay <- CreateChromatinAssay(thym_n.counts, fragments = frags.n, genome = 'hg38')
normal <- CreateSeuratObject(thym_n_assay, assay = "ATAC" , meta.data = md.thym_n)
#normal <- CreateSeuratObject(normal, assay = "peaks", meta.data = md.thtm_n)

thym_2391t_assay <- CreateChromatinAssay(thym_2391t.counts, fragments = frags.2391t, genome = 'hg38')
tumor_2391 <- CreateSeuratObject(thym_2391t_assay, assay = "ATAC" , meta.data = md.thym_2391t)
#tumor <- CreateSeuratObject(tumor, assay = "peaks", meta.data = md.thtm_t)

thym_3234t_assay <- CreateChromatinAssay(thym_3234t.counts, fragments = frags.3234t, genome = 'hg38')
tumor_3234 <- CreateSeuratObject(thym_3234t_assay, assay = "ATAC" , meta.data = md.thym_3234t)
#tumor <- CreateSeuratObject(tumor, assay = "peaks", meta.data = md.thtm_t)

##Merge objects
normal$dataset <- 'normal'
tumor_2391$dataset <- 'tumor_2391'
tumor_3234$dataset <- 'tumor_3234'

thym_merging <- merge(
  x = normal,
  y = list(tumor_2391, tumor_3234),
  add.cell.ids = c("N","2391", "3234T")
)

thym_merging[["ATAC"]]

saveRDS(thym_merging, file = "thym_2391_3234_hg38.RDS")

thym_merging <- RunTFIDF(thym_merging)
thym_merging <- FindTopFeatures(thym_merging, min.cutoff = 20)
thym_merging <- RunSVD(thym_merging)
thym_merging <- RunUMAP(thym_merging, dims = 2:50, reduction = 'lsi')

DimPlot(thym_merging, group.by = 'dataset', pt.size = 0.1)

CoveragePlot(
  object = thym_merging,
  group.by = 'dataset',
  region = "chr14-99700000-99760000"
)

[timoast]

The merge vignette just shows how to merge the objects, the results will vary depending on your data. It looks like two of the samples here have higher background. If you saw different results when looking at each dataset individually, that could indicate this is a problem with the software.

[Ping-lin14]

After viewing each data separately, normal and tumor 2391 have higher background values.
normal

tumor 2391

tumor3234

[timoast]

This appears to be an issue with the data, rather than any issue with the code or software

[Ping-lin14]

So you suggest that I continue to use and merge them?

[timoast]

Well that's up to you, all I'm saying is there's no bug here. You should try to work out if you think these differences are genuinely biological, or whether they're due to technical differences. You could check the QC metrics for those datasets and see if there's anything to indicate there's something wrong with them, and look at other regions of the genome to see if this difference is specific to this region or is characteristic of the genome-wide signal.

[Ping-lin14]

OK, thank you very much.