How can I make this Biopython program (to correct erroneous barcodes) run faster, and is there any alternative method?

Question

This question has also been asked on Biostars

I am looking forward to getting a valuable suggestion for a bioinformatic problem.

Background:

Currently, I am performing a de novo whole genome assembly. At the stage of barcode correction, I lost nearly half of all the reads due to erroneous barcodes.

Barcodes, in this context, are 18-nucleotide-long DNA sequences that were added to the DNA insert during NGS library preparation. Now, I have a file that contains all the 18-base barcodes used in my projects. These barcodes are important because they are unique for reads that came from a DNA molecule and therefore helpful in generating long-range information.

While library construction, 18-base molecular barcodes were used for the linked-read generation, and the barcodes have a high incidence of "N" in positions between 12 and 18, as you may see below.

There is a very large occurrence of Ns in positions 12 to 18 of the barcodes and for this reason, a large pool of erroneous barcodes were filtered.

Now, I have been asked by the company to do one of two things, to recover the lost reads:

1. Either replace the barcodes(barcode fastq file from sequencing core) with the correct ones from a pool of barcodes (This is a text file I received from the company that contains the pool of all 18 base barcode sequences)

2. Since most "N"s are seen in positions 12 to 18 of the 18-base barcode sequences (barcode fastq file from the sequencing core), replace the 11-mers in the barcode file, and replace them with the 11-mers from the correct barcode file (text file).

I have decided to use Python code to execute the first method.

#Define the path where fastq file is stored.
import os
from Bio.SeqIO.QualityIO import FastqGeneralIterator 
path=os.path.join(os.path.expanduser("~"), "Desktop", "SVA1_S1_R2_001.fastq")

#create a new file handle to write the file
path_new=os.path.join(os.path.expanduser("~"), "Desktop", "trimmedfastq.fastq")
path_ust_barcode=os.path.join(os.path.expanduser("~"), "Desktop", "UST40MBarcodes.txt")
write_file=open(path_new, "w")

#open the text file containg the correct pool of barcodes

with open(path_ust_barcode, "r") as fh:
   barcodes = fh.readlines()

handle=open(path, "r")  
output = []

#===================================================
# Return the Hamming distance between string1 and string2.
def hamming_distance(string1, string2): 
    # Start with a distance of zero, and count up
    distance = 0
    # Loop over the indices of the string
    L = len(string1)
    for i in range(L):
        # Add 1 to the distance if these two characters are not equal
        if string1[i] != string2[i]:
            distance += 1
    # Return the final count of differences
    return distance
#===================================================

#set correct barcode to "None"
corr_barcode=None
#set the minimum hamming distance to "None"
min_distance=19
#open the fastq file

for title, sequence, qual in FastqGeneralIterator(handle):

    # if a perfect match exists, just output as is
    if sequence in barcodes:
        output.append("@%s\n%s+\n%s\n" % (title, corr_barcode, qual))
        continue

    for read in barcodes:

        h_distance=hamming_distance(sequence, read)

        # since if there was a perfect match we won't be here, 
        # if distance is 1, we'll never do better
        if h_distance == 1:
            corr_barcode = read
            break

        if  h_distance <= min_distance:
            min_distance=h_distance 
            corr_barcode=read

    output.append("@%s\n%s+\n%s\n" % (title, corr_barcode, qual))


with open(path_new, "a") as wf:
    for line in output:
        wf.write(line)

handle.close()
print("Done")

Problem

I would like to know if this is the only Pythonic way to tackle this problem. This method is dead slow. My erroneous barcode fastq file is 30+ Gb, and the UST10x barcode file (the text file that contains all the correct barcode pool) is of size 750Mb. I have split this text file into 1000 smaller files, each about 750kb and having ~40,000 lines per file. It's been the third day and the program is still running without outputting anything. I am looking forward to getting valuable suggestions here, with respect to the code, especially how I can make the process faster.

NB: I've tried replacing

with open(path_ust_barcode, "r") as fh:
   barcodes = fh.readlines()

with

barcodes = list()
with open(path_ust_barcode, "r") as fh:
   for b in fh:
      barcodes.append(b.rstrip())

Still, the code runs extremely slowly.

I am using an IBM server with a 197GM RAM, and an octa-core processor. In my initial code, it was iterating through all the correct barcodes of the text file (UST10x file), for each erroneous barcode in the fastq file, at the end of every iteration, it opened a text file in append mode and wrote the correct version of fastq reads. Since my server has high memory capacity, I rewrote the code in such a way that it stored all the correct barcodes in the memory, and finally wrote it in a text file, by bringing the write function outside of the loop.

Python is not absolutely necessary. I'm not familiar with C or C++, but that being said, I would really love to get this process done. I'm using an IBM server with Ubuntu (Jammy), and 197GB Ram, and the current local storage has only 150GM free space. if you can provide me suggestions as to how to execute it in C or C++, I can do it.

Answer 1

Try profiling (benchmarking) your code to work out where it is slow.

My other suggestions are:

• replace your simple Pure Python hamming_distance function with something optimised like RapidFuzz - see https://github.com/maxbachmann/RapidFuzz

• cache the corrections you are making - you have tonnes of RAM so a building a dictionary mapping the fragments with N to the corrected fragment ought to be possible. And then for each new sequence, check if you already solved it and if so reuse the same choice.