LRSIM hangs during manifest generation step

[aquaskyline]

LRSIM hangs during manifest generation step

[CodingKaiser]

I think we realized what the issue is. It would appear that the hanging can be traced back to the loop on lines 509-520 and the fact that, if we run out of barcodes, there is no check to see if we have run out of barcodes, resulting in the program helplessly searching for an available barcode in the same empty barcodes list forever. The number of barcodes in my list of parameters did not exceed the 4.7M limit however, so it was surprising that we were running out. The culprit ended up being the size of the input genome, which was ~160M bp in length, and the fact that I was asking for 400M reads to be simulated.

This was a problem because it would appear that LRSim only uses as many of the simulated DWGSIM reads as there are positions in the genome, storing a mapping of genomic locations to reads in the .fp file (see readToPos and writeToPos). In the whole genome case, this is not an issue since the genome is about 3B in size and most whole genome simulation with 30X coverage requires only about 600M reads, so there are enough reads to go around. In my case however, since I subsampled the human genome to a total of ~160M bp, there weren't enough reads to meet the ~400M number I was feeding to the -x parameter.

Here is how we think this is resulting in the infinite-loop hanging I mentioned earlier. The first thing the loop on line 553 of the main program does is look for a read position at a given genomic location. If all of the reads in the immediate vicinity of the given genomic location have already been used, we skip this read and go to the next read (lines 555-556). Since virtually all of the available reads are exhausted once my simulation reaches close to 160M reads, almost every iteration of this loop will mean we don't decrement the read counter but do erase the barcode from the available barcodes list, therefore quickly exhausting all remaining barcodes without actually assigning any reads to them. Once all barcodes are used up in this manner, the loop on lines 509-520 will run infinitely.

Introducing a check such as the following in lines 504-521 at least lets the user know that we have run out of barcodes and exits rather than looping infinitely:

#Pick a barcode
my $selectedBarcode;
{
  my $idx = int(rand($numBarcodes));
  lock($barcodesMutexLock);
  my $wentToZero = 0;
  while(1)
  {
    if($barcodes[$idx] eq "")
    {
      ++$idx;
      if($idx == $numBarcodes && not($wentToZero)) {
          $idx = 0;
          $wentToZero = 1;
      } elsif($idx == $numBarcodes && $wentToZero) {
          &LogAndDie("Reached end of barcodes list. No more barcodes. Last read processed: $readsCountDown. Exiting.");
      }
      next;
    }
    $selectedBarcode = $barcodes[$idx];
    $barcodes[$idx] = "";
    last;
  }
}

It would also be great if the number of reads wasn't limited to the size of the input genome, but that would probably require a more extensive rewrite.

[aquaskyline]

@CodingKaiser thanks for the detailed analysis, that helps other users a lot. For the fix, do you want me to do the integration or you want to make a pull request? The latter better records your contribution to the code.

[CodingKaiser]

Thanks for getting back to me. I can submit a pull request shortly. Is there any chance that the rest of the code will be rewritten to remove the limitation of 1 read per genomic location?

[aquaskyline]

That requires a revamp algorithmic-wise. Would it possible for you to superimpose two simulated datasets to imitate reads duplication?