text: user1911: I have a fastq file of 400,000 reads (so speed is important). In the sequences there are barcodes integrated that should be present twice. Given a barcode, I want to find the sequences that have the barcode present twice with <= 2 mismatches. So, with a barcode 'ATTCGACCGATAGG', I would like to retrieve all of the following sequences-
  
  TATCTTGTGGAAAGGACGAAACACCGAACACAAAGCATAGATGCGTTTAAGAGCTATGCTGGAAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACCGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA
  TATCTTGTGGAAAGGACGAAACACCGGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACCGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA
  TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACCGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA
  TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA
  
  Note that the first barcode in the fourth sequence is short of one character. I have tried with biopython and regex but it's just too slow given I have 5k barcodes. I am wondering if there is a fast solution available in python or in something like grep, awk or anything else. Thanks.
  
  user7153: This isn't perfect but you can build on it to deal with any issues (ie overalpping hits). Also if you need it faster just chunk it and run in parallel threads.
  I just copy pasted your examples seqs into files to read (called them seqs and code).
  #need blast installed 
  import subprocess
  import pandas as pd
  
  path = '/Users/ur/dir/'
  #make a blast db
  with open (path + 'seqs.txt', 'r') as fin:
      with open (path + 'seqs.fa', 'w') as fout:
          for i, line in enumerate(fin):
              fout.write('>'+str(i)+'\n')
              fout.write(line)
              
  subprocess.run([
      'makeblastdb',
      '-in', path + 'seqs.fa',
      '-dbtype', 'nucl'
  ])
  
  #use blast
  subprocess.run([
      'blastn',
      '-task', 'blastn-short',
      '-query', path + 'code.txt',
      '-db', path + 'seqs.fa',
      '-out', path + 'out.txt',
      '-outfmt', '6'
  ])
  
  #grab seqs with 2 hits (ie ~2x len(query)
  df = pd.read_csv(
      path + 'out.txt',
      sep='\t',
      header=None,
      names=['query', 'db_name', 'percent','len'],
      usecols = [0,1,2,3]
  )
  df = df.groupby('db_name').agg('sum')
  print(*list(df.query('len > 27').index))
  
  I this case they are all hits. If you look at the blast output you'll see one hit was fragmented (blast isn't perfect with small seqs). But total len works to resolve this.
  
original_text: 1911: <p>I have a fastq file of 400,000 reads (so speed is important). In the sequences there are barcodes integrated that should be present twice. Given a barcode, I want to find the sequences that have the barcode present twice with &lt;= 2 mismatches. So, with a barcode 'ATTCGACCGATAGG', I would like to retrieve all of the following sequences-</p>
  <blockquote>
  <p>TATCTTGTGGAAAGGACGAAACACCGAACACAAAGCATAGATGCGTTTAAGAGCTATGCTGGAAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTT<strong>ATTCGACCGATAGG</strong>GGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAG<strong>ATTCGACCGATAGG</strong>TGGCGTAACTAGATCTTGAGACAAA
  TATCTTGTGGAAAGGACGAAACACCGGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTT<strong>ATTCGACCGATAGG</strong>GGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAG<strong>ATTCGACCGATAGG</strong>TGGCGTAACTAGATCTTGAGACAAA
  TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTT<strong>ATTCGACCGATAGG</strong>GGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAG<strong>ATTCGACCGATAGG</strong>TGGCGTAACTAGATCTTGAGACAAA
  TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTT<strong>ATTCGACGATAGG</strong>GGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAG<strong>ATTCGACCGATAGG</strong>TGGCGTAACTAGATCTTGAGACAAA</p>
  </blockquote>
  <p>Note that the first barcode in the fourth sequence is short of one character. I have tried with biopython and regex but it's just too slow given I have 5k barcodes. I am wondering if there is a fast solution available in python or in something like grep, awk or anything else. Thanks.</p>
   7153: <p>This isn't perfect but you can build on it to deal with any issues (ie overalpping hits). Also if you need it faster just chunk it and run in parallel threads.</p>
  <p>I just copy pasted your examples seqs into files to read (called them seqs and code).</p>
  <pre><code>#need blast installed 
  import subprocess
  import pandas as pd
  
  path = '/Users/ur/dir/'
  #make a blast db
  with open (path + 'seqs.txt', 'r') as fin:
      with open (path + 'seqs.fa', 'w') as fout:
          for i, line in enumerate(fin):
              fout.write('&gt;'+str(i)+'\n')
              fout.write(line)
              
  subprocess.run([
      'makeblastdb',
      '-in', path + 'seqs.fa',
      '-dbtype', 'nucl'
  ])
  
  #use blast
  subprocess.run([
      'blastn',
      '-task', 'blastn-short',
      '-query', path + 'code.txt',
      '-db', path + 'seqs.fa',
      '-out', path + 'out.txt',
      '-outfmt', '6'
  ])
  
  #grab seqs with 2 hits (ie ~2x len(query)
  df = pd.read_csv(
      path + 'out.txt',
      sep='\t',
      header=None,
      names=['query', 'db_name', 'percent','len'],
      usecols = [0,1,2,3]
  )
  df = df.groupby('db_name').agg('sum')
  print(*list(df.query('len &gt; 27').index))
  </code></pre>
  <p>I this case they are all hits. If you look at the blast output you'll see one hit was fragmented (blast isn't perfect with small seqs). But total len works to resolve this.</p>
  