Scanbam Error

Question

[e.karasmani]

Dear All,

I have a problem and I would like your help....

I got some published ChIP-seq data....So I got the Bam and the Bai files.....and I do the following....

            chr   <- 'chr10' 

            fileName<-c("/data/lena/GG/file_1.bam")

            fragment.size=200
            chr.size<-seqlengths(Mmusculus)[chr]                 
            which <- GRanges(seqnames=chr,IRanges(1,chr.size))
            what <- c("rname", "strand", "pos", "qwidth","seq","qual")
            param <- ScanBamParam(which = which,what=what)

            chip.bam <- scanBam(fileName, param=param)

and I get the following error message

 record: 0
 error: 0
  file: /data/lena/GG/file_1.bam
  index: /data/lena/GG/file_1.bam
In addition: Warning message:
In .Call2(func, .extptr(file), space, flag, simpleCigar, ..., PACKAGE = "Rsamtools") :
 space 'chr10' not in BAM header

I can understand that has something to do with the header BUT how can I fix that?

Could you please help me?

Thank you in advance

Best regards Lena

[Niek De Klein]

I don't know anything about scanBam, but it is only a warning message and it says you have 0 errors, so it should work. Does it not give you no output or not the output taht you expect?

[e.karasmani]

there is no output at all! I have a for loop initially which scans the BAM file for all the chromosomes....so the loop stops and I get the message....My script is correct because it is working perfectly with other BAM files....

[Niek De Klein]

Can you post the header of the non-working BAM file and a working BAM file?

[Jeremy Leipzig]

yeah or run scanBamHeader to see if that chromosome really exists (it might just be '10')

[e.karasmani]

thank you very much! it works! however I have now another problem....if you can see the post below....can you please help me? thank you in advance. best regards. Lena

Answer 1

[Steve Lianoglou]

It looks like the error is telling you that there is no "chromosome" named chr10 in your BAM file. To verify, try this:

library(Rsamtools)
bf <- BamFile("/data/lena/GG/file_1.bam")
seqnames(seqinfo(bf))

Do you see chr10 in the name of the chromosomes there, or?

Update to address comment

Sot it looks like the index that your reads were aligned to used chromosome names without the chr prefix -- you must to refer to the chromosomes in the same way, so instead of asking for reads on "chr10", you need to get them from "10".

Try changing your for loop, from this:

for (chr in paste('chr',c(seq(1,19),'X','Y'),sep='')) {
...
}

to this:

for (chr in c(1:19, "X", "Y")) {
...
}

That should do the trick.

[e.karasmani]

Thank you very much!

It works!

Although now when I try the following

              chip.qual    <-QualityScaledBStringSet(chip.bam[[name]]$seq, chip.bam[[name]]$qual)

i get that error

Error in XStringSet("B", x, start = start, end = end, width = width, use.names = use.names) :
error in evaluating the argument 'x' in selecting a method for function 'XStringSet': Error in chip.bam[[name]] : no such index at level 2

do you know what should I do?

Thank you very much! I would appreciate your help!

best regards Lena

[Steve Lianoglou]

First: please don't update your questions via the comments -- there was a previous comment here you posted about that I responded to in my answer, but now the comment has changed to this one about quality scores.

Since this is now a different problem, please post a new question that gives us enough information to diagnose the problem -- for example: what is stored in name? What types of objects are chip.bam[[name]]$seq and chip.bam[[name]]$qual and are they what you expect them to be?

[e.karasmani]

Thank you very much for your answer and my sincere apologies for the problem, I changed by accident....I wanted just to put my new question....

So I still have a problem.....

So in total here is my script.....

        chr   <- '19' 

        fileName<-c("/data/lena/GG/file_1.bam")

        fragment.size=200
        chr.size<-seqlengths(Mmusculus)[]                 
        which <- GRanges(seqnames=chr,IRanges(1,chr.size))
        what <- c("rname", "strand", "pos", "qwidth","seq","qual")
        param <- ScanBamParam(which = which,what=what)

        chip.bam <- scanBam(fileName, param=param)
        name <- names(chip.bam)
        name


 [1] "19:1-197195432" "19:1-181748087" "19:1-159599783" "19:1-155630120"
 [5] "19:1-152537259" "19:1-149517037" "19:1-152524553" "19:1-131738871"
 [9] "19:1-124076172" "19:1-129993255" "19:1-121843856" "19:1-121257530"
[13] "19:1-120284312" "19:1-125194864" "19:1-103494974" "19:1-98319150"
[17] "19:1-95272651"  "19:1-90772031"  "19:1-61342430"  "19:1-166650296"
[21] "19:1-15902555"  "19:1-16299"     "19:1-1231697"   "19:1-41899"
[25] "19:1-160594"    "19:1-357350"    "19:1-362490"    "19:1-849593"
[29] "19:1-449403"    "19:1-400311"    "19:1-3994"      "19:1-628739"
[33] "19:1-1785075"   "19:1-58682461"  "19:1-5900358"


 chip.qual    <-QualityScaledBStringSet(chip.bam[[name]]$seq, chip.bam[[name]]$qual)

 Error in XStringSet("B", x, start = start, end = end, width = width, use.names = use.names) :
 error in evaluating the argument 'x' in selecting a method for function 'XStringSet': Error in chip.bam[[name]] : no such index at level 2

So what am I doing wrong? could you please help me?

[Steve Lianoglou]

As I said, please submit a new question to biostar since this is now different -- this site works better as a Q/A site, not long winded "support" threads. You should think about if your name vector is sane ... you're doing something crazy there w/ your which variable you've built from chr.size

[e.karasmani]

done! Thank you!

[Steve Lianoglou]

Can you provide a link to your question? I can't seem to find it.