search

stuart-lab / signac

R toolkit for the analysis of single-cell chromatin data
https://stuartlab.org/signac/
Other
327 stars 88 forks source link

CreatChromatinAssay changes number of peaks! #207

Closed fkeramati closed 4 years ago

fkeramati commented 4 years ago

Hi,

Last week, unintentionally I upgraded my Signac version from 0.2 to 1.0. Now the older version script is not working. I started changing the code from last version to the newer one. I have 3 human bone-marrow samples from the same donor over the time course. In the previous version, I could merge them very well without the need for any batch correction or anchoring. But, now in the new version, when I do the merging without any correction, each time point sample makes a separate population of cells!!! What I do is that I have a unified peak set that I used "FeatureMatrix" function to count fragments in each peak for each sample and then merge them. This peak set has 179266 peaks. Surprisingly, when I generate the seurat object from this set each sample gets a new number of peaks. Sample 1 now has 179210, sample 2 has 179255 and sample 3 has 179223 peaks!!! What is wrong here and why the code changes the number of peaks without notifying or giving any information! 

# This is an example for one of the objects (I use cellRanger peak set here (70214 peaks)
# When I use ChromatinAssay object the number of peaks decrease to 69033
# While using the original counts matrix (like in the previous version of Signac) it remains the same number.
counts = Read10X_h5(filename = "raw_peak_bc_matrix.h5")
counts = counts[,colnames(counts) %in% verified_cells] # I have a table of verified cells from ArchR and want to use those.
metadata = read.csv(file = "singlecell.csv", header = T, row.names = 1)
metadata = metadata[rownames(metadata) %in% verified_cells]
chrom_obj = CreateChromatinAssay(counts = counts, sep = c(":", "-"), genome = "hg38", fragments = "fragments.tsv.gz")
bm_1 = CreateSeuratObject(counts = chrom_obj, assay = 'ATAC', meta.data = metadata) # 69033 features 5173 samples
bm_2 = CreateSeuratObject(counts = counts, assay = "ATAC", meta.data = metadata) # 70214 features 5173 samples
fkeramati commented 4 years ago

I could find the answer actually. It is coming from "min.cells = 0". When I change it to "min.cells = -1", it solves the problem.