Answers.Body: <p>Try profiling (benchmarking) your code to work out where it is slow.</p>
  <p>My other suggestions are:</p>
  <ul>
  <li><p>replace your simple Pure Python hamming_distance function with something optimised like RapidFuzz - see <a href="https://github.com/maxbachmann/RapidFuzz" rel="nofollow noreferrer">https://github.com/maxbachmann/RapidFuzz</a></p>
  </li>
  <li><p>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.</p>
  </li>
  </ul>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-05-06T07:13:32.450Z
Answers.Text: Thank you for adding a response. I have posted the same question on Biostar, where you can find the code aligned with your logic. A modified version of the code was suggested by @IanSudbery, and you may take a look at it. https://www.biostars.org/p/9561109/#9561603
Answers.UserDisplayName: VIJITH KUMAR
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-05-04T21:16:49.017Z
Answers.LastActivityDate: 2023-05-04T21:16:49.017Z
Answers.UserDisplayName: Peter Cock
Body: <p><em>This question has also been asked on <a href="https://www.biostars.org/p/9561109" rel="nofollow noreferrer">Biostars</a></em></p>
  <p>I am looking forward to getting a valuable suggestion for a bioinformatic problem.</p>
  <p><strong>Background:</strong></p>
  <p>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.</p>
  <p>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.</p>
  <p><a href="https://i.sstatic.net/hdwIJ.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/hdwIJ.png" alt="Barcode statistics" /></a></p>
  <p>While library construction, 18-base molecular barcodes were used for the linked-read generation, and the barcodes have a high incidence of &quot;N&quot; in positions between 12 and 18, as you may see below.</p>
  <p><a href="https://i.sstatic.net/W70ih.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/W70ih.png" alt="Erroneous barcodes" /></a></p>
  <p>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.</p>
  <p>Now, I have been asked by the company to do one of two things, to recover the lost reads:</p>
  <ol>
  <li><p>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)</p>
  </li>
  <li><p>Since most &quot;N&quot;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).</p>
  </li>
  </ol>
  <p>I have decided to use Python code to execute the first method.</p>
  <pre class="lang-py prettyprint-override"><code>#Define the path where fastq file is stored.
  import os
  from Bio.SeqIO.QualityIO import FastqGeneralIterator 
  path=os.path.join(os.path.expanduser(&quot;~&quot;), &quot;Desktop&quot;, &quot;SVA1_S1_R2_001.fastq&quot;)
  
  #create a new file handle to write the file
  path_new=os.path.join(os.path.expanduser(&quot;~&quot;), &quot;Desktop&quot;, &quot;trimmedfastq.fastq&quot;)
  path_ust_barcode=os.path.join(os.path.expanduser(&quot;~&quot;), &quot;Desktop&quot;, &quot;UST40MBarcodes.txt&quot;)
  write_file=open(path_new, &quot;w&quot;)
  
  #open the text file containg the correct pool of barcodes
  
  with open(path_ust_barcode, &quot;r&quot;) as fh:
     barcodes = fh.readlines()
  
  handle=open(path, &quot;r&quot;)  
  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 &quot;None&quot;
  corr_barcode=None
  #set the minimum hamming distance to &quot;None&quot;
  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(&quot;@%s\n%s+\n%s\n&quot; % (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 &lt;= min_distance:
              min_distance=h_distance 
              corr_barcode=read
  
      output.append(&quot;@%s\n%s+\n%s\n&quot; % (title, corr_barcode, qual))
  
  
  with open(path_new, &quot;a&quot;) as wf:
      for line in output:
          wf.write(line)
  
  handle.close()
  print(&quot;Done&quot;)
  </code></pre>
  <p><strong>Problem</strong></p>
  <p>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.</p>
  <p>NB:  I've tried replacing</p>
  <pre class="lang-py prettyprint-override"><code>with open(path_ust_barcode, &quot;r&quot;) as fh:
     barcodes = fh.readlines()
  </code></pre>
  <p>with</p>
  <pre class="lang-py prettyprint-override"><code>barcodes = list()
  with open(path_ust_barcode, &quot;r&quot;) as fh:
     for b in fh:
        barcodes.append(b.rstrip())
  </code></pre>
  <p>Still, the code runs extremely slowly.</p>
  <p>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.</p>
  <p>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.</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-04-25T23:53:53.933Z
Comments.Text: **Comments have been [moved to chat](https://chat.stackexchange.com/rooms/145591/discussion-on-question-by-vijith-kumar-how-can-i-make-this-biopython-program-to); please do not continue the discussion here.** Before posting a comment below this one, please review the [purposes of comments](/help/privileges/comment). Comments that do not request clarification or suggest improvements usually belong as an [answer](/help/how-to-answer), on [meta], or in [chat]. Comments continuing discussion may be removed.
Comments.UserDisplayName: gringer
Tags: <sequence-analysis><biopython>
Title: How can I make this Biopython program (to correct erroneous barcodes) run faster, and is there any alternative method?
ThreadText: # How can I make this Biopython program (to correct erroneous barcodes) run faster, and is there any alternative method?
  
  (Asked by: username_0 on 2023-04-25)
  
  *This question has also been asked on [Biostars](https://www.biostars.org/p/9561109)*
  
  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.
  
  [![Barcode statistics](https://i.sstatic.net/hdwIJ.png)](https://i.sstatic.net/hdwIJ.png)
  
  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.
  
  [![Erroneous barcodes](https://i.sstatic.net/W70ih.png)](https://i.sstatic.net/W70ih.png)
  
  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.
  
  ```py
  #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
  
  ```py
  with open(path_ust_barcode, "r") as fh:
     barcodes = fh.readlines()
  
  ```
  
  with
  
  ```py
  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.
  
  **Tags:** sequence-analysis, biopython
  
  > **Comments:**
  > **username_1**: **Comments have been [moved to chat](https://chat.stackexchange.com/rooms/145591/discussion-on-question-by-username_0-username_0-how-can-i-make-this-biopython-program-to); please do not continue the discussion here.** Before posting a comment below this one, please review the [purposes of comments](/help/privileges/comment). Comments that do not request clarification or suggest improvements usually belong as an [answer](/help/how-to-answer), on [meta], or in [chat]. Comments continuing discussion may be removed.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2023-05-04. Score: 2
  
  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.
  
  > **Comments:**
  > **username_0**: Thank you for adding a response. I have posted the same question on Biostar, where you can find the code aligned with your logic. A modified version of the code was suggested by @IanSudbery, and you may take a look at it. https://www.biostars.org/p/9561109/#9561603  
  
  ---
  
  