1.1 --- a/allpy/homology.py	Wed Jul 06 22:28:11 2011 +0400
     1.2 +++ b/allpy/homology.py	Fri Jul 08 16:06:59 2011 +0400
     1.3 @@ -10,7 +10,7 @@
     1.4  
     1.5      DATA:
     1.6          monomer_ids = {monomer_id:class_id}
     1.7 -        classes     = {class_id:monomer_id)
     1.8 +        classes     = {class_id:[monomer_id1, monomer_id2]}
     1.9              monomer_id = (sequence_name, monomer_number) 
    1.10              class_id = string 
    1.11          columns = {class_id:column_number} optional
    1.12 @@ -19,7 +19,7 @@
    1.13  
    1.14      METHODS:
    1.15         + .read(file_name) 
    1.16 -       + .wright_monomer(file,monomer_id,monomer_homology_class_id)
    1.17 +       + .write_monomer(file,monomer_id,monomer_homology_class_id)
    1.18         + .write_class(file,list_of_monomer_ids,monomer_homology_class_id)
    1.19         + .write_block(file, block)
    1.20         + .compare_with(class)

     2.1 --- a/sequence_based_blocks_search/blocks_finder.py	Wed Jul 06 22:28:11 2011 +0400
     2.2 +++ b/sequence_based_blocks_search/blocks_finder.py	Fri Jul 08 16:06:59 2011 +0400
     2.3 @@ -3,6 +3,7 @@
     2.4  import sys
     2.5  sys.path.append("../")
     2.6  from allpy import protein
     2.7 +from allpy markups import *
     2.8  import copy
     2.9  import math
    2.10  from functional_groups import functional_groups
    2.11 @@ -306,6 +307,7 @@
    2.12                 
    2.13                  
    2.14  def find_connected_component(sequences, column):
    2.15 +    """output is a list of sequences"""
    2.16      #find one sequence that has blocks
    2.17      queue=[]
    2.18      output=[]
    2.19 @@ -334,9 +336,39 @@
    2.20  
    2.21  #[gblock for gblock in blocks_in_previous_column for gcolumn in gblock.columns if alignment.columns.index(gcolumn)==147]   
    2.22  
    2.23 -def create_table_of_homology(alignment, output_file_name):
    2.24 -    """"""
    2.25 +def create_monomer_homology(alignment):
    2.26 +    """returns MonomerHomology object, given alignment
    2.27 +    with links created by create_links function"""
    2.28 +    #creating markup
    2.29 +    for sequence in alignment:
    2.30 +        sim=SequenceIndexMarkup(sequence)
    2.31 +        sim.refresh()
    2.32 +    #inferring classes_of_equivalence = homologous monomers = connected_components from links
    2.33 +    output = MonomerHomology()
    2.34 +    class_number = 0
    2.35 +    for column in alignment.columns:
    2.36 +        remaining_sequences = copy.copy(column.sequences)
    2.37 +        while remaining_sequences!=0:
    2.38 +            connected_component = find_connected_component(remaining_sequences, column)
    2.39 +            if connected_component = []:
    2.40 +                #create class for each of the remaining monomers
    2.41 +                for sequence in remaining_sequences:
    2.42 +                    subset_of_sequences_in_column2monomer_homology_data(output, class_number, [sequence], column)
    2.43 +                    class_number+=1
    2.44 +                    remaining_sequences.remove(sequence)
    2.45 +            else:
    2.46 +                subset_of_sequences_in_column2monomer_homology_data(output, class_number, connected_component, column)
    2.47 +                class_number+=1
    2.48 +                for sequence in connected_component:
    2.49 +                    remaining_sequences.remove(sequence)
    2.50 +    return output
    2.51  
    2.52  
    2.53 +def subset_of_sequences_in_column2monomer_homology_data(monomer_homology, class_number, sequences, column):
    2.54 +    monomer_homology.classes[class_number] = []
    2.55 +    for sequence in sequences:
    2.56 +        monomer_homology.classes[class_number].append((sequence.name, column[sequence].index+1))
    2.57 +        monomer_homology.monomer_ids[(sequence.name, column[sequence].index+1)] = class_number
    2.58 +
    2.59  if __name__== '__main__':
    2.60      main(open(sys.argv[1]))