ATACseqQC - Cut-site probability greater than 1?

Question

[ldyer2006]

Hi guys,

I was just wondering if anybody here had any experience working with the ATACseqQC R package?

The package seems to function well for me over all, and gives reasonable results when I apply it to the entire mouse genome (Using the following command):

sigs.Irx3.0 <- factorFootprints(bamfiles = "Day0.shifted.bam", pfm = pwm,
                                                genome = Mmusculus, index = "Day0.shifted.bam",
                                                min.score="95%", seqlev=paste0("chr", c(1:19, "X", "Y")),
                                                upstream=100,downstream=100)

The problem then comes when I try to examine an individual chromosome, say chr1:

sigs.Irx3.0 <- factorFootprints(bamfiles = "Day0.shifted.bam", pfm = pwm,
                                                genome = Mmusculus, index = "Day0.shifted.bam",
                                                min.score="95%", seqlev=paste0("chr1"),
                                                upstream=100,downstream=100)

The function still runs, however it produces a graph with many spikes of cut-site probability greater than 1. I'm assuming that this is erroneous. If it's not, I can't find any details on the package's page or FAQs which would help me to interpret a cut-site probability higher than 1.

The package is designed for use with humans, as I understand it. I wonder if that has anything to do with the incorrect calculation of the cut-site. If so, are there any known solutions to this issue?

[GenoMax]

Please use the formatting bar (especially the code option) to format parts of you post that contain code snippets. I've done it for you this time.

[ldyer2006]

Thanks a lot, that is very helpful.

[Friederike]

can you post the graph?

[Lizzy]

How did you input you bam files using ATACseqQC? I used this commands in R: setwd("~/Documents/ATAC-Seq") where is the folder of my bam files, and the: bamfile <- "Test_ATAC-seq.bam", test <- readBamFile(bamfile, bigFile = TRUE). I checked the test and it always reported 0 sequences. While there are reads in my bamfile, it seemed that my bam file was not read in, could you please show me how to input my bam files? Thanks a lot.

Answer 1

[jianhong.ou]

Hi, Thank you for trying ATACseqQC. I am the author of ATACseqQC. If the probability greater than 1, there are two reasons: 1. most of the reads located in only one strand. 2. the reads depth are too low. If you don't mind, could you send me the files then I can improve the function to avoid this errors.

Jianhong

Answer 2

[haibol2017]

I am co-author of ATACseqQC paper (https://www.ncbi.nlm.nih.gov/pubmed/29490630). I have used 25% of aligned, duplicate-removed BAM file for one human chromosome (chr 1) to demonstrate the robustness of footprint analysis using ATACseqQC. I did see very spiky plots, but not encounter cutting probability greater than 1. Would you mind to share your chromosome-specific BAM file with us so we can try to reproduce the error and fix the potential bugs? Thanks!!