--- /dev/null
+###########################################################################
+# #
+# C O P Y R I G H T N O T I C E #
+# Copyright (c) 2003-10 by: #
+# * California Institute of Technology #
+# #
+# All Rights Reserved. #
+# #
+# Permission is hereby granted, free of charge, to any person #
+# obtaining a copy of this software and associated documentation files #
+# (the "Software"), to deal in the Software without restriction, #
+# including without limitation the rights to use, copy, modify, merge, #
+# publish, distribute, sublicense, and/or sell copies of the Software, #
+# and to permit persons to whom the Software is furnished to do so, #
+# subject to the following conditions: #
+# #
+# The above copyright notice and this permission notice shall be #
+# included in all copies or substantial portions of the Software. #
+# #
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, #
+# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF #
+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND #
+# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS #
+# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN #
+# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN #
+# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE #
+# SOFTWARE. #
+###########################################################################
+#
+import string
+import copy
+from math import log
+
+AA = {"TTT": "F",
+ "TTC": "F",
+ "TTA": "L",
+ "TTG": "L",
+ "TCT": "S",
+ "TCC": "S",
+ "TCA": "S",
+ "TCG": "S",
+ "TAT": "Y",
+ "TAC": "Y",
+ "TAA": "*",
+ "TAG": "*",
+ "TGT": "C",
+ "TGC": "C",
+ "TGA": "*",
+ "TGG": "W",
+ "CTT": "L",
+ "CTC": "L",
+ "CTA": "L",
+ "CTG": "L",
+ "CCT": "P",
+ "CCC": "P",
+ "CCA": "P",
+ "CCG": "P",
+ "CAT": "H",
+ "CAC": "H",
+ "CAA": "Q",
+ "CAG": "Q",
+ "CGT": "R",
+ "CGC": "R",
+ "CGA": "R",
+ "CGG": "R",
+ "ATT": "I",
+ "ATC": "I",
+ "ATA": "I",
+ "ATG": "M",
+ "ACT": "T",
+ "ACC": "T",
+ "ACA": "T",
+ "ACG": "T",
+ "AAT": "N",
+ "AAC": "N",
+ "AAA": "K",
+ "AAG": "K",
+ "AGT": "S",
+ "AGC": "S",
+ "AGA": "R",
+ "AGG": "R",
+ "GTT": "V",
+ "GTC": "V",
+ "GTA": "V",
+ "GTG": "V",
+ "GCT": "A",
+ "GCC": "A",
+ "GCA": "A",
+ "GCG": "A",
+ "GAT": "D",
+ "GAC": "D",
+ "GAA": "E",
+ "GAG": "E",
+ "GGT": "G",
+ "GGC": "G",
+ "GGA": "G",
+ "GGG": "G",
+ "---": "-",
+ "NNN": "X"
+}
+
+
+def getAA(codon):
+ """ returns one-letter AA symbol corresponding to codon, if existing.
+ returns X for unknown codon, '-' for a complete gap, and '*' for a
+ stop codon.
+ """
+ codon = codon.upper()
+ codon = string.replace(codon, "U", "T")
+ try:
+ aa = AA[codon]
+ except KeyError:
+ aa = "X"
+
+ return aa
+
+
+def translate(mRNA, frame=1):
+ """ translate a sequence into protein based on frame (1, 2, 3 only)
+ """
+ prot = ""
+ if frame < 1:
+ frame = 1
+ elif frame > 3:
+ frame = 3
+
+ for nucPos in range(frame - 1, len(mRNA) - 2, 3):
+ theCodon = mRNA[nucPos:nucPos+3]
+ print theCodon
+ theAA = getAA(theCodon)
+ prot += theAA
+
+ return prot
+
+
+def iCodon(codon):
+ """ returns the number of possible synonymous changes at
+ each site of the codon, e.g (1, 0, 3) for CTG (Leucine)
+ """
+ isynonyms = [0, 0, 0]
+ origAA = getAA(codon)
+ if origAA == "X" or origAA == "-":
+ return [-1, -1, -1]
+
+ for pos in range(len(codon)):
+ isyn = 0
+ bases = ["A", "C", "G", "T"]
+ site = codon[pos]
+ bases.remove(site)
+ for nt in bases:
+ newcodon = []
+ for nuc in codon:
+ newcodon.append(nuc)
+
+ newcodon[pos] = nt
+ newcodon = string.join(newcodon, "")
+ newAA = getAA(newcodon)
+ if newAA == origAA:
+ isyn += 1
+
+ isynonyms[pos] = isyn
+
+ return isynonyms
+
+
+def calcMutationTypes(codonList1, codonList2, verbose=False):
+ """ returns the number of (synonymous, nonsynonymous)
+ mutations in informative codons.
+ """
+ synonymous = 0.0
+ nonsynonymous = 0.0
+ for pos in range(len(codonList1)):
+ diffList = []
+ codon1 = codonList1[pos]
+ codon2 = codonList2[pos]
+ for isite in range(len(codon1)):
+ isite1 = codon1[isite]
+ isite2 = codon2[isite]
+ if isite1 != isite2:
+ diffList.append(isite)
+
+ aa1 = getAA(codon1)
+ aa2 = getAA(codon2)
+ if aa1 == aa2:
+ if aa1 == "S":
+ if codon1[0] != codon2[0]:
+ if verbose:
+ print "nonsynonymous Serines"
+
+ nonsynonymous += 1
+
+ synonymous += 1
+ else:
+ synonymous += len(diffList)
+
+ else:
+ if len(diffList) == 1:
+ nonsynonymous += 1
+ else:
+ if verbose:
+ print "parsimonious mutation path estimator - assuming 1 parameter"
+ print "diffList = %s\t%s (%s) \t%s (%s)" % (diffList, codon1, aa1, codon2, aa2)
+
+ if 1 in diffList:
+ if verbose:
+ print "middle site is nonsynonymous"
+
+ nonsynonymous += 1
+ diffList.remove(1)
+
+ if 2 in diffList:
+ codon3 = codon1[:-1] + codon1[2]
+ aa3 = getAA(codon3)
+ if aa3 == aa1:
+ if verbose:
+ print "last site is synonymous"
+
+ synonymous += 1
+ diffList.remove(2)
+
+ if 0 in diffList:
+ codon3 = codon2[0] + codon1[1:]
+ aa3 = getAA(codon3)
+ if aa3 == aa1:
+ if verbose:
+ print "first site is synonymous"
+
+ synonymous += 1
+ diffList.remove(0)
+
+ if len(diffList) > 0:
+ if verbose:
+ print "%s must be non-synonymous" % str(diffList)
+
+ nonsynonymous += len(diffList)
+
+ return (synonymous, nonsynonymous)
+
+
+def calcSubstitutionSites(codonList):
+ """ returns the number of (synonymous, nonsynonymous) sites
+ in a list of codons, which should be filtered for gaps and X's.
+ """
+ synonymous = 0.0
+ nonsynonymous = 0.0
+ for codon in codonList:
+ (site0, site1, site2) = iCodon(codon)
+ if site0 < 0:
+ continue
+
+ synonymous += (site0 + site1 + site2) / 3.0
+ nonsynonymous += (9.0 - site0 - site1 - site2) / 3.0
+
+ return (synonymous, nonsynonymous)
+
+
+def calcSubstitutionsPerSite(codonList1, codonList2):
+ """ returns the number of substitutions per site as
+ a triplet in comparable or informative sites.
+ """
+ site0 = 0.0
+ site1 = 0.0
+ site2 = 0.0
+ for pos in range(len(codonList1)):
+ codon1 = codonList1[pos]
+ codon2 = codonList2[pos]
+ if codon1[0] != codon2[0]:
+ site0 +=1
+
+ if codon1[1] != codon2[1]:
+ site1 += 1
+
+ if codon1[2] != codon2[2]:
+ site2 += 1
+
+ return (site0, site1, site2)
+
+
+def calcKs(Ms, Ns):
+ """ returns a Ks calculated using Ms and Ns and adjusted using
+ Jukes and Cantor's formula.
+ """
+ Ks = -0.75 * log(1 - ((4.0/3.0) * Ms / Ns))
+ return Ks
+
+
+def calcKa(Ma, Na):
+ """ returns a Ka calculated using Ma and Na and adjusted using
+ Jukes and Cantor's formula.
+ """
+ Ka = -0.75 * log(1 - ((4.0/3.0) * Ma / Na))
+ return Ka
+
+
+def printCDSdict(cdsDict, printAA=True):
+ """ Prints every locus in a cdsDict. Optionally truns off AA translation.
+ """
+ for locus in cdsDict:
+ printCDSlocus(cdsDict, locus, printAA)
+
+
+def printCDSlocus(cdsDict, locus, printAA=True):
+ """ Prints a locus in a given cdsDict. Optionally turns off AA translation.
+ """
+ cdsOutLines = []
+ aaOutLines = []
+ cdsOutLine = locus + "\t"
+ aaOutLine = " " * len(locus) + "\t"
+ for pos in range(len(cdsDict[locus])):
+ if len(cdsOutLine) > 69:
+ cdsOutLines.append(cdsOutLine + cdsDict[locus][pos] + "\n")
+ cdsOutLine = locus + "\t"
+ aaOutLines.append(aaOutLine + " " + getAA(cdsDict[locus][pos]) + " \n")
+ aaOutLine = " " * len(locus) + "\t"
+ else:
+ cdsOutLine += cdsDict[locus][pos] + " "
+ aaOutLine += " " + getAA(cdsDict[locus][pos]) + " "
+
+ cdsOutLines.append(cdsOutLine + "\n")
+ aaOutLines.append(aaOutLine + " \n")
+ for index in range(len(cdsOutLines)):
+ print cdsOutLines[index]
+ if printAA:
+ print aaOutLines[index]
+
+
+def getComparableSites(cdsDict, loci=[]):
+ """ given a cdsDict and a list of loci, returns a new cdsDict with positions with
+ gaps or X codons in any one sequence deleted in all sequences.
+ """
+ newDict = {}
+ seqArray = []
+ locArray = []
+ deleteList = []
+ index = 0
+ for locus in loci:
+ locArray.append(locus)
+ seqArray.append(copy.deepcopy(cdsDict[locus]))
+
+ for index in range(len(locArray)):
+ for pos in range(len(seqArray[index])):
+ aa = getAA(seqArray[index][pos])
+ if aa == "X" or aa == "-":
+ if pos not in deleteList:
+ deleteList.append(pos)
+
+ deleteList.sort()
+ deleteList.reverse()
+ for pos in deleteList:
+ for index in range(len(locArray)):
+ del seqArray[index][pos]
+
+ for index in range(len(locArray)):
+ newDict[locArray[index]] = seqArray[index]
+
+ return newDict
+
+
+def getInformativeSites(cdsDict, loci=[]):
+ """ given a cdsDict of comparable sites and a list of loci, returns a new
+ cdsDict with positions with only codons that differ in one or more
+ sequences and which are therefore (possibly) informative.
+ """
+ newDict = {}
+ seqArray = []
+ locArray = []
+ deleteList = []
+ index = 0
+ for locus in loci:
+ locArray.append(locus)
+ seqArray.append(copy.deepcopy(cdsDict[locus]))
+
+ for pos in range(len(seqArray[0])):
+ deleteCodon = True
+ refcodon = seqArray[0][pos]
+ for index in range(1, len(locArray)):
+ seqcodon = seqArray[index][pos]
+ if seqcodon != refcodon:
+ deleteCodon = False
+
+ if deleteCodon:
+ deleteList.append(pos)
+
+ deleteList.sort()
+ deleteList.reverse()
+ for pos in deleteList:
+ for index in range(len(locArray)):
+ del seqArray[index][pos]
+
+ for index in range(len(locArray)):
+ newDict[locArray[index]] = seqArray[index]
+
+ return newDict
+
+
+def buildCDSdict(cdsFileName):
+ """ imports a set of *ALIGNED* sequences in a fasta-format file and splits
+ each sequence into its individual codons. Returns a Dictionary of the
+ sequences with the fasta ID as the key and the codons in a list.
+ """
+ cdsfile = open(cdsFileName, "r")
+ cdslines = cdsfile.readlines()
+ cdsfile.close()
+ cdsDict = {}
+ locus = ""
+ for line in cdslines:
+ partialCodon = ""
+ inFrame = True
+ line = line[:-1]
+ if line[0] == ">":
+ fields = line.split(" ")
+ if len(fields[0]) > 1 and fields[0][0] == ">":
+ locus = fields[0][1:]
+ else:
+ locus = fields[1]
+
+ cdsDict[locus] = []
+ else:
+ for pos in range(0, len(line), 3):
+ codon = line[pos:pos+3]
+ if codon == "---":
+ cdsDict[locus].append("---")
+ elif "-" in codon and inFrame:
+ inFrame = False
+ partialCodon = string.replace(codon, "-", "")
+ cdsDict[locus].append("---")
+ elif not inFrame:
+ partialCodon += string.replace(codon, "-", "")
+ if len(partialCodon) == 3:
+ cdsDict[locus].append(partialCodon)
+ inFrame = True
+ partialCodon = ""
+
+ if len(partialCodon) > 3:
+ cdsDict[locus].append(partialCodon[:3])
+ partialCodon = partialCodon[3:]
+ else:
+ cdsDict[locus].append(codon)
+
+ return cdsDict
\ No newline at end of file
projects / erange.git / blobdiff