<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>

CC BY-SA 4.0
2021-01-27T07:19:08.207Z
2021-01-27T07:19:08.207Z
Liam McIntyre
<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>

CC BY-SA 4.0
2021-01-26T04:47:47.017Z
Do you have any additional assumptions you can make? For example, are the pairs of barcodes within a specific distance from each other  or the reads ends? I ask because with an edit distance of two (including indels) the chances of random matches will start to increase
GWW
<phylogenetics><fastq><long-reads><bioawk><barcode>
Find pattern that is present twice and allow <=2 mismatches on each
# Find pattern that is present twice and allow <=2 mismatches on each

(Asked by: username_0 on 2021-01-25)

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-

> TATCTTGTGGAAAGGACGAAACACCGAACACAAAGCATAGATGCGTTTAAGAGCTATGCTGGAAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTT**ATTCGACCGATAGG**GGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAG**ATTCGACCGATAGG**TGGCGTAACTAGATCTTGAGACAAA
> TATCTTGTGGAAAGGACGAAACACCGGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTT**ATTCGACCGATAGG**GGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAG**ATTCGACCGATAGG**TGGCGTAACTAGATCTTGAGACAAA
> TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTT**ATTCGACCGATAGG**GGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAG**ATTCGACCGATAGG**TGGCGTAACTAGATCTTGAGACAAA
> TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTT**ATTCGACGATAGG**GGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAG**ATTCGACCGATAGG**TGGCGTAACTAGATCTTGAGACAAA

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.

**Tags:** phylogenetics, fastq, long-reads, bioawk, barcode
> **Comments:**
> **username_1**: Do you have any additional assumptions you can make? For example, are the pairs of barcodes within a specific distance from each other  or the reads ends? I ask because with an edit distance of two (including indels) the chances of random matches will start to increase  

---

## 1 answers

### Answer by username_2 on 2021-01-27. Score: 1

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.


---

