This question has also been asked on Biostars
I am looking forward to getting a valuable suggestion for a bioinformatic problem.
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:
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)
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")
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()
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.