BAM Header for ATAC-seq

Question

[sophia.dingg]

Hello! I am currently trying to follow ATACseqQC in adjusting the read start sites using shiftGAlignmentsList (https://bioconductor.org/packages/release/bioc/vignettes/ATACseqQC/inst/doc/ATACseqQC.html) on my own files; I generated BAM files from fast.qc files in linux.

I have the lines of code in R (following what the author did):

library(BSgenome.Hsapiens.NCBI.GRCh38)
bamfile <- file.path("P.bam")
possibleTag <- combn(LETTERS, 2)
possibleTag <- c(paste0(possibleTag[1, ], possibleTag[2, ]),
              paste0(possibleTag[2, ], possibleTag[1, ]))
bamTop100 <- scanBam(BamFile(bam_SS123, yieldSize = 100),
                     param = ScanBamParam(tag = possibleTag))[[1]]$tag
tags <- names(bamTop100)[lengths(bamTop100)==100]

outPath <- "splited"
dir.create(outPath)
which <- as(seqinfo(Hsapiens)['1'], "GRanges")
gal <- readBamFile(bam_SS123, tag=tags, which=which, asMates=TRUE, bigFile=TRUE)
shiftedBam_SS123 <- file.path(outPath, "shifted_SS123.bam")
gal1 <- shiftGAlignmentsList(gal, outbam = shiftedBam_SS123)

When I do so, I am getting the error that "seq levels(param) not in BAM header". I checked the header of my BAM file, and they starkly differ from the seqlevels that I got from the library of NCBI GRCh38, even though that's what I used to create my indexes for generating my bam files.

I would greatly appreciate if anyone had any advice on how to tackle this!