scripts for supplementary
authorFabio Zanini <fabio.zanini@tuebingen.mpg.de>
Mon, 4 Feb 2013 14:28:07 +0000 (15:28 +0100)
committerFabio Zanini <fabio.zanini@tuebingen.mpg.de>
Mon, 4 Feb 2013 14:28:07 +0000 (15:28 +0100)
scripts/codon_entropy_synonymous_subtypeB.py [new file with mode: 0644]
scripts/conservation_syn_nonsyn_subtypeB.py [new file with mode: 0644]
scripts/modules/__init__.py
scripts/modules/__init__.pyc
scripts/modules/alignment.py [new file with mode: 0644]
scripts/modules/alignment.pyc [new file with mode: 0644]
scripts/modules/codon_usage.py [new file with mode: 0644]
scripts/modules/codon_usage.pyc [new file with mode: 0644]
scripts/modules/parser_reference.py [new file with mode: 0644]
scripts/modules/parser_reference.pyc [new file with mode: 0644]

diff --git a/scripts/codon_entropy_synonymous_subtypeB.py b/scripts/codon_entropy_synonymous_subtypeB.py
new file mode 100644 (file)
index 0000000..db5548f
--- /dev/null
@@ -0,0 +1,139 @@
+# vim: fdm=indent
+'''
+author:     Fabio Zanini
+date:       07/12/12
+content:    Measure the level of conservation of synonymous mutations in subtype
+            B.
+'''
+# Modules
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.patches import Rectangle
+
+# Custom modules
+import sys
+sys.path.insert(0, '.')
+import modules.alignment as alignment
+from modules.alphabet import alpha
+from modules.codon_usage import number_of_codons
+from modules.helper import translate
+from modules.parser_reference import seq_patterns, find_with_gaps
+
+
+
+# Globals
+genes = ['gag', 'pol', 'env', 'nef']
+reference = 'SHAPE'
+seq_patterns = seq_patterns[reference]
+
+
+
+# Functions
+def smoothen(x, f):
+    l = len(x)
+    x = x[:l - l % f]
+    x = x.reshape((len(x) / f, f)).mean(axis=1)
+    return x
+
+
+
+# Script
+if __name__ == '__main__':
+
+    # Load data for a certain gene
+    gene = genes[0]
+    ali = alignment.MSA(gene=gene, reference=reference)
+    msa = np.array(ali.MSA)
+    afs = ali.allele_frequencies()
+
+    # Filter only part of pol
+    if gene == 'pol':
+        msa = msa[:, :2718]
+        afs = afs[:, :2718]
+
+    # Find consensus sequence
+    consind = afs.argmax(axis=0)
+    consensus = alpha[consind]
+
+    # Some seqs are not viable and include frameshifts that mess up the
+    # translation, hence we restrict to positions for which gaps are a minority
+    is_gap = consensus == '-'
+    # Exclude full codons
+    tmp = np.unique(is_gap.nonzero()[0] / 3)
+    is_gap[np.concatenate([tmp * 3, tmp * 3 +1, tmp * 3 + 2])] = True
+
+    # Exclude stop codons
+    is_stop = np.zeros_like(is_gap)
+    tmp = (translate(consensus) == '*').nonzero()[0]
+    is_stop[np.concatenate([tmp * 3, tmp * 3 +1, tmp * 3 + 2])] = True
+
+    ## Plot base prevalence
+    #for i in xrange(4):
+    #    plt.plot(np.arange(len(consensus)), afs[i],
+    #             lw=1.5, alpha=0.5)
+    #plt.xlim(0, 2600)
+    #plt.ylim(-0.05, 1.25)
+    #plt.xlabel('position in '+gene)
+    #plt.ylabel('allele frequency')
+    #plt.legend(alpha, loc=9)
+
+    # Good are codons with no gaps and no stops
+    is_good = (-is_gap) & (-is_stop)
+
+    # For each codon, calculate the entropy
+    msa = msa[:, is_good]
+    consaa = translate(consensus[is_good])
+    entropy = np.zeros(len(consaa))
+    from collections import Counter
+    for i, aa in enumerate(consaa):
+        tmp = msa[:, i * 3: (i + 1) * 3]
+        count = Counter(map(''.join, tmp))
+        abus = []
+        for (cod, abu) in count.iteritems():
+            if translate(np.array(list(cod))) == aa:
+                abus.append(abu)
+        abus = np.array(abus)
+        freqcod = 1.0 * abus / abus.sum()
+        entropy[i] = -np.sum(freqcod * np.log(freqcod))
+
+    # Maximal entropy per codon
+    ncod = np.array(map(number_of_codons, consaa))
+    entropymax = np.log(ncod)
+    entropyfrac = entropy / (1e-5 + entropymax)
+
+    # Plot the entropy
+    plt.figure()
+    plt.plot(np.arange(len(entropy)), entropyfrac, c='b', lw=1)
+    plt.xlabel('Position in '+gene)
+    plt.ylabel('Codon entropy / maximal entropy')
+    plt.ylim(0, 1)
+
+    # Plot smoothed curve
+    smoothrange = 15
+    xsmooth = smoothen(np.arange(len(entropy)), smoothrange)
+    ysmooth = smoothen(entropyfrac, smoothrange)
+    plt.plot(xsmooth, ysmooth, c='r', lw=2, ls='-')
+    plt.title(gene, fontsize=18)
+
+    # Mark V loops
+    if gene == 'env':
+        refstr = str(ali.reference_seq.seq)
+        tmp = (-is_good).nonzero()[0]
+        Vind = {}
+        for key in ['V1', 'V2', 'V3', 'V4', 'V5']:
+            i0 = find_with_gaps(refstr, seq_patterns[key][0])
+            i0 -= (tmp < i0).sum()
+            i0 //= 3
+            i1 = find_with_gaps(refstr, seq_patterns[key][1])
+            i1 -= (tmp < i1).sum()
+            i1 //= 3
+            Vind[key] = (i0, i1)
+            h = Rectangle((i0, 0), i1 - i0, 1, color='yellow', alpha=0.5)
+            plt.gca().add_patch(h)
+            plt.text(i0 + 0.2 * (i1 - i0), 0.95, key, fontsize=12)
+        
+
+    plt.ion()
+    plt.show()
+
+    
diff --git a/scripts/conservation_syn_nonsyn_subtypeB.py b/scripts/conservation_syn_nonsyn_subtypeB.py
new file mode 100644 (file)
index 0000000..1e37fb4
--- /dev/null
@@ -0,0 +1,168 @@
+# vim: fdm=indent
+'''
+author:     Fabio Zanini
+date:       07/12/12
+content:    Measure the level of conservation of synonymous and nonsynonymous 
+            mutations in subtype B. Nonsynonymous mutations are expected to be
+            more broadly conserved, because they impair protein function.
+'''
+# Modules
+import sys
+from collections import Counter
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.patches import Rectangle
+
+# Custom modules
+import modules.alignment as alignment
+from modules.alphabet import alpha, alphal_ng
+from modules.codon_usage import all_codons
+from modules.codon_usage import translate as transstrict
+from modules.helper import translate
+from modules.parser_reference import seq_patterns, find_with_gaps
+
+
+
+# Globals
+gene = 'gag'
+reference = 'SHAPE'
+seq_patterns = seq_patterns[reference]
+
+
+
+# Functions
+def smooth(x, f):
+    l = len(x)
+    x = x[:l - l % f]
+    x = x.reshape((len(x) / f, f)).mean(axis=1)
+    return x
+
+
+def codon_single_mutants_synnonsyn(codon):
+    codon = list(codon)
+    aa = transstrict(''.join(codon))
+    syn = []
+    nonsyn = []
+    for pos in xrange(3):
+        for mut in alphal_ng:
+            if mut != codon[pos]:
+                tmp = list(codon)
+                tmp[pos] = mut
+                tmp = ''.join(tmp)
+                if transstrict(tmp) == aa:
+                    syn.append(tmp)
+                else:
+                    nonsyn.append(tmp)
+    return {'syn': syn, 'nonsyn': nonsyn}
+
+
+
+
+# Script
+if __name__ == '__main__':
+
+    # Load data for a certain gene
+    ali = alignment.MSA(gene=gene, reference=reference)
+    msaraw = np.array(ali.MSA)
+    afs = ali.allele_frequencies()
+
+    # Filter only part of pol
+    if gene == 'pol':
+        msaraw = msaraw[:, :2718]
+        afs = afs[:, :2718]
+
+    # Find consensus sequence
+    consind = afs.argmax(axis=0)
+    consensus = alpha[consind]
+
+    # Some seqs are not viable and include frameshifts that mess up the
+    # translation, hence we restrict to positions for which gaps are a minority
+    is_gap = consensus == '-'
+    # Exclude full codons
+    tmp = np.unique(is_gap.nonzero()[0] / 3)
+    is_gap[np.concatenate([tmp * 3, tmp * 3 +1, tmp * 3 + 2])] = True
+
+    # Exclude stop codons
+    is_stop = np.zeros_like(is_gap)
+    tmp = (translate(consensus) == '*').nonzero()[0]
+    is_stop[np.concatenate([tmp * 3, tmp * 3 +1, tmp * 3 + 2])] = True
+
+    ## Plot base prevalence
+    #for i in xrange(4):
+    #    plt.plot(np.arange(len(consensus)), afs[i],
+    #             lw=1.5, alpha=0.5)
+    #plt.xlim(0, 2600)
+    #plt.ylim(-0.05, 1.25)
+    #plt.xlabel('position in '+gene)
+    #plt.ylabel('allele frequency')
+    #plt.legend(alpha, loc=9)
+
+    # Good are codons with no gaps and no stops
+    is_good = (-is_gap) & (-is_stop)
+
+    # Get all syn and nonsyn single mutants of all codons
+    sms = {cod: codon_single_mutants_synnonsyn(cod) for cod in all_codons()}
+
+    # For each codon, calculate the number of syn and nonsyn changes from the
+    # consensus
+    msa = msaraw[:, is_good]
+    conscods = consensus[is_good].reshape((is_good.sum() / 3, 3))
+    # Count the number of sequences that enter the statistics for each codon, to
+    # counterbalance for ambiguous sequences
+    nseqs = msa.shape[0] * np.ones(len(conscods))
+    # Prepare the structure to save the counts (columns: not changed, syn, nonsyn)
+    counts3 = np.zeros((len(conscods), 3), int)
+    
+    # Run over consensus codons and count
+    for i, conscod in enumerate(conscods):
+        conscod = ''.join(conscod)
+        smscod = sms[conscod]
+        codsall = map(''.join, msa[:, i * 3: (i + 1) * 3])
+        counts = Counter(codsall)
+        for (cod, abundance) in counts.iteritems():
+            if cod == conscod:
+                counts3[i, 0] += abundance
+            elif cod in smscod['syn']:
+                counts3[i, 1] += abundance
+            elif cod in smscod['nonsyn']:
+                counts3[i, 2] += abundance
+
+    # Calculate the number of syn and nonsyn for each codon along the consensus
+    nsmscons = {'syn': [],
+                'nonsyn': []}
+    for cod in conscods:
+        tmp = sms[''.join(cod)]
+        for key in ['syn', 'nonsyn']:
+            nsmscons[key].append(len(tmp[key]))
+    for key in ['syn', 'nonsyn']:
+        nsmscons[key] = np.array(nsmscons[key])
+
+    # Take e.g. four-fold degenerate codons (as SMs only) and calculate frequencies
+    is_deg = nsmscons['syn'] == 3
+    nsmscons2 = np.vstack([nsmscons['syn'], nsmscons['nonsyn']]).T
+    freqs2 = 1.0 * counts3[:, 1:][is_deg] / nsmscons2[is_deg] / msa.shape[0]
+
+    ## Plot a histogram of conservation levels
+    #plt.figure()
+    #plt.hist(freqs2,color=('r', 'b'), bins=np.linspace(0, 0.04, 20),
+    #         normed=True, label=('syn', 'nonsyn'))
+    #plt.xlabel('diversity')
+    #plt.ylabel('histogram')
+    #plt.title('diversity at 4-fold degen codons, '+gene)
+
+    # Plot cumulative
+    freqsyn = np.sort(freqs2[:,0])
+    freqnonsyn = np.sort(freqs2[:,1])
+    fig = plt.figure(figsize=[6,4.6])
+    plt.plot(freqsyn, np.linspace(0, 1, len(freqsyn)), c='r', lw=2, label='syn')
+    plt.plot(freqnonsyn, np.linspace(0, 1, len(freqnonsyn)), c='b', lw=2, label='nonsyn')
+    plt.xlabel('diversity')
+    plt.ylabel('cumulative')
+    plt.xlim(0, 0.04)
+    plt.legend(loc=4)
+    plt.title('diversity at 4-fold degen codons, '+gene)
+    line = plt.plot([0.003] * 2, [0.08, 0.60], c='k', lw=1.5, ls='--')
+    txt = plt.text(0.0035, 0.5, r'$\Delta \approx 0.52$', fontsize=14)
+
+    plt.ion()
+    plt.show()
index e583ee6..b37840c 100644 (file)
@@ -1,7 +1,9 @@
-__all__ = [
-           'helper',
+__all__ = ['alignment',
            'alphabet',
+           'codon_usage',
+           'helper',
            'secondary',
+           'parser_reference',
            'parser_Shankarappa',
            'parser_Bunnik2008',
            'parser_LANL',
index 61fe17a..1db5d44 100644 (file)
Binary files a/scripts/modules/__init__.pyc and b/scripts/modules/__init__.pyc differ
diff --git a/scripts/modules/alignment.py b/scripts/modules/alignment.py
new file mode 100644 (file)
index 0000000..f70660f
--- /dev/null
@@ -0,0 +1,82 @@
+# vim: fdm=indent
+'''
+author:     Fabio Zanini
+date:       12/11/12
+content:    Tools for managing sequences and alignments of general nature.
+'''
+# Modules
+import os
+import glob
+import numpy as np
+from Bio import SeqIO, AlignIO
+import general.classes.alphabet as abc
+
+
+# Globals
+datadir = '/home/fabio/university/phd/general/data/'
+genes = ['gag', 'pol', 'env', 'nef', 'vpu', 'vpr', 'vif', 'tat']
+
+
+
+# Classes
+class MSA(object):
+    '''Multiple sequence alignment with general properties'''
+    def __repr__(self):
+        s = 'MSA(kind="'+self.kind+'", gene="'+self.gene+'"'
+        if hasattr(self, 'reference'):
+            s += ', reference="'+self.reference+'"'
+        s += ')'
+        return s
+
+    def __init__(self, kind='subtypeB', reference=None, gene=None):
+
+        kinddir = (datadir + kind).rstrip('/') + '/'
+        self.kind = kind.rstrip('/')
+        if reference is not None:
+            kinddir = (kinddir + 'to_reference/')
+            self.reference = reference
+        if not os.path.isdir(kinddir):
+            raise IOError('Directory not found')
+
+        found = False
+        for fn in glob.glob(kinddir+'*.fasta'):
+            if gene in fn:
+                MSA = AlignIO.read(fn, 'fasta')
+                if reference is None:
+                    self.MSA = MSA
+                else:
+                    self.reference_seq = MSA[0]
+                    if reference not in self.reference_seq.id:
+                        raise ValueError('Reference name not found in first sequence of the MSA')
+                    self.MSA = MSA[1:]
+                self.gene = gene
+                found = True
+                break
+        if not found:
+            raise IOError('Alignment file not found')
+
+
+    def allele_frequencies(self, alpha=abc.alpha):
+        '''Get the allele frequencies for an alphabet
+        
+        Parameters:
+            - alpha: an alphabet
+        '''
+        alpha = list(alpha)
+        L = len(self.MSA[0])
+        afs = np.zeros((len(alpha), L))
+        mat = np.array(self.MSA)
+        for i, a in enumerate(alpha):
+            afs[i] = (mat == a).sum(axis=0)
+        afs /= afs.sum(axis=0)
+
+        return afs
+
+
+
+
+# Script
+if __name__ == '__main__':
+
+    msa = MSA('subtypeB', gene='nef')
+
diff --git a/scripts/modules/alignment.pyc b/scripts/modules/alignment.pyc
new file mode 100644 (file)
index 0000000..5162f3a
Binary files /dev/null and b/scripts/modules/alignment.pyc differ
diff --git a/scripts/modules/codon_usage.py b/scripts/modules/codon_usage.py
new file mode 100644 (file)
index 0000000..06c10fd
--- /dev/null
@@ -0,0 +1,143 @@
+# vim: fdm=indent
+'''
+author:     Fabio Zanini
+date:       16/10/12
+content:    Read codon usage tables
+'''
+# Modules
+import re
+from longitudinal.classes.parser_generic import alpha
+from Bio.Seq import translate
+
+
+
+# Globals
+datadir = '/home/fabio/university/phd/general/data/codon_usage/'
+alpha = list(alpha)
+
+
+
+# Classes
+class codon_table(object):
+
+    def __init__(self, organism='homo sapiens'):
+
+        self.organism = organism
+        self.datafile = datadir + re.sub(' ', '_', organism) +'.dat'
+        
+        table_within = {}
+        table_abs = {}
+        with open(self.datafile, 'r') as f:
+            for line in f:
+                if line[0] != '#':
+                    fields = line.rstrip('\n').split('\t')
+                    key = fields[0]
+                    aa = fields[1]
+                    freqwithin = float(fields[2])
+                    freqabs = 1e-3 * float(fields[3])
+
+                    table_within[key] = (aa, freqwithin)
+                    table_abs[key] = freqabs
+
+        self.within = table_within
+        self.absolute = table_abs
+
+
+
+# Functions
+def all_codons():
+    '''Generate a list with all codons'''
+    return ('AAA', 'AAC', 'AAG', 'AAT', 'ACA', 'ACC',
+            'ACG', 'ACT', 'AGA', 'AGC', 'AGG', 'AGT',
+            'ATA', 'ATC', 'ATG', 'ATT', 'CAA', 'CAC',
+            'CAG', 'CAT', 'CCA', 'CCC', 'CCG', 'CCT',
+            'CGA', 'CGC', 'CGG', 'CGT', 'CTA', 'CTC',
+            'CTG', 'CTT', 'GAA', 'GAC', 'GAG', 'GAT',
+            'GCA', 'GCC', 'GCG', 'GCT', 'GGA', 'GGC',
+            'GGG', 'GGT', 'GTA', 'GTC', 'GTG', 'GTT',
+            'TAA', 'TAC', 'TAG', 'TAT', 'TCA', 'TCC',
+            'TCG', 'TCT', 'TGA', 'TGC', 'TGG', 'TGT',
+            'TTA', 'TTC', 'TTG', 'TTT')
+
+
+def all_amino_acids():
+    '''All amino acids'''
+    return ('A', 'C', 'E', 'D', 'G', 'F',
+            'I', 'H', 'K', '*', 'M', 'L',
+            'N', 'Q', 'P', 'S', 'R', 'T',
+            'W', 'V', 'Y')
+
+
+def number_of_codons(amino_acid):
+    '''Get the number of codons per amino acid'''
+    d = {'L': 6, 'S': 6, 'R': 6, 'A': 4,
+         'G': 4, 'P': 4, 'T': 4, 'V': 4,
+         'I': 3, '*': 3, 'C': 2, 'E': 2,
+         'D': 2, 'F': 2, 'H': 2, 'K': 2,
+         'N': 2, 'Q': 2, 'Y': 2, 'M': 1,
+         'W': 1}
+    return d[amino_acid]
+
+
+def volatility(codon):
+    '''Number of nonsilent single mutants'''
+    codon = ''.join(codon)
+    if codon not in all_codons():
+        raise ValueError('Codon not recognized. Remember we are in a DNA world.')
+
+    aa = translate(codon)
+    v = 0
+    for pos in xrange(3):
+        for m in alpha:
+            if m not in ('-', codon[pos]):
+                tmp = list(codon)
+                tmp[pos] = m
+                tmp = ''.join(tmp)
+                v += translate(tmp) != aa
+    return v
+
+
+def effective_number_of_codons(codonlist):
+    '''Calculate the ENC from a list of codons'''
+    def F(n, codonlist):
+        '''The probaibility that two codons are the same in the list'''
+        import numpy as np
+        from collections import Counter
+
+        Ps = []
+        for aa in all_amino_acids():
+            if number_of_codons(aa) == n:
+                cods = filter(lambda x: translate(x) == aa, codonlist)
+                if cods:
+                    P = 1.0 * np.array(Counter(cods).values()) / len(cods)
+                    Ps.append(np.dot(P, P))
+        if Ps:
+            return np.mean(Ps)
+        else:
+            return None
+
+    codonlist = map(''.join, codonlist)
+    F2 = F(2, codonlist)
+    F3 = F(3, codonlist)
+    F4 = F(4, codonlist)
+    F6 = F(6, codonlist)
+    ENC = 2
+    if F2 is not None:
+        ENC += 9 / F2
+    if F3 is not None:
+        ENC += 1 / F3
+    if F4 is not None:
+        ENC += 5 / F4
+    if F6 is not None:
+        ENC += 3 / F6
+
+    return ENC
+
+
+                
+
+
+# Script
+if __name__ == '__main__':
+
+    table = codon_table()
diff --git a/scripts/modules/codon_usage.pyc b/scripts/modules/codon_usage.pyc
new file mode 100644 (file)
index 0000000..4e70adc
Binary files /dev/null and b/scripts/modules/codon_usage.pyc differ
diff --git a/scripts/modules/parser_reference.py b/scripts/modules/parser_reference.py
new file mode 100644 (file)
index 0000000..4019aaf
--- /dev/null
@@ -0,0 +1,109 @@
+# vim: fdm=indent
+'''
+author:     Fabio Zanini
+date:       04/09/12
+content:    Implement some globals used by many parsers that interact with a
+            reference.
+'''
+# Modules
+import numpy as np
+import Bio.AlignIO as AlignIO
+from itertools import permutations
+
+from longitudinal.classes.helper import find_with_gaps
+from longitudinal.classes.parser_generic import alpha
+
+
+
+# Globals
+# these are derived from each other using MSA.corresponding_pattern
+# (and ultimately by hand from the horrible XLS file for HXB2)
+seq_patterns = {'SHAPE':{'V1': ['TGCACTGATTTGAA', 'TGCTCTTTCAATATCAGC'],
+                         'V2': ['TGCTCTTTCAATATCAGC', 'TGTAACACCTCAGTCAT'],
+                         'V3': ['TGTACAAGACCCAACAACA', 'TGTAACATTAGTAGA'],
+                         'V4': ['TGTAATTCAACACAACTGTTTAA', 'TGCAGAATAAAACAATT'],
+                         'V5': ['TAATAACAACAATGGGTCCGAGAT', 'CCTGGAGGA'],
+                         'env_signal': ['ATGAGAGTGAAGGAGAA', 'ACAGAAAAATTGTGGGTC']},
+                'NL4-3':{'V1': ['TGCACTGATTTGAA', 'TGCTCTTTCAATATCAGC'],
+                         'V2': ['TGCTCTTTCAATATCAGC', 'TGTAACACCTCAGTCAT'],
+                         'V3': ['TGTACAAGACCCAACAACA', 'TGTAACATTAGTAGA'],
+                         'V4': ['TGTAATTCAACACAACTGTTTAA', 'TGCAGAATAAAACAATT'],
+                         'V5': ['TAATAACAACAATGGGTCCGAGAT', 'CCTGGAGGA']},
+                'HXB2': {'V1': ['TGCACTGATTTGAA', 'TGCTCTTTCAATATCAGC'],
+                         'V2': ['TGCTCTTTCAATATCAGC', 'TGTAACACCTCAGTCAT'],
+                         'V3': ['TGTACAAGACCCAACAACA', 'TGTAACATTAGTAGA'],
+                         'V4': ['TGTAATTCAACACAACTGTTTAA', 'TGCAGAATAAAACAAAT'],
+                         'V5': ['TAATAGCAACAATGAGTCCGAGAT', 'CCTGGAGGA'],
+                         'env_signal': ['ATGAGAGTGAAGGAGAA', 'ACAGAAAAATTGTGGGTC']}
+               }
+
+genes = {'SHAPE': {'gag': ['ATGGGTGCGAGAGCGTCGGTA', 'AGATAGGGGGGCAATTAA'],
+                   'pol': ['CCTCAGATCACTCTTTGGCAGCGACCCC', 'CACATGGAAAAGATTAGTA'],
+                   'env': ['ATGAGAGTGAAGGAG', 'GATGGGTGGCAAGTG'],
+                   'nef': ['ATGGGTGGCAAGTG', 'CATCGAGCTTGCTACAAGGG']},
+         'HXB2': {'gag': ['ATGGGTGCGAGAGCGTCAGTA', 'AGATAGGGGGGCAACTAA'],
+                  'pol': ['CCTCAGGTCACTCTTTGGCAACGACCCC', 'AACATGGAAAAGTTTAGTA'],
+                  'env': ['ATGAGAGTGAAGGAG', 'GATGGGTGGCAAGTG'],
+                  'nef': ['ATGGGTGGCAAGTG', 'CATCGAGCTTGCTACAAG']},
+        }
+
+datadir = '/home/fabio/university/phd/longitudinal/data/reference'
+
+
+
+class MSA(object):
+    '''An MSA of the given references.
+
+    Note: the MSA must already exist, we are not calling MUSCLE here for reasons
+    of clarity (although we could).
+    '''
+
+    def __init__(self, references):
+        '''Initialize'''
+        for tmplist in permutations(references):
+    
+            datafile = datadir+'/'+'_'.join(tmplist)
+            if len(references) > 1:
+                datafile = datafile+'_aligned'
+            datafile = datafile+'.fasta'
+    
+            try:
+                self.MSA = AlignIO.read(datafile, format='fasta')
+                break
+            except IOError:
+                pass
+    
+        if not hasattr(self, 'MSA'):
+            raise IOError('Reference alignment not found')
+
+
+    def __repr__(self):
+        return self.MSA.__repr__()
+
+
+    def __str__(self):
+        return self.MSA.__str__()
+
+
+    def corresponding_pattern(self, pattern, origin, target):
+        L = len(pattern)
+        fo = lambda x: origin in self.MSA[x].name
+        ft = lambda x: target in self.MSA[x].name
+        iseqo = filter(fo, xrange(len(self.MSA)))[0]
+        iseqt = filter(ft, xrange(len(self.MSA)))[0]
+
+        index = find_with_gaps(str(self.MSA[iseqo].seq), pattern)
+        if index == (-1):
+            raise ValueError('Pattern not found in origin')
+
+        tmplist = []
+        while (len(tmplist) < L) and (index < self.MSA.get_alignment_length()):
+            base = self.MSA[iseqt, index]
+            if base != '-':
+                tmplist.append(base)
+            index += 1
+
+        return ''.join(tmplist)
+
+
+    
diff --git a/scripts/modules/parser_reference.pyc b/scripts/modules/parser_reference.pyc
new file mode 100644 (file)
index 0000000..c6c4418
Binary files /dev/null and b/scripts/modules/parser_reference.pyc differ