cn.MOPs: Error in value[[3L]](cond) :
3
0
Entering edit mode
@laurabuggiotti-11728
Last seen 5.6 years ago

Hi,

I have been using cn.mops successfully using the following commands:

bamDataRanges <- getReadCountsFromBAM(BAMFiles, refSeqName=readLines("Chr_list.txt"), WL=1000)

results <- calcIntegerCopyNumbers(cn.mops(bamDataRanges))

segm <- as.data.frame(segmentation(results))
CNVs <- as.data.frame(cnvs(results))
CNVRegions <- as.data.frame(cnvr(results))

However, im now running it on additional bam files and i get the following error:

Error in value[[3L]](cond) : 
  failed to open BamFile: failed to load BAM index
  file: A27.realigned.rmdup.bam
Calls: getReadCountsFromBAM ... tryCatch -> tryCatchList -> tryCatchOne -> <Anonymous>
In addition: Warning message:
In doTryCatch(return(expr), name, parentenv, handler) :
  [bam_index_load] fail to load BAM index.
Execution halted

 

i have checked the header of my bam file and looks fine, would you please help me in trying to fix the error?

Thanks for your time and support,

Laura

 

 

cn.mops • 4.3k views
ADD COMMENT
0
Entering edit mode
@gunter-klambauer-5426
Last seen 3.2 years ago
Austria

I think this is not a problem of cn.mops but of an internally used function (Rsamtools?).

  • Could you run "traceback" right after the error occurs and post the result here.
  • Please check if for each BAM file, there exists an BAI (index) file. This can be generated using samtools

Regards,
Günter

ADD COMMENT
0
Entering edit mode

Hi again,

i did try various bam files and it seems that the error occurs when using recalibrated bam files (after the BQSR step)...do you reckon there might be a link with the issue? Did you ever experience the same issue?

Thanks a lot for your time,

Laura

ADD REPLY
0
Entering edit mode
@laurabuggiotti-11728
Last seen 5.6 years ago

Thanks for your prompt answer, here is the traceback:

> traceback()
11: stop("failed to open BamFile: ", conditionMessage(err))
10: value[[3L]](cond)
9: tryCatchOne(expr, names, parentenv, handlers[[1L]])
8: tryCatchList(expr, classes, parentenv, handlers)
7: tryCatch({
       .io_check_exists(path(con))
       index <- sub("\\.bai$", "", index(con, asNA = FALSE))
       con$.extptr <- .Call(.bamfile_open, path(con), index, "rb")
   }, error = function(err) {
       stop("failed to open BamFile: ", conditionMessage(err))
   })
6: open.BamFile(BamFile(file, index), "rb")
5: open(BamFile(file, index), "rb")
4: scanBam(bam.file, param = ScanBamParam(what = scan.what, which = range(granges.subset)))
3: scanBam(bam.file, param = ScanBamParam(what = scan.what, which = range(granges.subset)))
2: exomeCopy::countBamInGRanges(bam.file = BAMFiles[i], granges = GR,
       ...)
1: getReadCountsFromBAM(BAMFiles, refSeqName = readLines("/users/lbuggiotti/WGresRussia/cn.MOPS/Chr_list.txt"),
       WL = 1000)

Moreover i have BAI file for each bam created by using samtools.

Thanks again,

Laura

ADD COMMENT
0
Entering edit mode
@gunter-klambauer-5426
Last seen 3.2 years ago
Austria

Ok, this is a problem either with your BAM files or with Rsamtools with the function "scanBam". Please contact the maintainers of Rsamtools.

ADD COMMENT
0
Entering edit mode

Is the bam file index '.bai' in the same directory as the bam file, and following a standard naming convention? You could simplify debugging with

library(Rsamtools)
open(BamFile("path/to/file.bam"))

(updated to the correct path to one of your bam files as presented in the input to getReadCountsFromBAM().

ADD REPLY
0
Entering edit mode

yes, it is in the same directory, here one example of the name:

A27.realigned.rmdup.bai
A27.realigned.rmdup.bam

and got no error when opening the bam file as you suggested (via Rsamtools).

anything else you could think of?

 

ADD REPLY
0
Entering edit mode

Look carefully at BAMFiles, and make sure that you can open each element via open(BamFile(BAMFiles[i])). Programmatically, you could do

trace(Rsamtools:::open.BamFile, quote(print(path(con))))
bamDataRanges <- getReadCountsFromBAM(BAMFiles, refSeqName=readLines("Chr_list.txt"), WL=1000)

Each time Rsamtools tries to open a BAM file, it will print the path to the file. The last path printed is the one that causes problems.

ADD REPLY
0
Entering edit mode

i have tried with one bam file as you suggested:

Tracing open.BamFile(BamFile(file, character(0))) on entry
[1] "/users/lbuggiotti/bamKorea/A27.realigned.rmdup.bam"

Identified the following reference sequences:  Chr1,Chr2......
.....
Missing "refSeqNames"! Selecting all identified reference sequence for analysis.

PLEASE BE PATIENT... this might take a while. Consider using the parallel version of this function

Error in if (WL >= sl[i]) { : missing value where TRUE/FALSE needed

Now it seems to be the window length the issue...

Thanks a lot for your help!

ADD REPLY
0
Entering edit mode

I think I did not communicate correctly. Try just opening each element of BAMFile with

for (i in seq_along(BAMfiles)
    open(BamFile(BAMFiles[i]))

If that doesn't help then try

trace(Rsamtools:::open.BamFile, quote(print(path(con))))
bamDataRanges <- getReadCountsFromBAM(BAMFiles, refSeqName=readLines("Chr_list.txt"), WL=1000)

Your original error was about opening BAM files, so let's solve that first before moving on to errors in the cn.mops package (or errors that only occur while trying to debug other errors).

ADD REPLY
0
Entering edit mode

Ok..this is what i did:

> bamFile <- open(BamFile("bamKorea/A27.realigned.rmdup.bam"))
> bamFile
class: BamFile
path: bamKorea/A27.realigned.rmdup.bam
index: bamKorea/A27.realigned.rmdup.bai
isOpen: TRUE
yieldSize: NA
obeyQname: FALSE
asMates: FALSE

> BAMFiles = ("bamKorea/SAMN02225733.realigned.rmdup.bam")
> trace(Rsamtools:::open.BamFile, quote(print(path(con))))
Tracing function "open.BamFile" in package "Rsamtools (not-exported)"
[1] "open.BamFile"
> bamDataRanges <- getReadCountsFromBAM(BAMFiles, refSeqName=readLines("/users/lTracing open.BamFile(BamFile(file, character(0))) on entry
[1] "bamKorea/SAMN02225733.realigned.rmdup.bam"
Identified the following reference sequences:  Chr1,Chr2,Chr3,Chr4,Chr5,Chr6...

 

Using Chr1, Chr2, Chr3, Chr4, Chr5, Chr6, .... as reference.

PLEASE BE PATIENT... this might take a while. Consider using the parallel version of this function

Tracing open.BamFile(BamFile(file, character(0))) on entry
[1] "bamKorea/SAMN02225733.realigned.rmdup.bam"
Tracing open.BamFile(BamFile(file, index), "rb") on entry
[1] "bamKorea/SAMN02225733.realigned.rmdup.bam"
Error in value[[3L]](cond) :
  failed to open BamFile: failed to load BAM index
  file: bamKorea/SAMN02225733.realigned.rmdup.bam
In addition: Warning message:
In doTryCatch(return(expr), name, parentenv, handler) :
  [bam_index_load] fail to load BAM index.

Im back to square one...the problem is the index? I used samtools to produce the .bai for each of the bam; plus when i open bamFile it does see that there is a correspondent .bai.

thanks a lot for your time!

 

Hi again, i did try various bam files and it seems that the error occurs when using recalibrated bam files (after the BQSR step)...do you reckon there might be a link with the issue?

Thanks,

Laura

ADD REPLY
0
Entering edit mode

So what happens when you try to open the problematic BAM file

open(BamFile("bamKorea/SAMN02225733.realigned.rmdup.bam""))

?

ADD REPLY
0
Entering edit mode

it doesnt complain:

​> open(BamFile("bamKorea/SAMN02225733.realigned.rmdup.bam"))

>

Im running the scan Bam at the moment and will need sometime...

Rsamtools::scanBam("bamKorea/SAMN02225733.realigned.rmdup.bam")

 

ADD REPLY
0
Entering edit mode

Have you tried using the full path to the BAM files in the getReadCountsFromBAM() command? By full path I mean from the root of the directory, along the lines of /home/mtmorgan/projects/bamKorea/SAMN02225733.realigned.rmdup.bam . You might also modify the trace command to

trace(Rsamtools:::open.BamFile, quote({ print(path(con)); print(index(con)) }))

so that it prints both the bam file and index path. I do not think the problem is with scanBam.

ADD REPLY
0
Entering edit mode
> BAMFile="bamKorea/SAMN02225733.realigned.rmdup.bam"

>bamDataRanges <- getReadCountsFromBAM(BAMFile, refSeqName=readLines("/users/lbuggiotti/WGresRussia/cn.MOPS/Chr_list.txt"), WL=1000)

I get the same error:

Error in value[[3L]](cond) :

  failed to open BamFile: failed to load BAM index

  file: bamKorea/SAMN02225733.realigned.rmdup.bam

In addition: Warning message:

In doTryCatch(return(expr), name, parentenv, handler) :

  [bam_index_load] fail to load BAM index.
ADD REPLY
0
Entering edit mode

I meant the full path to the BAM file. Also please use the updated trace command.

ADD REPLY
0
Entering edit mode
Tracing open.BamFile(BamFile(file, character(0))) on entry

[1] "/users/lbuggiotti/bamKorea/SAMN02225733.realigned.rmdup.bam"

[1] NA

Tracing open.BamFile(BamFile(file, index), "rb") on entry

[1] "/users/lbuggiotti/bamKorea/SAMN02225733.realigned.rmdup.bam"

[1] "/users/lbuggiotti/bamKorea/SAMN02225733.realigned.rmdup.bam"

Error in value[[3L]](cond) :

  failed to open BamFile: failed to load BAM index

  file: /users/lbuggiotti/bamKorea/SAMN02225733.realigned.rmdup.bam

In addition: Warning message:

In doTryCatch(return(expr), name, parentenv, handler) :

  [bam_index_load] fail to load BAM index.

if i open BamFile:

> open(BamFile("/users/lbuggiotti/bamKorea/SAMN02225733.realigned.rmdup.bam"))

Tracing open.BamFile(BamFile("/users/lbuggiotti/bamKorea/SAMN02225733.realigned.rmdup.bam")) on entry

[1] "/users/lbuggiotti/bamKorea/SAMN02225733.realigned.rmdup.bam"

[1] "/users/lbuggiotti/bamKorea/SAMN02225733.realigned.rmdup.bai"
ADD REPLY
0
Entering edit mode

Can you provide the output of the commands

packageVersion("Rsamtools")
sessionInfo()

A workaround might rename the index files from

bamKorea/SAMN02225733.realigned.rmdup.bai"

to

bamKorea/SAMN02225733.realigned.rmdup.bam.bai
ADD REPLY
0
Entering edit mode

when renaming the index file i dont get the error...thank you!
Ill modify all indexes and i hope all the samples will be analysed with cn.MOPS.

Thank you so much for your help!
cheers,

Laura

ADD REPLY

Login before adding your answer.

Traffic: 950 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6