There are some issues when running ATACseqQC.

1. In out-pdf functions such as fragSizeDist, PTscore, only when I exit R session can I get non-empty pdf document, and if I don't exit R, the out-pdf is always empty, so that run a function, I have to go into R, and then exit to get the target file. Is this the way it was done? Or there other way to get out files?
2. When running splitGAlignmentsByCut function, I got an error.

library(ATACseqQC)
library(GenomicRanges)
bamfile <- '/simm/huangyulab/liufei/Owl/seq_data/ATACseq/GYP1.sort.clean.bam'
tags <- c("AS", "XN", "XM", "XO", "XG", "NM", "MD", "YS", "YT")
gal1 <- readBamFile(bamFile=bamfile, tag=tags,which=GRanges("chr1", IRanges(1, 1e6)),asMates=FALSE)
This is a big BAM file. Do you want to set bigFile=TRUE to save memory? (Y/n)? y
names(gal1) <- mcols(gal1)$qname
library(BSgenome.Hsapiens.UCSC.hg19)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txs <- transcripts(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(phastCons100way.UCSC.hg19)
objs <- splitGAlignmentsByCut(gal1,txs=txs,genome=genome,outPath='/simm/huangyulab/liufei/Owl/seq_data/ATACseq/GYP1.splitbam/',conservation=phastCons100way.UCSC.hg19)
  

Error in splitGAlignmentsByCut(gal1, txs = txs, genome = genome, outPath = "/simm/huangyulab/liufei/Owl/seq_data/ATACseq/GYP1.splitbam/",  :
  length(names(obj)) == length(obj) is not TRUE

and then try splitBam

> objs <- splitBam(bamfile, tags=tags, outPath='/simm/huangyulab/liufei/Owl/seq_data/ATACseq/RESULT/GYP1.splitBam/',txs=txs,genome=genome,conservation=phastCons100way.UCSC.hg19)

  Error in splitBam(bamfile, tags = tags, outPath = "/simm/huangyulab/liufei/Owl/seq_data/ATACseq/RESULT/GYP1.splitBam/", :   is(genome, "BSgenome") is not TRUE

then I set :

genome <- BSgenome.Hsapiens.UCSC.hg19

try again :

> objs <- splitBam(bamfile, tags=tags, outPath='/simm/huangyulab/liufei/Owl/seq_data/ATACseq/RESULT/GYP1.splitBam/',txs=txs,genome=genome,conservation=phastCons100way.UCSC.hg19)

error :

[E::sam_parse1] unrecognized type N 
[W::sam_read1] Parse error at line2780 
[E::sam_parse1] unrecognized type N 
[W::sam_read1] Parse error at line 2780 
[E::sam_parse1] unrecognized type N 
[W::sam_read1] Parse error at line 2780 
[E::sam_parse1] unrecognized type N 
[W::sam_read1] Parse error at line 2780 
[E::sam_parse1] unrecognized type N
[W::sam_read1] Parse error at line 2780 
[E::sam_parse1] unrecognized type N 
[W::sam_read1] Parse error at line 2780 
[E::sam_parse1] unrecognized type N 
[W::sam_read1] Parse error at line 2780
[E::sam_parse1] unrecognized type N 
[W::sam_read1] Parse error at line 2780 
Error in shiftGAlignmentsList(chunk0, positive = positive,
    negative = negative) :   all(elementNROWS(gal) < 3) is not TRUE

What does that mean and how to deal with the error when running splitGAlignmentsByCut and splitBam?

1. In out-pdf functions such as fragSizeDist, PTscore, only when I exit R session can I get non-empty pdf document, and if I don't exit R, the out-pdf is always empty, so that run a function, I have to go into R, and then exit to get the target file. Is this the way it was done? Or there other way to get out files?

A: A better way to generate plots and save to a PDF is like this: pdf(file = "XXXXX.pdf", width = 10, height = 12)

plot script here

dev.off() # this will close the graphical device for plotting, the plots in the buffer will be written to the PDF file.

1. This problem is not clear to me. I need know more about your data. 1) Are you analyze paired end ATAC-seq data from human? 2) what is the reference genome sequence you used for mapping? and where was it downloaded from? To use library(BSgenome.Hsapiens.UCSC.hg19) library(TxDb.Hsapiens.UCSC.hg19.knownGene) You must align your sequences to hg19 from UCSC. Otherwise you need make some change to the code. And son'g use library(phastCons100way.UCSC.hg19), which will trigger a computation-intensive machine-learning algorithm. 3) The third error is due the abortion of calling the function splitGAlignmentsByCut(). You should make sure the previous steps run correctly before running the next step. My suggestion is to run the example code in the ATACseqQC vignette first to see if you can get it work. Then you can try your own data. Also read our BMC Genomics Paper about the ATACseqQC package, especially the Supplementary File 1, where more detailed scripts for preprocessing of ATAC-seq data is documented.

Haibo

Hi Fei,

As Haibo suggested, to output figures, you'd better plot the figure to a device and close it after plot.

length(names(obj)) == length(obj) is not TRUE

means the length of name of gal1 is not identical to the length of gal. Please double check the names of gal1 to make sure that there is no issue with the name of the reads.

And

all(elementNROWS(gal) < 3) is not TRUE

I think this should be fixed, see https://support.bioconductor.org/p/120707/. Could you please share your sessionInfo and packageVersion of ATACseqQC.

Jianhong.

After testing, the splitGAlignmentByCut function may depend on the shiftGAlignmentsList function and use the same parameters, so I obtained adapt output successfully.

But splitBam function still not work.

Thank you all.