This is a test version of Biostars. For the public version, visit https://www.biostars.org.
string search in fastq file

I have a list of 20 bp sequences stored in a list called known. I have a fastq file called A.fastq, from which I cut out the relevant 20 bp region and compare against the known list. I would like to use Hamming distance as a metric for choosing the correct match i.e for every fastq read, I want to find the string in the known list which has is the least Hamming distance apart and update a counter good_count only if the least Hamming distance is <= 3. This is my code, but this does not account for Hamming distance criteria. Can you suggest what I should do?

import Bio
from Bio import SeqIO

known = set()
for s in Bio.SeqIO.parse("lib.fa","fasta"):
        known.add(str(s.seq[10:20]))

good_count = 0
for r in Bio.SeqIO.parse("read1.fastq","fastq"):
        if str(r.seq[10:20]) in known:
                good_count += 1
print good_count

Specifically, I want help at the "if" statement to compare string by string from the list "known"

sequencing
• 2,557 views • link
updated by Alex Reynolds Alex Reynolds • written by sumithrasank75 sumithrasank75

Look into the python distance package.

• 0 views • link
written by st.ph.n st.ph.n

1 answer

The itertools.imap() method runs the operator.ne function on each pair of characters in the strings a and b. The sum of mismatches (True -> 1) and matches (False -> 0) gives the distance score.

import itertools
import operator

def hamming(a, b):
    return sum( itertools.imap( operator.ne, a, b ) )
• 0 views • link
written by Alex Reynolds Alex Reynolds

Log in to answer this question.

More posts like this