<p>The easiest solution is <code>finditer</code></p>
<pre><code>use re

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

iterator = p.finditer(useTF)
for match in iterator:
    print(match.span())
</code></pre>
<p>** Example Output**</p>
<pre><code>(0, 2)
(22, 24)
(29, 31)
</code></pre>
<p>Easy does it. Its worth looking at the documentation for <code>re</code> here, <a href="https://docs.python.org/3/howto/regex.html" rel="nofollow noreferrer">https://docs.python.org/3/howto/regex.html</a> for alternative approaches using <code>match</code> or <code>search</code> below. <code>finditer</code> is definitely the solution in context.</p>
<hr />
<p>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.</p>
<p>In summary, I kinda like your substrings and do like the count method. I'm not a fan of your <code>==</code> in this instance and your C-style matching should be substituted for reg-ex which in Python is implemented via <code>re</code>.</p>
<p>The basic issues are,</p>
<ul>
<li>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]</li>
<li>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]</li>
<li>Chopping of -2 from the length ... I dunno I'd check out -1 [minor]</li>
<li>C-style matching. I don't like this. I would definitely use <code>re</code>. One solution uses <code>match</code> or <code>search</code>, but the solution above uses <code>finditer</code>, Reg-ex will allow to search for NGG and TTTN, which you are not doing</li>
<li>I think you are looping a locus of 3 and 4 bp and I'm not sure how the C-style <code>==</code> 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.</li>
<li>I like your <code>count</code> method, but I think <code>finditer</code> is better because it gives you location information</li>
<li>The dictionary is sort of fine. Personally, I would dump the dictionary into a <code>pandas</code> dataframe This way the results are output as e.g. CSV, spreadsheet flat-files for your colleagues to view.</li>
<li>the <code>if</code> and <code>elif</code> would be placed in a subroutine, this is good practice to convert to an object orientated approach.</li>
</ul>

CC BY-SA 4.0
2020-11-06T16:07:35.090Z
Great, thank you @Michael! This already gives me lots of stuff to research.
codexero
CC BY-SA 4.0
2020-11-06T16:39:45.400Z
Upvotes and 'accept' preferred.
M__
CC BY-SA 4.0
2020-11-04T13:31:21.460Z
2020-11-06T16:37:46.250Z
2020-11-06T16:37:46.250Z
M__
<p>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:</p>
<pre><code>TargetFile = open('file.fasta')
readTF = TargetFile.read()
useTF = &quot;&quot;.join(readTF[66:].split())


#no_of_basepairs = len(useTF)


sites = {&quot;ngg_strand&quot; : [],
         &quot;ngg_antistrand&quot; : [],     #location counted from reverse direction
         &quot;tttn_strand&quot; : [],
         &quot;tttn_antistrand&quot; : []}    #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 == &quot;GG&quot;:
        sites[&quot;ngg_strand&quot;].append(i)
    elif gg_mer == &quot;CC&quot;:
        #look for reverse complement, record location on antistrand from reverse direction
        sites[&quot;ngg_antistrand&quot;].append(len(useTF)-i-2)
    elif ttt_mer == &quot;TTT&quot;:
        sites[&quot;tttn_strand&quot;].append(i+1)
    elif ttt_mer == &quot;AAA&quot;:
        # look for reverse complement, record location on antistrand from reverse direction
        sites[&quot;tttn_antistrand&quot;].append(len(useTF)-i-2)


print(len(sites[&quot;ngg_strand&quot;]))
print(useTF.count(&quot;GG&quot;))
</code></pre>
<p>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 &quot;NNNGGGGNNN&quot; my code would count 3 occurences of the &quot;GG&quot; substring, whereas the .count() method returns 2 occurences.</p>
<p>My Questions:</p>
<ol>
<li><p>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 &quot;triplet-steps&quot;. However, I am very much new to this and might be missing something.</p>
</li>
<li><p>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)?</p>
</li>
<li><p>Like I said, I am new to this so if you see any other glaring problems I would appreciate any and all help. :-)</p>
</li>
</ol>
<p>Thank you!</p>

CC BY-SA 4.0
2020-11-03T10:33:52.323Z
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?
M__
CC BY-SA 4.0
2020-11-03T10:49:01.183Z
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
codexero
CC BY-SA 4.0
2020-11-03T11:36:25.320Z
@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?
codexero
CC BY-SA 4.0
2020-11-03T15:35:45.490Z
Okay @Michael, thank you for your help.
codexero
CC BY-SA 4.0
2020-11-03T15:35:56.367Z
@ATpoint it is an informal assignment given to me by one of the postdocs at the department. I'll have a look at the specialized software out of curiosity anyway, thank you!
codexero
<crispr><python>
Find CRISPR PAM-sites with python
# Find CRISPR PAM-sites with python

(Asked by: username_0 on 2020-11-03)

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:

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.

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

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!

**Tags:** crispr, python
> **Comments:**
> **username_1**: 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?  
> **username_0**: 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  
> **username_0**: @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?  
> **username_0**: Okay @Michael, thank you for your help.  
> **username_0**: @ATpoint it is an informal assignment given to me by one of the postdocs at the department. I'll have a look at the specialized software out of curiosity anyway, thank you!  

---

## 1 answers

### Answer by username_1 on 2020-11-04. Score: 1 | Accepted answer

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.

> **Comments:**
> **username_0**: Great, thank you @Michael! This already gives me lots of stuff to research.  
> **username_1**: Upvotes and 'accept' preferred.  

---

