I have a Python program to find mutations in my sequence and load them into a dictionary. I search both strands by comparing the complement mutation to my mutation list. I use comp = Seq(str(ref_seq+alt_seq)).complement()
from Biopython to generate the complement. When I run this on a box that only has Python 3.3.3 available I get thousands fewer mutation than when I run with Python 3.4. It is the Biopython complement function that is causing this because I wrote my own function to make the complements and the error is gone. Has anyone seen this before? The results are frighteningly different. The highlighted line of code is the Biopython complement that gives errors. I have also listed my code for generating the complement. The burning is, "is this a bug?" The error is not random but does not occur with each call, maybe 3% of the time.
elif ref[i+4] != qseq[i] and qseq[i] != "N" and qseq[i] != ".":
segment_mut_count += 1 # Base substitutions in this segment.
interval_mutation_count += 1
alt_seq = qseq[i]
ref_seq = ref[i+4]
minus3 = ref[i+1:i+4]
plus3 = ref[i+5:i+8]
read_mutation_data.append(["mut", a.reference_start + i, ref_seq, alt_seq, minus3, plus3, mut_dep])
if "M" in o.output and int(o.matrix_type) != 3:
qmutation = ref_seq+alt_seq
comp = complement(ref_seq+alt_seq)
`#comp = Seq(str(ref_seq+alt_seq)).complement()`
for key in key1: # Count mutation patterns in the query sequence.
if qmutation == key:
mutation_matrix[qmutation][key2][key3][0] += 1
elif comp == key:
try:
mutation_matrix[comp][key2][key3][0] += 1
except TypeError:
pass
def complement(seq):
comp_seq = ''
for I in range(len(seq)):
if seq[i] is "G":
nuc = "C"
elif seq[i] is "A":
nuc = "T"
elif seq[i] is "T":
nuc = "A"
elif seq[i] is "C":
nuc = "G"
else:
nuc = "N"
comp_seq += nuc
return comp_seq
Hello dennis!
It appears that your post has been cross-posted to another site: Duplicate of issue opened about same time at Biopython
https://github.com/biopython/biopython/issues/518
This is typically not recommended as it runs the risk of annoying people in both communities.
Sorry. When I tracked it into Biopython I thought I better post it as a bug closer to the source.