Question: Can anybody work out why my python script is not working properly? I want to search through a FASTA looking for a specific sequence (**CAGTGGG..GCAA[TA]AA** and the reverse complement), then print the location of the sequence. However, the output only prints the location of the "forward" sequence and the location is incorrect (in both examples I'm sure it should be 2/3) and doesn't seem to identify both locations. I can't work out why. So far I have the following code:

#!/bin/python

#My Mimp_finder

import re

from Bio import SeqIO

for seq in SeqIO.parse("Focub_mimps.fasta", "fasta"):
print("\n\n")
print(seq)
hit = re.search(r"CAGTGGG..GCAA[TA]AA", str(seq))
hit_rc = re.search(r"TT[TA]TTGC..CCCACTG", str(seq))
if hit and hit_rc:
print(f'complete mimp found at {hit.start()} to {hit_rc.end()}')
elif hit:
print(f'Partial forward mimp found at {hit.start()}')
elif hit_rc:
print(f'Partial rc mimp found at {hit_rc.end()}')
else:
print("No Mimps")

Here is an example of the FASTA (matching sequences in bold):

>Focub_II5_mimp_1__contig_1.16(656599:656809)
TA**CAGTGGGATGCAAAAA**GTATTCGCAGGTGTGTAGAGAGATTTGTTGCTCGGAAGCTAGTTAGGTGTAGCTTGTCAGGTTCTCAGTACCCTATATTACACCGAGATCAGCGGGATAATCTAGTCTCGAGTACATAAGCTAAGTTAAGCTACTAACTAGCGCAGCTGACACAACTTACACACCTGCAAATAC**TTTTTGCATCCCACTG**TA
>Focub_II5_mimp_2__contig_1.47(41315:41540)
TA**CAGTGGGAGGCAATAA**GTATGAATACCGGGCGTGTATTGTTTTCTGCCGCTAGCCCATTTTAACAGCTAGAGTGTGTATATTAACCTCACACATAGCTATCTCTTATACTAATTGGTTAGGGAAAACCTCTAACCAGGATTAGGAGTCAACATAGCTTGTTTTAGGCTAAGAGGTGTGTGTCAGTACACCAAAGGGTATTCATAC**TTATTGCCCCCCACTG**TA

The output for these sequences:

ID: Focub_II5_mimp_1__contig_1.16(656599:656809)
Name: Focub_II5_mimp_1__contig_1.16(656599:656809)
Description: Focub_II5_mimp_1__contig_1.16(656599:656809)
Number of features: 0
Seq('TACAGTGGGATGCAAAAAGTATTCGCAGGTGTGTAGAGAGATTTGTTGCTCGGA...GTA', SingleLetterAlphabet())
Partial forward mimp found at 187

ID: Focub_II5_mimp_2__contig_1.47(41315:41540)
Name: Focub_II5_mimp_2__contig_1.47(41315:41540)
Description: Focub_II5_mimp_2__contig_1.47(41315:41540)
Number of features: 0
Seq('TACAGTGGGAGGCAATAAGTATGAATACCGGGCGTGTATTGTTTTCTGCCGCTA...GTA', SingleLetterAlphabet())
Partial forward mimp found at 181

Why are only the "forward" mimp locations printed?

P.S. any recommendations for improvement or helpful tips will be greatly appreciated - I am just starting out.

Thanks,

Jamie

Answer: ##NB, this post is heavily edited from its original content due to length.

I've revised the code with the clarification from the OP, and I think this (and their) code should work - though not for the reasons they thought.

import re, sys
from Bio import SeqIO

def red(text):
return "\033[0;31m{}\033[0m".format(text)

def blue(text):
return "\033[0;34m{}\033[0m".format(text)

regex = re.compile("CAGTGGG..GCAA[TA]AA")
regex_rc = re.compile("TT[TA]TTGC..CCCACTG")

for record in SeqIO.parse(sys.argv[1], "fasta"):
hit = re.search(regex, str(record.seq))
hit_rc = re.search(regex_rc, str(record.seq))

print("="*80 + "\n", record)
if hit and hit_rc:
record.seq = str(record.seq).replace(hit.group(0), red(hit.group(0))).replace(hit_rc.group(0), blue(hit_rc.group(0)))
print(f"MIMPs found on forward strand at {hit.start()}:{hit.end()}, and reverse strand at {hit_rc.start()}:{hit_rc.end()}")
print(record.seq)
elif hit and not hit_rc:
record.seq = str(record.seq).replace(hit.group(0), red(hit.group(0)))
print(f"Partial MIMP found on forward strand only at {hit.start()}:{hit.end()}.")
print(record.seq)
elif hit_rc and not hit:
record.seq = str(record.seq).replace(hit_rc.group(0), blue(hit_rc.group(0)))
print(f"Partial MIMP found on reverse strand only at {hit_rc.start()}:{hit_rc.end()}")
print(record.seq)
elif not hit and not hit_rc:
print(f"No MIMPs found on either strand")

I still think a more robust method might be to search the reverse strand specifically and adjust the coordinates relative to the forward strand, but the recognition sequence is long enough that the odds of it appearing by chance are relatively low I think.

For clarity of output, the code above colours the result:

I believe this is correct. The code to colour the matches can be thrown away. It seems my previous code searching the reverse strand was correct, and the hits coming up at 2:18 all the time, was because the target sites are (in this test data) always 2 bases in from the 5' end of the sequence.

#### At least for this sample set, searching the forward strand with the regex and its reverse complement appears to work satisfactorily, however it is worth being aware that the greedy/conservative nature of regex capture will affect which sequence is captured in the event that there are more than one per query sequence, since the current regex strategy does not use `findall`.