MBoemo / DNAscent

Software for detecting regions of BrdU and EdU incorporation in Oxford Nanopore reads.
https://www.boemogroup.org/
GNU General Public License v3.0
26 stars 13 forks source link

Resolution Issue #11

Closed henricheng closed 3 years ago

henricheng commented 3 years ago

Hi, I'm using version 2. The resolution keeps on being 1000 even when I change it to -r 250 when using regions. I'm not sure what the issue is. Is there another option I need to include?

MBoemo commented 3 years ago

I just took a look and wasn't able to reproduce this issue - can you share the full command you ran and a few lines from the output file so we can get this resolved?

henricheng commented 3 years ago

DNAscent detect -b alignment.sorted.bam -r w303.fasta -i index.dnascent -o output.detect -t 200 --GPU 0 DNAscent regions -d output.detect -o output.regions -r 250

MBoemo commented 3 years ago

Command looks fine, can I see a few lines from output.regions where you're having the issue? Also, does this problem persist on every region across all reads, or just a few? If you're looking in a very GC-rich region, note this from the documentation:

Note that the region width may sometimes vary slightly from the value specified. The region width is designated as the coordinate of the first thymidine greater than the window width (100 bp by default) from the starting coordinate. In order to guard against assigning a score to regions with very few thymidines, DNAscent regions will also extend the region until at least 10 calls are considered.

henricheng commented 3 years ago

The regions script isn't an issue per say but I just used the -l 250 on detect and it fixed the program.

henricheng commented 3 years ago

I'll getting this error when I change the --length for the detect portion less than 250.

DNAscent: src/alignment.cpp:1600: std::pair<bool, std::shared_ptr > eventalign_detect(read&, unsigned int, double): Assertion `scaledEvent > 0.0' failed.

MBoemo commented 3 years ago

I haven't managed to reproduce this issue yet (thus far I've tried it on a few thousand 50 bp and 100 bp reads and everything worked as expected) so I'm going to need a little more information. Does this happen right away (e.g., on the first read you analyse) or does it make it through some amount of reads first before the issue crops up?

henricheng commented 3 years ago

I don't know the exact number of reads but it happens after a few hours into the detect portion.

MBoemo commented 3 years ago

Okay thanks. My best guess at this point is that there's something unusual about the signal from a particular read (probably something rare, hence it's hard to reproduce) that's not being caught correctly. Personally I rarely run detect on reads less than 10 kb or so, so it's possible that some sort of rare issue specifically pertaining to short reads slipped under the radar. Might take me a day or two to get it sorted, but I'll give a shout with either a resolution or if I need extra information.

henricheng commented 3 years ago

I'm not sure if this has anything to do with it, but for this specific run, it is creating a temporary file when doing the detect program that didn't do the same thing with the other runs.

MBoemo commented 3 years ago

The temporary file does sound unusual - do you have any more information about it? Did you analyse any other runs with similar length where this wasn't an issue?

Unfortunately I still haven't managed to reproduce this, but I did make you a debugging build so we can get to the bottom of it. Please can you do the following:

git pull git checkout debugging make clean && make

Make sure you're on commit fdfabb8 (git log -1) and then use DNAscent detect on the run that's causing the issue, but redirect stderr to a file. Copying your command from above, this would look something like:

DNAscent detect -b alignment.sorted.bam -r w303.fasta -i index.dnascent -o output.detect -t 200 --GPU 0 2> debugging.stderr

With this build, some diagnostic information about the read should be printed to stderr right before it trips on the problematic read. If you either attach debugging.stderr here or Email it to me, we should be able to figure out what's going on.

MBoemo commented 3 years ago

For the git commands, just navigate to the DNAscent directory and run them there. Make sure you run it with whatever length options caused the issue though.