1
$\begingroup$

I am a college student and just started bioinformatics. I am trying to write a script (not for school) that finds all potential NGG and TTTN protospacer adjacent motif sequences in a genome string. I wrote the following script:

TargetFile = open('file.fasta')
readTF = TargetFile.read()
useTF = "".join(readTF[66:].split())


#no_of_basepairs = len(useTF)


sites = {"ngg_strand" : [],
         "ngg_antistrand" : [],     #location counted from reverse direction
         "tttn_strand" : [],
         "tttn_antistrand" : []}    #location counted from reverse direction


for i in range(0, len(useTF)-2):
    gg_mer = useTF[i:i + 2]
    ttt_mer = useTF[i:i + 3]

    if gg_mer == "GG":
        sites["ngg_strand"].append(i)
    elif gg_mer == "CC":
        #look for reverse complement, record location on antistrand from reverse direction
        sites["ngg_antistrand"].append(len(useTF)-i-2)
    elif ttt_mer == "TTT":
        sites["tttn_strand"].append(i+1)
    elif ttt_mer == "AAA":
        # look for reverse complement, record location on antistrand from reverse direction
        sites["tttn_antistrand"].append(len(useTF)-i-2)


print(len(sites["ngg_strand"]))
print(useTF.count("GG"))

At the bottom I used the .count() method to compare results. I noticed that the method returns a smaller number than my code. I believe this is due to the .count() method not counting overlap, whereas my code does. I.e. in the string "NNNGGGGNNN" my code would count 3 occurences of the "GG" substring, whereas the .count() method returns 2 occurences.

My Questions:

  1. If I want to find all possible PAM sites, which way of counting is correct? I believe my way would be correct, since as far as I know the RNA molecule doesn't advance in "triplet-steps". However, I am very much new to this and might be missing something.

  2. If the reverse complement is found, is the location usually recorded from the same direction as on the strand, or in the reverse direction (as I implemented in the code)?

  3. Like I said, I am new to this so if you see any other glaring problems I would appreciate any and all help. :-)

Thank you!

$\endgroup$
6
  • $\begingroup$ Hi @heibai thanks for your question. If you have a moment could you please clarify "NGG and TTTN protospacer"? Its just question 2 is difficult to understand and I suspect relates to the definitinon of NGG and TTTN protospacer motifs. Do you mean in question 2 that AAA and TTT are reverse complements? $\endgroup$
    – M__
    Nov 3 '20 at 10:33
  • $\begingroup$ Hi @Michael! So I am trying to find all possible editing sites for Cas9 and Cas12, which in this case have PAM sequences of NGG and TTTN. As far as I understand the protein complex reads the strands in opposite directions, meaning that if the .fasta file contains a substring "NAAA", the antistrand at that position would be "NTTT". Read in the reverse direction it becomes "TTTN", i.e. a potential editing site. Since the "N" can be any nucleotide, I am not including it in the search. I hope this clarifies things and doesn't just further confuse what I mean. :-D $\endgroup$
    – heibai
    Nov 3 '20 at 10:49
  • $\begingroup$ @Michael If you have time could you expand on that? What kind of situations would result in false hits? What do you mean by defining forward and reverse direction, and how would I do that? $\endgroup$
    – heibai
    Nov 3 '20 at 11:36
  • $\begingroup$ If this is not an assignment but for productive use then use specialized software such as ChopChop of CCTop (or many others) as you also have to check for genome-wide similarity regions (off-targets) plus there are many factors that can influence efficacy of the gRNA to actually recruit the Cas protein to the target site which is implemented in these tools. It is not as simple as a mattern match. For example: cctop.cos.uni-heidelberg.de:8043 It has a stand-alone version. $\endgroup$
    – ATpoint
    Nov 3 '20 at 12:13
  • $\begingroup$ Okay @Michael, thank you for your help. $\endgroup$
    – heibai
    Nov 3 '20 at 15:35
1
$\begingroup$

The easiest solution is finditer

use re

# Compile you search term where [A|T|G|C] means A or G or C or T 
p = re.compile("[A|T|G|C]GG")

iterator = p.finditer(useTF)
for match in iterator:
    print(match.span())

** Example Output**

(0, 2)
(22, 24)
(29, 31)

Easy does it. Its worth looking at the documentation for re here, https://docs.python.org/3/howto/regex.html for alternative approaches using match or search below. finditer is definitely the solution in context.


I support your coding approach because it provides lots of flexibility. For example, you might want to target GG or TTT around a given locus and you can easily modify the code to accomodate this. Someone elses code would be less flexible.

In summary, I kinda like your substrings and do like the count method. I'm not a fan of your == in this instance and your C-style matching should be substituted for reg-ex which in Python is implemented via re.

The basic issues are,

  • Biopython would be used via AlignIO or SeqIO. This will allow any given format to be analysed. Yours approach is okay but prone to bugs (below). [probably a minor issue]
  • Your sequence is in some sort of weird format, like extended phylip format and it is not fasta (that wouldn't work). It needs verifying because removing the first 65 characters might remove bits of sequence [probably a minor issue]
  • Chopping of -2 from the length ... I dunno I'd check out -1 [minor]
  • C-style matching. I don't like this. I would definitely use re. One solution uses match or search, but the solution above uses finditer, Reg-ex will allow to search for NGG and TTTN, which you are not doing
  • I think you are looping a locus of 3 and 4 bp and I'm not sure how the C-style == would match this because thats 2 and 3 bp respectively. Just thoughts. You will have 'edge effects' at the beginning and end of sequences, e.g. TTT 3' will match but that isn't TTTN.
  • I like your count method, but I think finditer is better because it gives you location information
  • The dictionary is sort of fine. Personally, I would dump the dictionary into a pandas dataframe This way the results are output as e.g. CSV, spreadsheet flat-files for your colleagues to view.
  • the if and elif would be placed in a subroutine, this is good practice to convert to an object orientated approach.
$\endgroup$
2
  • $\begingroup$ Great, thank you @Michael! This already gives me lots of stuff to research. $\endgroup$
    – heibai
    Nov 6 '20 at 16:07
  • $\begingroup$ Upvotes and 'accept' preferred. $\endgroup$
    – M__
    Nov 6 '20 at 16:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.