splitGAlignmentsByCut returns empty .bam files

[chklee30]

The function runs and creates the NucleosomeFree.bam, mononucleosome.bam, dinucleosome.bam, and trinucleosome.bam files, but these files only contain some header info. There are no actual alignment records.

The only message I get during execution of the function is the following:

[W::sam_read1_sam] Parse error at line 196
[W::sam_read1_sam] Parse error at line 196
[W::sam_read1_sam] Parse error at line 196
[W::sam_read1_sam] Parse error at line 196
[W::sam_read1_sam] Parse error at line 196
[W::sam_read1_sam] Parse error at line 196
[W::sam_read1_sam] Parse error at line 196
[W::sam_read1_sam] Parse error at line 196

There are no other warnings or errors. I'm not sure how to diagnose this problem.

P.S - what is line 196 here referring to? The object I am passing into splitGAlignmentsByCut is a GAlignments vector-like object of length N. Is it somehow accessing the original bam file that was read in to create the GAlignments object? In either case, I did not see anything weird with line 196 in the bam file.

Thanks for any help.

[jianhong]

which version are you using?Best!Your sincerely,Jianhong OuOn Oct 25, 2024, at 5:28 PM, chklee30 ***@***.***> wrote: The function runs and creates the NucleosomeFree.bam, mononucleosome.bam, dinucleosome.bam, and trinucleosome.bam files, but these files only contain some header info. There are no actual alignment records. The only message I get during execution of the function is the following: [W::sam_read1_sam] Parse error at line 196 [W::sam_read1_sam] Parse error at line 196 [W::sam_read1_sam] Parse error at line 196 [W::sam_read1_sam] Parse error at line 196 [W::sam_read1_sam] Parse error at line 196 [W::sam_read1_sam] Parse error at line 196 [W::sam_read1_sam] Parse error at line 196 [W::sam_read1_sam] Parse error at line 196 There are no other warnings or errors. I'm not sure how to diagnose this problem. P.S - what is line 196 here referring to? The object I am passing into splitGAlignmentsByCut is a GAlignments vector-like object of length N. Is it somehow accessing the original bam file that was read in to create the GAlignments object? In either case, I did not see anything weird with line 196 in the bam file. Thanks for any help. —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[chklee30]

Hello @jianhong ,

I am using ATACseqQC version 1.26.0

[jianhong]

Do you have trouble to install the newest version? I think this error is due to the NA values in some of the tags. There are several solutions, 1, updates the package. 2, remove some of the reads tags that not necessary; 3, load the patches for the export function.Best!Your sincerely,Jianhong OuOn Oct 25, 2024, at 7:53 PM, chklee30 ***@***.***> wrote: Hello @jianhong , I am using ATACseqQC version 1.26.0 —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>

[chklee30]

Thank you for the reply @jianhong

I have updated the package to version 1.29.0 but unfortunately get the same error.

For your second suggestion, how can I code in R to remove the NA tags in the GAlignments object?

Thank you

[jianhong]

Thank you for letting me know that. a new bug. Could you please share me the files to repeat the error?Best!Your sincerely,Jianhong OuOn Oct 26, 2024, at 9:43 AM, chklee30 ***@***.***> wrote: Thank you for the reply @jianhong I have updated the package to version 1.29.0 but unfortunately get the same error. For your second suggestion, how can I code in R to remove the NA tags in the GAlignments object? Thank you —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>

[chklee30]

Hi @jianhong ,

How should I send you the file? The RDS file for the GAlignments object is 220 MB. I cannot attach it here on github.

[jianhong]

just subset the object to one small region of chromosome to repeat the error will be OK.Best!Your sincerely,Jianhong OuOn Oct 26, 2024, at 10:00 AM, chklee30 ***@***.***> wrote: Hi @jianhong , How should I send you the file? The RDS file for the GAlignments object is 220 MB. I cannot attach it here on github. —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>

[chklee30]

sample_galignments_obj.rds.gz

Hi @jianhong ,

Here is the file.

[jianhong]

Thank for the file. I pushed a update to the development version. Please try it by BiocManager::install('jianhong/ATACseqQC') the version should be 1.29.1

[chklee30]

Thank you @jianhong the splitGAlignmentsByCut() function works now. However, when I move on to the downstream steps as follows, I get another error running the enrichedFragments() function:

objs <- splitGAlignmentsByCut(gal1, txs=txs, outPath = outdir)

bamfiles <- file.path(outdir,
                     c("NucleosomeFree.bam",
                     "mononucleosome.bam",
                     "dinucleosome.bam",
                     "trinucleosome.bam"))
TSS <- promoters(txs, upstream=0, downstream=1)
TSS <- unique(TSS)
## estimate the library size for normalization
(librarySize <- estLibSize(bamfiles))

NTILE <- 101
dws <- ups <- 1010
sigs <- enrichedFragments(gal=objs[c("NucleosomeFree", 
                                     "mononucleosome",
                                     "dinucleosome",
                                     "trinucleosome")], 
                          TSS=TSS,
                          librarySize=librarySize,
                          TSS.filter=0.5,
                          n.tile = NTILE,
                          upstream = ups,
                          downstream = dws)

The error is:

Error in .normarg_f(f, x) : 
  split factor has length 0 but 'NROW(x)' is > 0

Could you help me with this? I very much appreciate your time on this. Thank you.

[jianhong]

Could you please run traceback() to show more error message?

[chklee30]

@jianhong

12: stop("split factor has length 0 but 'NROW(x)' is > 0")
11: .normarg_f(f, x)
10: IRanges:::default_splitAsList(x, f, drop = drop)
9: .local(x, f, drop, ...)
8: splitAsList(x, f, drop = drop, ...)
7: splitAsList(x, f, drop = drop, ...)
6: split(ga, qname)
5: split(ga, qname)
4: (function (ga, .fLen) 
   {
       if (is(ga, "GAlignmentPairs")) {
           granges(.ele)
       }
       else {
           qname <- mcols(ga)$qname
           if (pe == "auto") {
               if (length(qname) == 0) {
                   pe <- FALSE
               }
               else {
                   if (all(table(qname) == 1)) {
                     pe <- FALSE
                   }
                   else {
                     if (all(table(qname) < 3)) {
                       pe <- TRUE
                     }
                     else {
    ...
3: mapply(function(ga, .fLen) {
       if (is(ga, "GAlignmentPairs")) {
           granges(.ele)
       }
       else {
           qname <- mcols(ga)$qname
           if (pe == "auto") {
               if (length(qname) == 0) {
                   pe <- FALSE
               }
               else {
                   if (all(table(qname) == 1)) {
                     pe <- FALSE
                   }
                   else {
                     if (all(table(qname) < 3)) {
                       pe <- TRUE
                     }
                     else {
                       pe <- FALSE
    ...
2: featureAlignedExtendSignal(gal = gal, feature.gr = TSS, upstream = upstream, 
       downstream = downstream, n.tile = n.tile, fragmentLength = fragmentLength, 
       librarySize = librarySize, adjustFragmentLength = adjustFragmentLength, 
       pe = "PE")
1: enrichedFragments(gal = objs[c("NucleosomeFree", "mononucleosome", 
       "dinucleosome", "trinucleosome")], TSS = TSS, librarySize = librarySize, 
       TSS.filter = 0.5, n.tile = NTILE, upstream = ups, downstream = dws)

Here you go