Difference between revisions of "Findseq"

From PyMOLWiki
Jump to navigation Jump to search
m (add reference to Psico function select_pepseq)
 
(7 intermediate revisions by 2 users not shown)
Line 1: Line 1:
 +
{{Infobox script-repo
 +
|type      = script
 +
|filename  = findseq.py
 +
|author    = [[User:Inchoate|Jason Vertrees]]
 +
|license  = BSD
 +
}}
 +
 
= Overview & Motivation =
 
= Overview & Motivation =
 
Anyone ever give you a protein and then say, find the sequence "FLVEW"?  Well, this script will find any string or regular expression in a given object and return that selection for you.  Here's an example,
 
Anyone ever give you a protein and then say, find the sequence "FLVEW"?  Well, this script will find any string or regular expression in a given object and return that selection for you.  Here's an example,
−
<source lang="python">
+
<syntaxhighlight lang="python">
 +
reinitialize
 +
import findseq
 
# fetch two sugar-binding PDB
 
# fetch two sugar-binding PDB
−
fetch 1tvn
+
fetch 1tvn, async=0
 
# Now, find FATEW in 1tvn, similarly
 
# Now, find FATEW in 1tvn, similarly
−
findSeq FATEW, 1tvn
+
findseq FATEW, 1tvn
 
# lower-case works, too
 
# lower-case works, too
−
findSeq fatew, 1tvn
+
findseq fatew, 1tvn
 
# how about a regular expression?
 
# how about a regular expression?
−
findSeq F.*W, 1tvn
+
findseq F.*W, 1tvn
−
 
+
 
# Find the regular expression:
 
# Find the regular expression:
 
#  ..H[TA]LVWH
 
#  ..H[TA]LVWH
 
# in the few proteins loaded.
 
# in the few proteins loaded.
 
# I then showed them as sticks and colored them to highlight matched AAs
 
# I then showed them as sticks and colored them to highlight matched AAs
−
for x in cmd.get_names(): seqFinder( "..H[TA]LVWH", x, "sele_"+x, firstOnly=1 )
+
for x in cmd.get_names(): findseq.findseq("..H[TA]LVWH", x, "sele_"+x, firstOnly=1)
−
</source>
+
</syntaxhighlight>
  
 
[[Image:SeqFinder.png|center|thumb|400px|Red residues were those matching the regular expression '..H[TA]LVWH'.]]
 
[[Image:SeqFinder.png|center|thumb|400px|Red residues were those matching the regular expression '..H[TA]LVWH'.]]
Line 22: Line 31:
 
= Usage =
 
= Usage =
 
I built this to be rather flexible.  You call it as:
 
I built this to be rather flexible.  You call it as:
−
<source lang="python">findSeq needle, haystack[, selName[, het[, firstOnly ]]]</source> where the options are:
+
<source lang="python">findseq needle, haystack[, selName[, het[, firstOnly ]]]</source> where the options are:
 
:: '''needle''' the sequence of amino acids to find.  Should be a string of one letter amino acid abbreviations.  Can also be a string-style regular expression (eg. FW.*QQ).
 
:: '''needle''' the sequence of amino acids to find.  Should be a string of one letter amino acid abbreviations.  Can also be a string-style regular expression (eg. FW.*QQ).
 
:: '''haystack''' the PyMOL object or selection in which to search
 
:: '''haystack''' the PyMOL object or selection in which to search
Line 29: Line 38:
 
:: '''firstOnly''' 0/1 -- if 0 then all matches are selected and returned; if 1 then only the first is returned
 
:: '''firstOnly''' 0/1 -- if 0 then all matches are selected and returned; if 1 then only the first is returned
  
−
= The Code =
+
= See Also =
−
<source lang="python">
+
select_pepseq@[[Psico]]
−
# -*- coding: utf-8 -*-
 
−
#
 
−
# -- seqFinder.py
 
−
#
 
−
# Finds given amino acids in a protein & returns the selection.
 
−
#
 
−
#
 
−
import re
 
−
import types
 
−
import random
 
−
 
 
−
import pymol
 
−
from pymol import stored, cmd
 
−
 
 
−
def checkParams(needle,haystack,selName,het,firstOnly):
 
−
"""
 
−
This is just a helper function for checking the user input
 
−
"""
 
−
# check Needle
 
−
if len(needle)==0 or type(needle)!=types.StringType:
 
−
print "Error: Please provide a string 'needle' to search for."
 
−
print "Error: For help type 'help motifFinder'."
 
−
return False
 
−
 
−
# check Haystack
 
−
if len(haystack)==0 or type(haystack)!=types.StringType:
 
−
print "Error: Please provide valid PyMOL object or selection name"
 
−
print "Error: in which to search."
 
−
print "Error: For help type 'help motifFinder'."
 
−
return False
 
−
 
 
−
# check het
 
−
try:
 
−
het = bool(int(het))
 
−
except ValueError:
 
−
print "Error: The 'het' parameter was not 0 or 1."
 
−
return False
 
−
 
 
−
# check first Only
 
−
try:
 
−
firstOnly = bool(int(het))
 
−
except ValueError:
 
−
print "Error: The 'firstOnly' parameter was not 0 or 1."
 
−
return False
 
−
 
 
−
# check selName
 
−
if type(selName)!=types.StringType:
 
−
print "Error: selName was not a string."
 
−
return False
 
−
 
 
−
return True
 
−
 
 
−
 
 
−
def seqFinder(needle, haystack, selName=None, het=0, firstOnly=0):
 
−
"""
 
−
DESCRIPTION:
 
−
Given a sequence/regex to find, select those
 
−
matching amino acids in the protein.
 
−
 
 
−
USAGE:
 
−
seqFinder needle, haystack[, selName[, het[, firstOnly]]]
 
−
 
 
−
PARAMS:
 
−
needle (string)
 
−
the sequence of amino acids to match and select
 
−
in the haystack.  This can be a sequence of amino
 
−
acids, or a string-style regular expression.  See
 
−
examples.
 
−
 
−
hastack (string or PyMOL selection)
 
−
name of the PyMOL object/selection in which
 
−
to find the needle.
 
−
 
 
−
selName (string; defaults to None)
 
−
This is the name of the selection to return.  If selName
 
−
is left blank (None), then the selection name will be
 
−
foundSeqXYZ where XYZ is some random number; if selName is
 
−
"sele" the usual PyMOL "(sele)" will be used; and, lastly,
 
−
if selName is anything else, that name will be used verbatim.
 
−
 
 
−
het (0 or 1; defaults to 0)
 
−
This boolean flag allows (1) or disallows (0) heteroatoms
 
−
from being considered.
 
−
 
 
−
firstOnly (0 or 1; defaults to 0)
 
−
Subsequences or motifs might be repeated, this controls how we
 
−
consider multiple matches.  If firstOnly is False (0) then we return
 
−
all found subsequences; if firstOnly is True (1), then we just return
 
−
the first found sequence.
 
−
 
 
−
RETURNS:
 
−
a newly created selection with the atoms you sought.  If there are
 
−
more than two contiguous regions, then a newly created group is
 
−
returned with each contiguous segment its own selection.
 
−
 
 
−
EXAMPLE:
 
−
# find SPVI in 1h12, foundSeqXYZ as return name
 
−
seqFinder SPVI, 1h12
 
−
 
−
# find FATEW and make it (sele).
 
−
seqFinder FATEW, 1g01, sele
 
−
 
−
# find the regular expression GMS.*QWY in 1a3h
 
−
# and put the return value in (sele).
 
−
fetch 1a3h
 
−
# this ends up finding the sequence, GMSSHGLQWY
 
−
seqFinder GMS.*QWY, 1a3h, sele
 
−
 
 
−
NOTES:
 
−
Assumes we're using the ONE LETTER amino acid abbreviations.
 
−
 
 
−
AUTHOR:
 
−
Jason Vertrees, 2009.
 
−
 
 
−
"""
 
−
# input checking
 
−
if not checkParams(needle, haystack, selName, het, firstOnly):
 
−
print "There was an error with a parameter.  Please see"
 
−
print "the above error message for how to fix it."
 
−
return None
 
−
 
 
−
one_letter ={'VAL':'V', 'ILE':'I', 'LEU':'L', 'GLU':'E', 'GLN':'Q', \
 
−
'ASP':'D', 'ASN':'N', 'HIS':'H', 'TRP':'W', 'PHE':'F', 'TYR':'Y',    \
 
−
'ARG':'R', 'LYS':'K', 'SER':'S', 'THR':'T', 'MET':'M', 'ALA':'A',    \
 
−
'GLY':'G', 'PRO':'P', 'CYS':'C'}
 
−
 
−
# remove hetero atoms (waters/ligands/etc) from consideration?
 
−
if het:
 
−
cmd.select( "__h", "br. " + haystack )
 
−
else:
 
−
cmd.select( "__h", "br. " + haystack + " and not het" )
 
−
 
−
# get the AAs in the haystack
 
−
aaDict = { 'aaList' : [] }
 
−
cmd.iterate("(name ca) and __h","aaList.append((resi,resn))",space=aaDict)
 
−
 
−
IDs = map( lambda x: int(x[0]), aaDict['aaList'] )
 
−
AAs = ''.join(map( lambda x: one_letter[x[1]], aaDict['aaList'] ))
 
−
 
 
−
reNeedle = re.compile( needle.upper() )
 
−
it = reNeedle.finditer( AAs )
 
−
 
 
−
# set the name of the selection to return.
 
−
if selName == None:
 
−
rSelName = "foundSeq" + str(random.randint(0,32000))
 
−
elif selName == "sele":
 
−
rSelName = "sele"
 
−
else:
 
−
rSelName = selName
 
−
 
 
−
# make an empty selection to which we add residues
 
−
cmd.select( rSelName, 'None')
 
−
 
 
−
for i in it:
 
−
(start,stop) = i.span()
 
−
cmd.select( rSelName, rSelName + " __h and i. " + str(IDs[start]) + "-" + str(IDs[stop-1]))
 
−
if int(firstOnly):
 
−
break
 
−
 
 
−
cmd.delete("__h")
 
−
return rSelName
 
−
 
 
−
cmd.extend("seqFinder", seqFinder )
 
−
</source>
 
  
 
[[Category:Script_Library]]
 
[[Category:Script_Library]]
 
[[Category:ObjSel_Scripts]]
 
[[Category:ObjSel_Scripts]]
 +
[[Category:Pymol-script-repo]]

Latest revision as of 11:19, 3 January 2013

Type Python Script
Download findseq.py
Author(s) Jason Vertrees
License BSD
This code has been put under version control in the project Pymol-script-repo

Overview & Motivation

Anyone ever give you a protein and then say, find the sequence "FLVEW"? Well, this script will find any string or regular expression in a given object and return that selection for you. Here's an example, 

reinitialize
import findseq
# fetch two sugar-binding PDB
fetch 1tvn, async=0
# Now, find FATEW in 1tvn, similarly
findseq FATEW, 1tvn
# lower-case works, too
findseq fatew, 1tvn
# how about a regular expression?
findseq F.*W, 1tvn
 
# Find the regular expression:
#  ..H[TA]LVWH
# in the few proteins loaded.
# I then showed them as sticks and colored them to highlight matched AAs
for x in cmd.get_names(): findseq.findseq("..H[TA]LVWH", x, "sele_"+x, firstOnly=1)
Red residues were those matching the regular expression '..H[TA]LVWH'.

Usage

I built this to be rather flexible. You call it as: 

findseq needle, haystack[, selName[, het[, firstOnly ]]]

where the options are:

needle the sequence of amino acids to find. Should be a string of one letter amino acid abbreviations. Can also be a string-style regular expression (eg. FW.*QQ).
haystack the PyMOL object or selection in which to search
selName the name of the returned selection. If you leave this blank, it'll be foundSeqXYZ where XYZ is some random integer (eg. foundSeq1435); if you supply sele then the usual PyMOL (sele) is used; and, finally, if it's anything else, then that will be used verbatim. Defaults to foundSeqXYZ so as not to overwrite any selections you might have in sele.
het 0/1 -- if 0 then heteroatoms are not considered; if 1 then they are; defaults to 0
firstOnly 0/1 -- if 0 then all matches are selected and returned; if 1 then only the first is returned

See Also

select_pepseq@Psico

Retrieved from "http://pymolwiki.org/index.php?title=Findseq&oldid=10922"