I'm running a scATAC-seq analysis. However, the peaks I've called are somehow weird
0
0
Entering edit mode
9 weeks ago
Pallondyle • 0

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 • 221 views
ADD COMMENT
0
Entering edit mode

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
ADD REPLY

Login before adding your answer.

Traffic: 1924 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6