ATACSeqQC.r: fails for single end sequences

[PollyTikhonova]

I've got a error in shiftGAlignmentsList(gal) : is(gal, "GAlignmentsList") is not TRUE

As far as I understand, that's because readBamFile has a parameter asMates=FALSE and, at this case, readBamFile generates A GAlignments object but not A GAlignmentsList object.

So, this part of the code is invalid

# works only when the input is single end
if (opt$PE == 0) {
	# shift the BAM file - forward strand by +4 bp
	possibleTag <- c("AS", "XN", "XM", "XO", "XG", "NM", "MD", "YS", "YT")
	gal <- readBamFile(bamfile, asMates=FALSE, bigFile=TRUE)
	shiftedBamfile <- paste0(outdir, '/shifted.bam')
	gal1 <- shiftReads(gal)
	export(gal1, shiftedBamfile)
	cat(sprintf("\n *** shifted bam file ****  \n"))
}```

[PollyTikhonova]

And, by the way, I think it is better to set bigFile=TRUE parameter for readBamFile function. (But I could be wrong, cause I did not use ATACSeqQC that much)

[ay-lab]

HI @PollyTikhonova Thanks for using our pipeline. If you are using the latest code (updated 5 Months ago), it avoids using ATACseqQC since we found some problems in importing this package. We have used deepTools to extract nucleosome free and nucleosome containing regions.

Are you using the script "Sample_ATACseqQC_script.r" by any chance? In such a case, we recommend not to use them. Rather, please run the "pipeline.sh" script.

[PollyTikhonova]

I run this script on purpose, not as part of any of your pipelines. (just wanted to see how this tool works). But I understood you and will use your main script in future)