forester/java/src/org/forester/surfacing/MinimalDomainomeCalculator.java

   1
   2 package org.forester.surfacing;
   3
   4 import java.io.BufferedWriter;
   5 import java.io.File;
   6 import java.io.FileWriter;
   7 import java.io.IOException;
   8 import java.util.ArrayList;
   9 import java.util.HashMap;
  10 import java.util.HashSet;
  11 import java.util.List;
  12 import java.util.Map;
  13 import java.util.Map.Entry;
  14 import java.util.Set;
  15 import java.util.SortedMap;
  16 import java.util.SortedSet;
  17 import java.util.TreeMap;
  18 import java.util.TreeSet;
  19
  20 import org.forester.application.surfacing;
  21 import org.forester.phylogeny.Phylogeny;
  22 import org.forester.phylogeny.PhylogenyNode;
  23 import org.forester.phylogeny.iterators.PhylogenyNodeIterator;
  24 import org.forester.protein.Domain;
  25 import org.forester.protein.Protein;
  26 import org.forester.species.BasicSpecies;
  27 import org.forester.species.Species;
  28 import org.forester.util.ForesterUtil;
  29
  30 public final class MinimalDomainomeCalculator {
  31
  32     static final public void calcDomainome( final Phylogeny tre,
  33                                             final SortedMap<Species, List<Protein>> protein_lists_per_species,
  34                                             final double ie_cutoff ) {
  35         if ( protein_lists_per_species == null || tre == null ) {
  36             throw new IllegalArgumentException( "argument is null" );
  37         }
  38         if ( protein_lists_per_species.size() < 2 ) {
  39             throw new IllegalArgumentException( "not enough genomes" );
  40         }
  41         for( final PhylogenyNodeIterator iter = tre.iteratorPostorder(); iter.hasNext(); ) {
  42             final PhylogenyNode node = iter.next();
  43             if ( node.isInternal() ) {
  44                 System.out.println();
  45                 if ( node.getNodeData().isHasTaxonomy() ) {
  46                     System.out.println( node.getNodeData().getTaxonomy().getScientificName() + ":" );
  47                 }
  48                 else {
  49                     System.out.println( node.getName() + ":" );
  50                 }
  51                 final List<PhylogenyNode> external_descs = node.getAllExternalDescendants();
  52                 final List<Set<String>> domains_per_genome_list = new ArrayList<Set<String>>();
  53                 for( final PhylogenyNode external_desc : external_descs ) {
  54                     final String code = external_desc.getNodeData().getTaxonomy().getTaxonomyCode();
  55                     System.out.print( code + " " );
  56                     final List<Protein> proteins_per_species = protein_lists_per_species
  57                             .get( new BasicSpecies( code ) );
  58                     if ( proteins_per_species != null ) {
  59                         final SortedSet<String> domains_per_genome = new TreeSet<String>();
  60                         for( final Protein protein : proteins_per_species ) {
  61                             List<Domain> domains = protein.getProteinDomains();
  62                             for( final Domain domain : domains ) {
  63                                 if ( ( domain.getPerDomainEvalue() <= ie_cutoff ) || ( ie_cutoff <= -1 ) ) {
  64                                     domains_per_genome.add( domain.getDomainId() );
  65                                 }
  66                             }
  67                         }
  68                         if ( domains_per_genome.size() > 0 ) {
  69                             domains_per_genome_list.add( domains_per_genome );
  70                         }
  71                     }
  72                 }
  73                 System.out.println();
  74                 if ( domains_per_genome_list.size() > 0 ) {
  75                     Set<String> intersection = calcIntersection( domains_per_genome_list );
  76                     System.out.println( intersection );
  77                 }
  78             }
  79         }
  80     }
  81
  82     static final public void calcOme( final boolean use_domain_architectures,
  83                                       final Phylogeny tre,
  84                                       final SortedMap<Species, List<Protein>> protein_lists_per_species,
  85                                       final String separator,
  86                                       final double ie_cutoff,
  87                                       final String outfile_base )
  88             throws IOException {
  89         final SortedMap<String, SortedSet<String>> species_to_das_map = new TreeMap<String, SortedSet<String>>();
  90         if ( protein_lists_per_species == null || tre == null ) {
  91             throw new IllegalArgumentException( "argument is null" );
  92         }
  93         if ( protein_lists_per_species.size() < 2 ) {
  94             throw new IllegalArgumentException( "not enough genomes" );
  95         }
  96         final String x;
  97         if ( use_domain_architectures ) {
  98             x = "DA";
  99         }
 100         else {
 101             x = "domain";
 102         }
 103         final File outfile = new File( outfile_base + "_minimal_" + x + "ome.tsv" );
 104         final File outfile_table = new File( outfile_base + "_minimal_" + x + "ome_matrix.tsv" );
 105         SurfacingUtil.checkForOutputFileWriteability( outfile );
 106         SurfacingUtil.checkForOutputFileWriteability( outfile_table );
 107         final BufferedWriter out = new BufferedWriter( new FileWriter( outfile ) );
 108         final BufferedWriter out_table = new BufferedWriter( new FileWriter( outfile_table ) );
 109         out.write( "SPECIES\tCOMMON NAME\tCODE\tRANK\t#EXT NODES\tEXT NODE CODES\t#DA\tDA" );
 110         out.write( ForesterUtil.LINE_SEPARATOR );
 111         for( final PhylogenyNodeIterator iter = tre.iteratorPostorder(); iter.hasNext(); ) {
 112             final PhylogenyNode node = iter.next();
 113             final String species_name = node.getNodeData().isHasTaxonomy()
 114                     ? node.getNodeData().getTaxonomy().getScientificName() : node.getName();
 115             final String common = node.getNodeData().isHasTaxonomy() ? node.getNodeData().getTaxonomy().getCommonName()
 116                     : "";
 117             final String tcode = node.getNodeData().isHasTaxonomy() ? node.getNodeData().getTaxonomy().getTaxonomyCode()
 118                     : "";
 119             final String rank = node.getNodeData().isHasTaxonomy() ? node.getNodeData().getTaxonomy().getRank() : "";
 120             out.write( species_name );
 121             if ( !ForesterUtil.isEmpty( common ) ) {
 122                 out.write( "\t" + common );
 123             }
 124             else {
 125                 out.write( "\t" );
 126             }
 127             if ( !ForesterUtil.isEmpty( tcode ) ) {
 128                 out.write( "\t" + tcode );
 129             }
 130             else {
 131                 out.write( "\t" );
 132             }
 133             if ( !ForesterUtil.isEmpty( rank ) ) {
 134                 out.write( "\t" + rank );
 135             }
 136             else {
 137                 out.write( "\t" );
 138             }
 139             final List<PhylogenyNode> external_descs = node.getAllExternalDescendants();
 140             if ( node.isInternal() ) {
 141                 out.write( "\t" + external_descs.size() + "\t" );
 142             }
 143             else {
 144                 out.write( "\t\t" );
 145             }
 146             final List<Set<String>> das_per_genome_list = new ArrayList<Set<String>>();
 147             boolean first = true;
 148             for( final PhylogenyNode external_desc : external_descs ) {
 149                 final String code = external_desc.getNodeData().getTaxonomy().getTaxonomyCode();
 150                 if ( node.isInternal() ) {
 151                     if ( first ) {
 152                         first = false;
 153                     }
 154                     else {
 155                         out.write( ", " );
 156                     }
 157                     out.write( code );
 158                 }
 159                 final List<Protein> proteins_per_species = protein_lists_per_species.get( new BasicSpecies( code ) );
 160                 if ( proteins_per_species != null ) {
 161                     final SortedSet<String> das_per_genome = new TreeSet<String>();
 162                     for( final Protein protein : proteins_per_species ) {
 163                         if ( use_domain_architectures ) {
 164                             final String da = protein.toDomainArchitectureString( separator, ie_cutoff );
 165                             das_per_genome.add( da );
 166                         }
 167                         else {
 168                             List<Domain> domains = protein.getProteinDomains();
 169                             for( final Domain domain : domains ) {
 170                                 if ( ( ie_cutoff <= -1 ) || ( domain.getPerDomainEvalue() <= ie_cutoff ) ) {
 171                                     das_per_genome.add( domain.getDomainId() );
 172                                 }
 173                             }
 174                         }
 175                     }
 176                     if ( das_per_genome.size() > 0 ) {
 177                         das_per_genome_list.add( das_per_genome );
 178                     }
 179                 }
 180             }
 181             if ( das_per_genome_list.size() > 0 ) {
 182                 SortedSet<String> intersection = calcIntersection( das_per_genome_list );
 183                 out.write( "\t" + intersection.size() + "\t" );
 184                 first = true;
 185                 for( final String s : intersection ) {
 186                     if ( first ) {
 187                         first = false;
 188                     }
 189                     else {
 190                         out.write( ", " );
 191                     }
 192                     out.write( s );
 193                 }
 194                 out.write( ForesterUtil.LINE_SEPARATOR );
 195                 species_to_das_map.put( species_name, intersection );
 196             }
 197         }
 198         final SortedSet<String> all_species_names = new TreeSet<String>();
 199         final SortedSet<String> all_das = new TreeSet<String>();
 200         for( final Entry<String, SortedSet<String>> e : species_to_das_map.entrySet() ) {
 201             all_species_names.add( e.getKey() );
 202             for( final String das : e.getValue() ) {
 203                 all_das.add( das );
 204             }
 205         }
 206         out_table.write( '\t' );
 207         boolean first = true;
 208         for( final String species_name : all_species_names ) {
 209             if ( first ) {
 210                 first = false;
 211             }
 212             else {
 213                 out_table.write( '\t' );
 214             }
 215             out_table.write( species_name );
 216         }
 217         out_table.write( ForesterUtil.LINE_SEPARATOR );
 218         for( final String das : all_das ) {
 219             out_table.write( das );
 220             out_table.write( '\t' );
 221             first = true;
 222             for( final String species_name : all_species_names ) {
 223                 if ( first ) {
 224                     first = false;
 225                 }
 226                 else {
 227                     out_table.write( '\t' );
 228                 }
 229                 if ( species_to_das_map.get( species_name ).contains( das ) ) {
 230                     out_table.write( '1' );
 231                 }
 232                 else {
 233                     out_table.write( '0' );
 234                 }
 235             }
 236             out_table.write( ForesterUtil.LINE_SEPARATOR );
 237         }
 238         out.flush();
 239         out.close();
 240         out_table.flush();
 241         out_table.close();
 242         ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote minimal DAome data to           : " + outfile );
 243         ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote minimal DAome data to (as table): " + outfile_table );
 244     }
 245
 246     static final public void calcDAome( final Phylogeny tre,
 247                                         final SortedMap<Species, List<Protein>> protein_lists_per_species,
 248                                         final String separator,
 249                                         final double ie_cutoff,
 250                                         final String outfile_base )
 251             throws IOException {
 252         final SortedMap<String, SortedSet<String>> species_to_das_map = new TreeMap<String, SortedSet<String>>();
 253         if ( protein_lists_per_species == null || tre == null ) {
 254             throw new IllegalArgumentException( "argument is null" );
 255         }
 256         if ( protein_lists_per_species.size() < 2 ) {
 257             throw new IllegalArgumentException( "not enough genomes" );
 258         }
 259         final File outfile = new File( outfile_base + "_minimal_daome.txt" );
 260         final File outfile_table = new File( outfile_base + "_minimal_daome.tsv" );
 261         SurfacingUtil.checkForOutputFileWriteability( outfile );
 262         SurfacingUtil.checkForOutputFileWriteability( outfile_table );
 263         final BufferedWriter out = new BufferedWriter( new FileWriter( outfile ) );
 264         final BufferedWriter out_table = new BufferedWriter( new FileWriter( outfile_table ) );
 265         out.write( "SPECIES\tCOMMON NAME\tCODE\tRANK\t#EXT NODES\tEXT NODE CODES\t#DA\tDA" );
 266         out.write( ForesterUtil.LINE_SEPARATOR );
 267         for( final PhylogenyNodeIterator iter = tre.iteratorPostorder(); iter.hasNext(); ) {
 268             final PhylogenyNode node = iter.next();
 269             final String species_name = node.getNodeData().isHasTaxonomy()
 270                     ? node.getNodeData().getTaxonomy().getScientificName() : node.getName();
 271             final String common = node.getNodeData().isHasTaxonomy() ? node.getNodeData().getTaxonomy().getCommonName()
 272                     : "";
 273             final String tcode = node.getNodeData().isHasTaxonomy() ? node.getNodeData().getTaxonomy().getTaxonomyCode()
 274                     : "";
 275             final String rank = node.getNodeData().isHasTaxonomy() ? node.getNodeData().getTaxonomy().getRank() : "";
 276             out.write( species_name );
 277             if ( !ForesterUtil.isEmpty( common ) ) {
 278                 out.write( "\t" + common );
 279             }
 280             else {
 281                 out.write( "\t" );
 282             }
 283             if ( !ForesterUtil.isEmpty( tcode ) ) {
 284                 out.write( "\t" + tcode );
 285             }
 286             else {
 287                 out.write( "\t" );
 288             }
 289             if ( !ForesterUtil.isEmpty( rank ) ) {
 290                 out.write( "\t" + rank );
 291             }
 292             else {
 293                 out.write( "\t" );
 294             }
 295             final List<PhylogenyNode> external_descs = node.getAllExternalDescendants();
 296             if ( node.isInternal() ) {
 297                 out.write( "\t" + external_descs.size() + "\t" );
 298             }
 299             else {
 300                 out.write( "\t\t" );
 301             }
 302             final List<Set<String>> das_per_genome_list = new ArrayList<Set<String>>();
 303             boolean first = true;
 304             for( final PhylogenyNode external_desc : external_descs ) {
 305                 final String code = external_desc.getNodeData().getTaxonomy().getTaxonomyCode();
 306                 if ( node.isInternal() ) {
 307                     if ( first ) {
 308                         first = false;
 309                     }
 310                     else {
 311                         out.write( ", " );
 312                     }
 313                     out.write( code );
 314                 }
 315                 final List<Protein> proteins_per_species = protein_lists_per_species.get( new BasicSpecies( code ) );
 316                 if ( proteins_per_species != null ) {
 317                     final SortedSet<String> das_per_genome = new TreeSet<String>();
 318                     for( final Protein protein : proteins_per_species ) {
 319                         final String da = protein.toDomainArchitectureString( separator, ie_cutoff );
 320                         das_per_genome.add( da );
 321                     }
 322                     if ( das_per_genome.size() > 0 ) {
 323                         das_per_genome_list.add( das_per_genome );
 324                     }
 325                 }
 326             }
 327             if ( das_per_genome_list.size() > 0 ) {
 328                 SortedSet<String> intersection = calcIntersection( das_per_genome_list );
 329                 out.write( "\t" + intersection.size() + "\t" );
 330                 first = true;
 331                 for( final String s : intersection ) {
 332                     if ( first ) {
 333                         first = false;
 334                     }
 335                     else {
 336                         out.write( ", " );
 337                     }
 338                     out.write( s );
 339                 }
 340                 out.write( ForesterUtil.LINE_SEPARATOR );
 341                 species_to_das_map.put( species_name, intersection );
 342             }
 343         }
 344         final SortedSet<String> all_species_names = new TreeSet<String>();
 345         final SortedSet<String> all_das = new TreeSet<String>();
 346         for( final Entry<String, SortedSet<String>> e : species_to_das_map.entrySet() ) {
 347             all_species_names.add( e.getKey() );
 348             for( final String das : e.getValue() ) {
 349                 all_das.add( das );
 350             }
 351         }
 352         out_table.write( '\t' );
 353         boolean first = true;
 354         for( final String species_name : all_species_names ) {
 355             if ( first ) {
 356                 first = false;
 357             }
 358             else {
 359                 out_table.write( '\t' );
 360             }
 361             out_table.write( species_name );
 362         }
 363         out_table.write( ForesterUtil.LINE_SEPARATOR );
 364         for( final String das : all_das ) {
 365             out_table.write( das );
 366             out_table.write( '\t' );
 367             first = true;
 368             for( final String species_name : all_species_names ) {
 369                 if ( first ) {
 370                     first = false;
 371                 }
 372                 else {
 373                     out_table.write( '\t' );
 374                 }
 375                 if ( species_to_das_map.get( species_name ).contains( das ) ) {
 376                     out_table.write( '1' );
 377                 }
 378                 else {
 379                     out_table.write( '0' );
 380                 }
 381             }
 382             out_table.write( ForesterUtil.LINE_SEPARATOR );
 383         }
 384         out.flush();
 385         out.close();
 386         out_table.flush();
 387         out_table.close();
 388         ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote minimal DAome data to           : " + outfile );
 389         ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote minimal DAome data to (as table): " + outfile_table );
 390     }
 391
 392     private final static SortedSet<String> calcIntersection( final List<Set<String>> features_per_genome_list ) {
 393         final Set<String> first = features_per_genome_list.get( 0 );
 394         final SortedSet<String> my_first = new TreeSet<String>();
 395         for( final String s : first ) {
 396             my_first.add( s );
 397         }
 398         for( int i = 1; i < features_per_genome_list.size(); ++i ) {
 399             my_first.retainAll( features_per_genome_list.get( i ) );
 400         }
 401         return my_first;
 402     }
 403
 404     public static void main( final String[] args ) {
 405         Set<String> a = new HashSet<String>();
 406         Set<String> b = new HashSet<String>();
 407         Set<String> c = new HashSet<String>();
 408         Set<String> d = new HashSet<String>();
 409         a.add( "x" );
 410         a.add( "b" );
 411         a.add( "c" );
 412         b.add( "a" );
 413         b.add( "b" );
 414         b.add( "c" );
 415         c.add( "a" );
 416         c.add( "b" );
 417         c.add( "c" );
 418         c.add( "c" );
 419         c.add( "f" );
 420         d.add( "a" );
 421         d.add( "c" );
 422         d.add( "d" );
 423         List<Set<String>> domains_per_genome_list = new ArrayList<Set<String>>();
 424         domains_per_genome_list.add( a );
 425         domains_per_genome_list.add( b );
 426         domains_per_genome_list.add( c );
 427         domains_per_genome_list.add( d );
 428         Set<String> x = calcIntersection( domains_per_genome_list );
 429         System.out.println( x );
 430     }
 431 }