Issue with na.rm=FALSE

[arijita88]

Hi @ssvQC, I am trying to run ssvQC on my Cut&Run files and I am getting this error:

Error in quantile.default(frip_dt$frip, c(0.1, 0.96)) : missing values and NaN's not allowed if 'na.rm' is FALSE

Is there a way to set na.rm=TRUE while running the package?

Thank you

[jrboyd]

Hi, thanks for the question and sorry about the delay in getting back to you. Crazy times.

That error looks like it's coming from the plots related to FRIP. I think there's probably something going wrong upstream that ssvQC should maybe be catching.

I can only guess, but somehow you probably either have 0 total mapped reads in a bam file or 0 peaks for one or more groups.

Can you do me a favor and instead of either running ssvQC.runAll() or ssvQC.plotFRIP(), can you run ssvQC.prepFRIP()?

So:

#where sqc is a valid ssvQC object:
sqc = ssvQC.prepFRIP(sqc)
#then can you look at a few attributes:
lengths(sqc$FRIP)
range(sqc$FRIP[[1]][[1]]$mapped_reads)
range(sqc$FRIP[[1]][[1]]$frip)

If nothing sticks out to you, maybe you could send me the whole sqc$FRIP table.

The solution might be na.rm = TRUE like you've proposed but I'd like to catch the underlying error if there is one. Thanks!

[arijita88]

Thanks so much for your reply. This is really helpful. There aren't any mapped reads.

range(sqc$FRIP[[1]][[1]]$mapped_reads) [1] 0 0 range(sqc$FRIP[[1]][[1]]$frip) [1] NaN NaN

I didn't mention the bam read_mode to be PE which is leading to this problem may be?

Is there any way to assign the read_mode argument without using the config file? If not, can you please explain the details to be provided in the config file. I went through the example config file and I am not sure if I understand all the options provided there.

Thanks

[jrboyd]

I don't think read_mode has an impact on how mapped reads gets calculated.

How many bam files are you looking at? Do they all have 0 mapped reads?

Could you try running Rsamtools::idxstatsBam on all of them? You should get a table that looks like this:

   seqnames seqlength mapped unmapped
1      chr1 248956422 162980        0
2      chr2 242193529 156276        0
3      chr3 198295559 135520        0

The sum of the "mapped" column is what gets used as mapped_reads per bam file.

[arijita88]

I am running two BAM files. They do have a lot of mapped reads.

When I run Rsamtools::idxstatsBam one each of the files I get the following result:

Rsamtools::idxstatsBam("file1.bam")

| seqnames | seqlength | mapped | unmapped

1 | 1 | 249250621 | 353025 | 0 2 | 2 | 243199373 | 367821 | 0 3 | 3 | 198022430 | 253461 | 0 4 | 4 | 191154276 | 276273 | 0 5 | 5 | 180915260 | 231007 | 0 6 | 6 | 171115067 | 249537 | 0 7 | 7 | 159138663 | 242854 | 0 8 | 8 | 146364022 | 211374 | 0 9 | 9 | 141213431 | 182834 | 0 10 | 10 | 135534747 | 225268 | 0 11 | 11 | 135006516 | 189472 | 0 12 | 12 | 133851895 | 190579 | 0 13 | 13 | 115169878 | 141542 | 0 14 | 14 | 107349540 | 116443 | 0 15 | 15 | 102531392 | 115523 | 0 16 | 16 | 90354753 | 133882 | 0 17 | 17 | 81195210 | 111606 | 0 18 | 18 | 78077248 | 118658 | 0 19 | 19 | 59128983 | 88431 | 0 20 | 20 | 63025520 | 83148 | 0 21 | 21 | 48129895 | 58120 | 0 22 | 22 | 51304566 | 50862 | 0 23 | X | 155270560 | 139970 | 0 24 | Y | 59373566 | 50566 | 0 25 | MT | 16569 | 5807 | 0

Rsamtools::idxstatsBam("file2.bam")

| seqnames | seqlength | mapped | unmapped

1 | 1 | 249250621 | 583938 | 0 2 | 2 | 243199373 | 617374 | 0 3 | 3 | 198022430 | 462514 | 0 4 | 4 | 191154276 | 482485 | 0 5 | 5 | 180915260 | 455838 | 0 6 | 6 | 171115067 | 397892 | 0 7 | 7 | 159138663 | 387254 | 0 8 | 8 | 146364022 | 372774 | 0 9 | 9 | 141213431 | 287158 | 0 10 | 10 | 135534747 | 362260 | 0 11 | 11 | 135006516 | 316757 | 0 12 | 12 | 133851895 | 319331 | 0 13 | 13 | 115169878 | 238195 | 0 14 | 14 | 107349540 | 209815 | 0 15 | 15 | 102531392 | 187925 | 0 16 | 16 | 90354753 | 199700 | 0 17 | 17 | 81195210 | 167935 | 0 18 | 18 | 78077248 | 188446 | 0 19 | 19 | 59128983 | 123162 | 0 20 | 20 | 63025520 | 129372 | 0 21 | 21 | 48129895 | 96856 | 0 22 | 22 | 51304566 | 75070 | 0 23 | X | 155270560 | 317548 | 0 24 | Y | 59373566 | 30757 | 0 25 | MT | 16569 | 11583 | 0

The only thing is I don't have chromosomes as 'chr1,chr2...' instead I have '1,2...'. Not sure if this matters.

[jrboyd]

I found the bug! It was related to the pattern of chromosome names. Thanks for helping me track this down.

If you reinstall from github your issue should be fixed. You should be on version 1.02 of ssvQC after updating.

[arijita88]

I reinstalled from GitHub and my ssvQC version is 1.0.2. But I am still getting the same error:

Error in quantile.default(frip_dt$frip, c(0.1, 0.96)) : missing values and NaN's not allowed if 'na.rm' is FALSE

[jrboyd]

OK, sorry that didn't work right away. I think the issue is you have cached FRIP results that are getting loaded with the 0 mapped_reads counts.

Update again to ssvQC version 1.0.3 then before running force an overwrite of the cache like so: SQC_OPTIONS$SQC_FORCE_CACHE_OVERWRITE = TRUE

[arijita88]

Works now! Thank you for solving the issue so quickly. Appreciate it.