This is a test version of Biostars. For the public version, visit https://www.biostars.org.
I'm running a scATAC-seq analysis. However, the peaks I've called are somehow weird

I'm running a scATAC-seq analysis. The data was obtained from the published work: Single-Cell Chromatin Accessibility Analysis Reveals the Epigenetic Basis and Signature Transcription Factors for the Molecular Subtypes of Colorectal Cancers; https://ngdc.cncb.ac.cn/omix/; OMIX005759; OMIX005759-04.sort.fragments.tsv.gz. However, when I run the following code, I call these peaks:

peaks2 <- CallPeaks(frags2,macs2.path="/home/anaconda3/envs/lowpython/bin/macs2",combine.peaks=F)
peaks2@ranges
#IRanges object with 24 ranges and 0 metadata columns:
#           start       end     width
#       <integer> <integer> <integer>
#   [1]      9968 248946051 248936084
#   [2]     10191 133787023 133776833
#   [3]    146452 135076001 134929550
#   [4]     15323 133236135 133220813
#   [5]  18172397 114353789  96181393
#   ...       ...       ...       ...
#  [20]     19503 159335730 159316228
#  [21]     82336 145075845 144993510
#  [22]     10409 138233182 138222774
#  [23]   2781347 156030667 153249321
#  [24]    268144  57217121  56948978

They are obviously too wide to be "peaks"! My complete code is as follows:

## PATH
fragpath2 <- "/data/Tide/rawdata/TangLab-scATAC-seq/OMIX005759-04.sort.fragments.tsv.gz"
## require essential informations
total_counts <- CountFragments(fragpath)
cutoff <- 1000 # Change this number depending on your dataset!
barcodes <- total_counts[total_counts$frequency_count > cutoff, ]$CB
## create .tbi index
gunzip rawdata/TangLab-scATAC-seq/OMIX005759-04.sort.fragments.tsv.gz
bgzip rawdata/TangLab-scATAC-seq/OMIX005759-04.sort.fragments.tsv
tabix -p bed rawdata/TangLab-scATAC-seq/OMIX005759-04.sort.fragments.tsv.gz
## Create fragment object
frags2 <- CreateFragmentObject(path = fragpath2, cells = barcodes)
## peak calling
peaks2 <- CallPeaks(frags2,macs2.path="/home/anaconda3/envs/lowpython/bin/macs2")

The version of RStudio: 4.3.1 The version of Signac: 1.14.0 The version of MACS2: 2.1.4 How could I fix this problem?

macs2 scatac-seq signac

In addition, I tried the example data from 10x Genomics, and the result doesn't seem to have any problems.

peaks@ranges
#IRanges object with 102411 ranges and 0 metadata columns:
#               start       end     width
#           <integer> <integer> <integer>
#       [1]     10070     10467       398
#       [2]    565201    565409       209
#       [3]    569285    569496       212
#       [4]    713482    714649      1168
#       [5]    752323    752810       488
#       ...       ...       ...       ...
#  [102407]  59004755  59005111       357
#  [102408]  59015825  59016323       499
#  [102409]  59016964  59017177       214
#  [102410]  59029917  59030163       247
#  [102411]  59031616  59032089       474

0 answers

No answers yet.

Log in to answer this question.