forester/java/src/org/forester/surfacing/SurfacingUtil.java

   1 // $Id:
   2 //
   3 // FORESTER -- software libraries and applications
   4 // for evolutionary biology research and applications.
   5 //
   6 // Copyright (C) 2008-2009 Christian M. Zmasek
   7 // Copyright (C) 2008-2009 Burnham Institute for Medical Research
   8 // All rights reserved
   9 //
  10 // This library is free software; you can redistribute it and/or
  11 // modify it under the terms of the GNU Lesser General Public
  12 // License as published by the Free Software Foundation; either
  13 // version 2.1 of the License, or (at your option) any later version.
  14 //
  15 // This library is distributed in the hope that it will be useful,
  16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  18 // Lesser General Public License for more details.
  19 //
  20 // You should have received a copy of the GNU Lesser General Public
  21 // License along with this library; if not, write to the Free Software
  22 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
  23 //
  24 // Contact: phylosoft @ gmail . com
  25 // WWW: https://sites.google.com/site/cmzmasek/home/software/forester
  26
  27 package org.forester.surfacing;
  28
  29 import java.awt.Color;
  30 import java.io.BufferedWriter;
  31 import java.io.File;
  32 import java.io.FileWriter;
  33 import java.io.IOException;
  34 import java.io.Writer;
  35 import java.text.DecimalFormat;
  36 import java.text.NumberFormat;
  37 import java.util.ArrayList;
  38 import java.util.Arrays;
  39 import java.util.Collections;
  40 import java.util.Comparator;
  41 import java.util.HashMap;
  42 import java.util.HashSet;
  43 import java.util.Iterator;
  44 import java.util.List;
  45 import java.util.Map;
  46 import java.util.Map.Entry;
  47 import java.util.PriorityQueue;
  48 import java.util.Set;
  49 import java.util.SortedMap;
  50 import java.util.SortedSet;
  51 import java.util.TreeMap;
  52 import java.util.TreeSet;
  53 import java.util.regex.Matcher;
  54 import java.util.regex.Pattern;
  55
  56 import org.forester.application.surfacing;
  57 import org.forester.evoinference.distance.NeighborJoining;
  58 import org.forester.evoinference.matrix.character.BasicCharacterStateMatrix;
  59 import org.forester.evoinference.matrix.character.CharacterStateMatrix;
  60 import org.forester.evoinference.matrix.character.CharacterStateMatrix.BinaryStates;
  61 import org.forester.evoinference.matrix.character.CharacterStateMatrix.Format;
  62 import org.forester.evoinference.matrix.character.CharacterStateMatrix.GainLossStates;
  63 import org.forester.evoinference.matrix.distance.DistanceMatrix;
  64 import org.forester.go.GoId;
  65 import org.forester.go.GoNameSpace;
  66 import org.forester.go.GoTerm;
  67 import org.forester.go.PfamToGoMapping;
  68 import org.forester.io.parsers.nexus.NexusConstants;
  69 import org.forester.io.parsers.phyloxml.PhyloXmlUtil;
  70 import org.forester.io.parsers.util.ParserUtils;
  71 import org.forester.io.writers.PhylogenyWriter;
  72 import org.forester.phylogeny.Phylogeny;
  73 import org.forester.phylogeny.PhylogenyMethods;
  74 import org.forester.phylogeny.PhylogenyNode;
  75 import org.forester.phylogeny.PhylogenyNode.NH_CONVERSION_SUPPORT_VALUE_STYLE;
  76 import org.forester.phylogeny.data.BinaryCharacters;
  77 import org.forester.phylogeny.data.Confidence;
  78 import org.forester.phylogeny.data.Taxonomy;
  79 import org.forester.phylogeny.factories.ParserBasedPhylogenyFactory;
  80 import org.forester.phylogeny.iterators.PhylogenyNodeIterator;
  81 import org.forester.protein.BasicDomain;
  82 import org.forester.protein.BasicProtein;
  83 import org.forester.protein.BinaryDomainCombination;
  84 import org.forester.protein.Domain;
  85 import org.forester.protein.Protein;
  86 import org.forester.species.Species;
  87 import org.forester.surfacing.DomainSimilarity.PRINT_OPTION;
  88 import org.forester.surfacing.DomainSimilarityCalculator.Detailedness;
  89 import org.forester.surfacing.GenomeWideCombinableDomains.GenomeWideCombinableDomainsSortOrder;
  90 import org.forester.util.AsciiHistogram;
  91 import org.forester.util.BasicDescriptiveStatistics;
  92 import org.forester.util.BasicTable;
  93 import org.forester.util.BasicTableParser;
  94 import org.forester.util.CommandLineArguments;
  95 import org.forester.util.DescriptiveStatistics;
  96 import org.forester.util.ForesterUtil;
  97 import org.forester.util.TaxonomyColors;
  98 import org.forester.util.TaxonomyGroups;
  99
 100 public final class SurfacingUtil {
 101
 102     public final static Pattern              PATTERN_SP_STYLE_TAXONOMY        = Pattern.compile( "^[A-Z0-9]{3,5}$" );
 103     private final static Map<String, String> _TAXCODE_HEXCOLORSTRING_MAP      = new HashMap<String, String>();
 104     private final static Map<String, String> _TAXCODE_TAXGROUP_MAP            = new HashMap<String, String>();
 105     private static final Comparator<Domain>  ASCENDING_CONFIDENCE_VALUE_ORDER = new Comparator<Domain>() {
 106
 107                                                                                   @Override
 108                                                                                   public int compare( final Domain d1,
 109                                                                                                       final Domain d2 ) {
 110                                                                                       if ( d1.getPerDomainEvalue() < d2
 111                                                                                               .getPerDomainEvalue() ) {
 112                                                                                           return -1;
 113                                                                                       }
 114                                                                                       else if ( d1
 115                                                                                               .getPerDomainEvalue() > d2
 116                                                                                                       .getPerDomainEvalue() ) {
 117                                                                                           return 1;
 118                                                                                       }
 119                                                                                       else {
 120                                                                                           return d1.compareTo( d2 );
 121                                                                                       }
 122                                                                                   }
 123                                                                               };
 124     private final static NumberFormat        FORMATTER_3                      = new DecimalFormat( "0.000" );
 125
 126     private SurfacingUtil() {
 127         // Hidden constructor.
 128     }
 129
 130     public static void addAllBinaryDomainCombinationToSet( final GenomeWideCombinableDomains genome,
 131                                                            final SortedSet<BinaryDomainCombination> binary_domain_combinations ) {
 132         final SortedMap<String, CombinableDomains> all_cd = genome.getAllCombinableDomainsIds();
 133         for( final String domain_id : all_cd.keySet() ) {
 134             binary_domain_combinations.addAll( all_cd.get( domain_id ).toBinaryDomainCombinations() );
 135         }
 136     }
 137
 138     public static void addAllDomainIdsToSet( final GenomeWideCombinableDomains genome,
 139                                              final SortedSet<String> domain_ids ) {
 140         final SortedSet<String> domains = genome.getAllDomainIds();
 141         for( final String domain : domains ) {
 142             domain_ids.add( domain );
 143         }
 144     }
 145
 146     public static DescriptiveStatistics calculateDescriptiveStatisticsForMeanValues( final Set<DomainSimilarity> similarities ) {
 147         final DescriptiveStatistics stats = new BasicDescriptiveStatistics();
 148         for( final DomainSimilarity similarity : similarities ) {
 149             stats.addValue( similarity.getMeanSimilarityScore() );
 150         }
 151         return stats;
 152     }
 153
 154     public static void checkForOutputFileWriteability( final File outfile ) {
 155         final String error = ForesterUtil.isWritableFile( outfile );
 156         if ( !ForesterUtil.isEmpty( error ) ) {
 157             ForesterUtil.fatalError( surfacing.PRG_NAME, error );
 158         }
 159     }
 160
 161     public static void checkWriteabilityForPairwiseComparisons( final DomainSimilarity.PRINT_OPTION domain_similarity_print_option,
 162                                                                 final String[][] input_file_properties,
 163                                                                 final String automated_pairwise_comparison_suffix,
 164                                                                 final File outdir ) {
 165         for( int i = 0; i < input_file_properties.length; ++i ) {
 166             for( int j = 0; j < i; ++j ) {
 167                 final String species_i = input_file_properties[ i ][ 1 ];
 168                 final String species_j = input_file_properties[ j ][ 1 ];
 169                 String pairwise_similarities_output_file_str = surfacing.PAIRWISE_DOMAIN_COMPARISONS_PREFIX + species_i
 170                         + "_" + species_j + automated_pairwise_comparison_suffix;
 171                 switch ( domain_similarity_print_option ) {
 172                     case HTML:
 173                         if ( !pairwise_similarities_output_file_str.endsWith( ".html" ) ) {
 174                             pairwise_similarities_output_file_str += ".html";
 175                         }
 176                         break;
 177                 }
 178                 final String error = ForesterUtil
 179                         .isWritableFile( new File( outdir == null ? pairwise_similarities_output_file_str
 180                                 : outdir + ForesterUtil.FILE_SEPARATOR + pairwise_similarities_output_file_str ) );
 181                 if ( !ForesterUtil.isEmpty( error ) ) {
 182                     ForesterUtil.fatalError( surfacing.PRG_NAME, error );
 183                 }
 184             }
 185         }
 186     }
 187
 188     public static void collectChangedDomainCombinationsFromBinaryStatesMatrixAsListToFile( final CharacterStateMatrix<CharacterStateMatrix.GainLossStates> matrix,
 189                                                                                            final BinaryDomainCombination.DomainCombinationType dc_type,
 190                                                                                            final List<BinaryDomainCombination> all_binary_domains_combination_gained,
 191                                                                                            final boolean get_gains ) {
 192         final SortedSet<String> sorted_ids = new TreeSet<String>();
 193         for( int i = 0; i < matrix.getNumberOfIdentifiers(); ++i ) {
 194             sorted_ids.add( matrix.getIdentifier( i ) );
 195         }
 196         for( final String id : sorted_ids ) {
 197             for( int c = 0; c < matrix.getNumberOfCharacters(); ++c ) {
 198                 if ( ( get_gains && ( matrix.getState( id, c ) == CharacterStateMatrix.GainLossStates.GAIN ) )
 199                         || ( !get_gains
 200                                 && ( matrix.getState( id, c ) == CharacterStateMatrix.GainLossStates.LOSS ) ) ) {
 201                     if ( dc_type == BinaryDomainCombination.DomainCombinationType.DIRECTED_ADJACTANT ) {
 202                         all_binary_domains_combination_gained.add( AdjactantDirectedBinaryDomainCombination
 203                                 .obtainInstance( matrix.getCharacter( c ) ) );
 204                     }
 205                     else if ( dc_type == BinaryDomainCombination.DomainCombinationType.DIRECTED ) {
 206                         all_binary_domains_combination_gained
 207                                 .add( DirectedBinaryDomainCombination.obtainInstance( matrix.getCharacter( c ) ) );
 208                     }
 209                     else {
 210                         all_binary_domains_combination_gained
 211                                 .add( BasicBinaryDomainCombination.obtainInstance( matrix.getCharacter( c ) ) );
 212                     }
 213                 }
 214             }
 215         }
 216     }
 217
 218     public static Map<String, List<GoId>> createDomainIdToGoIdMap( final List<PfamToGoMapping> pfam_to_go_mappings ) {
 219         final Map<String, List<GoId>> domain_id_to_go_ids_map = new HashMap<String, List<GoId>>( pfam_to_go_mappings
 220                 .size() );
 221         for( final PfamToGoMapping pfam_to_go : pfam_to_go_mappings ) {
 222             if ( !domain_id_to_go_ids_map.containsKey( pfam_to_go.getKey() ) ) {
 223                 domain_id_to_go_ids_map.put( pfam_to_go.getKey(), new ArrayList<GoId>() );
 224             }
 225             domain_id_to_go_ids_map.get( pfam_to_go.getKey() ).add( pfam_to_go.getValue() );
 226         }
 227         return domain_id_to_go_ids_map;
 228     }
 229
 230     public static Map<String, Set<String>> createDomainIdToSecondaryFeaturesMap( final File secondary_features_map_file )
 231             throws IOException {
 232         final BasicTable<String> primary_table = BasicTableParser.parse( secondary_features_map_file, '\t' );
 233         final Map<String, Set<String>> map = new TreeMap<String, Set<String>>();
 234         for( int r = 0; r < primary_table.getNumberOfRows(); ++r ) {
 235             final String domain_id = primary_table.getValue( 0, r );
 236             if ( !map.containsKey( domain_id ) ) {
 237                 map.put( domain_id, new HashSet<String>() );
 238             }
 239             map.get( domain_id ).add( primary_table.getValue( 1, r ) );
 240         }
 241         return map;
 242     }
 243
 244     public static Phylogeny createNjTreeBasedOnMatrixToFile( final File nj_tree_outfile,
 245                                                              final DistanceMatrix distance ) {
 246         checkForOutputFileWriteability( nj_tree_outfile );
 247         final NeighborJoining nj = NeighborJoining.createInstance();
 248         final Phylogeny phylogeny = nj.execute( ( DistanceMatrix ) distance );
 249         phylogeny.setName( nj_tree_outfile.getName() );
 250         writePhylogenyToFile( phylogeny, nj_tree_outfile.toString() );
 251         return phylogeny;
 252     }
 253
 254     public static StringBuilder createParametersAsString( final boolean ignore_dufs,
 255                                                           final double ie_value_max,
 256                                                           final double fs_e_value_max,
 257                                                           final int max_allowed_overlap,
 258                                                           final boolean no_engulfing_overlaps,
 259                                                           final File cutoff_scores_file,
 260                                                           final BinaryDomainCombination.DomainCombinationType dc_type ) {
 261         final StringBuilder parameters_sb = new StringBuilder();
 262         parameters_sb.append( "iE-value: " + ie_value_max );
 263         parameters_sb.append( ", FS E-value: " + fs_e_value_max );
 264         if ( cutoff_scores_file != null ) {
 265             parameters_sb.append( ", Cutoff-scores-file: " + cutoff_scores_file );
 266         }
 267         else {
 268             parameters_sb.append( ", Cutoff-scores-file: not-set" );
 269         }
 270         if ( max_allowed_overlap != surfacing.MAX_ALLOWED_OVERLAP_DEFAULT ) {
 271             parameters_sb.append( ", Max-overlap: " + max_allowed_overlap );
 272         }
 273         else {
 274             parameters_sb.append( ", Max-overlap: not-set" );
 275         }
 276         if ( no_engulfing_overlaps ) {
 277             parameters_sb.append( ", Engulfing-overlaps: not-allowed" );
 278         }
 279         else {
 280             parameters_sb.append( ", Engulfing-overlaps: allowed" );
 281         }
 282         if ( ignore_dufs ) {
 283             parameters_sb.append( ", Ignore-dufs: true" );
 284         }
 285         else {
 286             parameters_sb.append( ", Ignore-dufs: false" );
 287         }
 288         parameters_sb.append( ", DC type (if applicable): " + dc_type );
 289         return parameters_sb;
 290     }
 291
 292     public static void createSplitWriters( final File out_dir,
 293                                            final String my_outfile,
 294                                            final Map<Character, Writer> split_writers )
 295             throws IOException {
 296         split_writers.put( 'a',
 297                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 298                                    + "_domains_A.html" ) ) );
 299         split_writers.put( 'b',
 300                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 301                                    + "_domains_B.html" ) ) );
 302         split_writers.put( 'c',
 303                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 304                                    + "_domains_C.html" ) ) );
 305         split_writers.put( 'd',
 306                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 307                                    + "_domains_D.html" ) ) );
 308         split_writers.put( 'e',
 309                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 310                                    + "_domains_E.html" ) ) );
 311         split_writers.put( 'f',
 312                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 313                                    + "_domains_F.html" ) ) );
 314         split_writers.put( 'g',
 315                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 316                                    + "_domains_G.html" ) ) );
 317         split_writers.put( 'h',
 318                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 319                                    + "_domains_H.html" ) ) );
 320         split_writers.put( 'i',
 321                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 322                                    + "_domains_I.html" ) ) );
 323         split_writers.put( 'j',
 324                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 325                                    + "_domains_J.html" ) ) );
 326         split_writers.put( 'k',
 327                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 328                                    + "_domains_K.html" ) ) );
 329         split_writers.put( 'l',
 330                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 331                                    + "_domains_L.html" ) ) );
 332         split_writers.put( 'm',
 333                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 334                                    + "_domains_M.html" ) ) );
 335         split_writers.put( 'n',
 336                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 337                                    + "_domains_N.html" ) ) );
 338         split_writers.put( 'o',
 339                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 340                                    + "_domains_O.html" ) ) );
 341         split_writers.put( 'p',
 342                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 343                                    + "_domains_P.html" ) ) );
 344         split_writers.put( 'q',
 345                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 346                                    + "_domains_Q.html" ) ) );
 347         split_writers.put( 'r',
 348                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 349                                    + "_domains_R.html" ) ) );
 350         split_writers.put( 's',
 351                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 352                                    + "_domains_S.html" ) ) );
 353         split_writers.put( 't',
 354                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 355                                    + "_domains_T.html" ) ) );
 356         split_writers.put( 'u',
 357                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 358                                    + "_domains_U.html" ) ) );
 359         split_writers.put( 'v',
 360                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 361                                    + "_domains_V.html" ) ) );
 362         split_writers.put( 'w',
 363                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 364                                    + "_domains_W.html" ) ) );
 365         split_writers.put( 'x',
 366                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 367                                    + "_domains_X.html" ) ) );
 368         split_writers.put( 'y',
 369                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 370                                    + "_domains_Y.html" ) ) );
 371         split_writers.put( 'z',
 372                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 373                                    + "_domains_Z.html" ) ) );
 374         split_writers.put( '0',
 375                            new BufferedWriter( new FileWriter( out_dir + ForesterUtil.FILE_SEPARATOR + my_outfile
 376                                    + "_domains_0.html" ) ) );
 377     }
 378
 379     public static Map<String, Integer> createTaxCodeToIdMap( final Phylogeny phy ) {
 380         final Map<String, Integer> m = new HashMap<String, Integer>();
 381         for( final PhylogenyNodeIterator iter = phy.iteratorExternalForward(); iter.hasNext(); ) {
 382             final PhylogenyNode n = iter.next();
 383             if ( n.getNodeData().isHasTaxonomy() ) {
 384                 final Taxonomy t = n.getNodeData().getTaxonomy();
 385                 final String c = t.getTaxonomyCode();
 386                 if ( !ForesterUtil.isEmpty( c ) ) {
 387                     if ( n.getNodeData().getTaxonomy() == null ) {
 388                         ForesterUtil.fatalError( surfacing.PRG_NAME, "no taxonomy id for node " + n );
 389                     }
 390                     final String id = n.getNodeData().getTaxonomy().getIdentifier().getValue();
 391                     if ( ForesterUtil.isEmpty( id ) ) {
 392                         ForesterUtil.fatalError( surfacing.PRG_NAME, "no taxonomy id for node " + n );
 393                     }
 394                     if ( m.containsKey( c ) ) {
 395                         ForesterUtil.fatalError( surfacing.PRG_NAME, "taxonomy code " + c + " is not unique" );
 396                     }
 397                     final int iid = Integer.valueOf( id );
 398                     if ( m.containsValue( iid ) ) {
 399                         ForesterUtil.fatalError( surfacing.PRG_NAME, "taxonomy id " + iid + " is not unique" );
 400                     }
 401                     m.put( c, iid );
 402                 }
 403             }
 404             else {
 405                 ForesterUtil.fatalError( surfacing.PRG_NAME, "no taxonomy for node " + n );
 406             }
 407         }
 408         return m;
 409     }
 410
 411     public static void decoratePrintableDomainSimilarities( final SortedSet<DomainSimilarity> domain_similarities,
 412                                                             final Detailedness detailedness ) {
 413         for( final DomainSimilarity domain_similarity : domain_similarities ) {
 414             if ( domain_similarity instanceof DomainSimilarity ) {
 415                 final DomainSimilarity printable_domain_similarity = domain_similarity;
 416                 printable_domain_similarity.setDetailedness( detailedness );
 417             }
 418         }
 419     }
 420
 421     public static void doit( final List<Protein> proteins,
 422                              final List<String> query_domain_ids_nc_order,
 423                              final Writer out,
 424                              final String separator,
 425                              final String limit_to_species,
 426                              final Map<String, List<Integer>> average_protein_lengths_by_dc )
 427             throws IOException {
 428         for( final Protein protein : proteins ) {
 429             if ( ForesterUtil.isEmpty( limit_to_species )
 430                     || protein.getSpecies().getSpeciesId().equalsIgnoreCase( limit_to_species ) ) {
 431                 if ( protein.contains( query_domain_ids_nc_order, true ) ) {
 432                     out.write( protein.getSpecies().getSpeciesId() );
 433                     out.write( separator );
 434                     out.write( protein.getProteinId().getId() );
 435                     out.write( separator );
 436                     out.write( "[" );
 437                     final Set<String> visited_domain_ids = new HashSet<String>();
 438                     boolean first = true;
 439                     for( final Domain domain : protein.getProteinDomains() ) {
 440                         if ( !visited_domain_ids.contains( domain.getDomainId() ) ) {
 441                             visited_domain_ids.add( domain.getDomainId() );
 442                             if ( first ) {
 443                                 first = false;
 444                             }
 445                             else {
 446                                 out.write( " " );
 447                             }
 448                             out.write( domain.getDomainId() );
 449                             out.write( " {" );
 450                             out.write( "" + domain.getTotalCount() );
 451                             out.write( "}" );
 452                         }
 453                     }
 454                     out.write( "]" );
 455                     out.write( separator );
 456                     if ( !( ForesterUtil.isEmpty( protein.getDescription() )
 457                             || protein.getDescription().equals( SurfacingConstants.NONE ) ) ) {
 458                         out.write( protein.getDescription() );
 459                     }
 460                     out.write( separator );
 461                     if ( !( ForesterUtil.isEmpty( protein.getAccession() )
 462                             || protein.getAccession().equals( SurfacingConstants.NONE ) ) ) {
 463                         out.write( protein.getAccession() );
 464                     }
 465                     out.write( SurfacingConstants.NL );
 466                 }
 467             }
 468         }
 469         out.flush();
 470     }
 471
 472     public static void domainsPerProteinsStatistics( final String genome,
 473                                                      final List<Protein> protein_list,
 474                                                      final DescriptiveStatistics all_genomes_domains_per_potein_stats,
 475                                                      final SortedMap<Integer, Integer> all_genomes_domains_per_potein_histo,
 476                                                      final SortedSet<String> domains_which_are_always_single,
 477                                                      final SortedSet<String> domains_which_are_sometimes_single_sometimes_not,
 478                                                      final SortedSet<String> domains_which_never_single,
 479                                                      final Writer writer ) {
 480         final DescriptiveStatistics stats = new BasicDescriptiveStatistics();
 481         for( final Protein protein : protein_list ) {
 482             final int domains = protein.getNumberOfProteinDomains();
 483             //System.out.println( domains );
 484             stats.addValue( domains );
 485             all_genomes_domains_per_potein_stats.addValue( domains );
 486             if ( !all_genomes_domains_per_potein_histo.containsKey( domains ) ) {
 487                 all_genomes_domains_per_potein_histo.put( domains, 1 );
 488             }
 489             else {
 490                 all_genomes_domains_per_potein_histo.put( domains,
 491                                                           1 + all_genomes_domains_per_potein_histo.get( domains ) );
 492             }
 493             if ( domains == 1 ) {
 494                 final String domain = protein.getProteinDomain( 0 ).getDomainId();
 495                 if ( !domains_which_are_sometimes_single_sometimes_not.contains( domain ) ) {
 496                     if ( domains_which_never_single.contains( domain ) ) {
 497                         domains_which_never_single.remove( domain );
 498                         domains_which_are_sometimes_single_sometimes_not.add( domain );
 499                     }
 500                     else {
 501                         domains_which_are_always_single.add( domain );
 502                     }
 503                 }
 504             }
 505             else if ( domains > 1 ) {
 506                 for( final Domain d : protein.getProteinDomains() ) {
 507                     final String domain = d.getDomainId();
 508                     // System.out.println( domain );
 509                     if ( !domains_which_are_sometimes_single_sometimes_not.contains( domain ) ) {
 510                         if ( domains_which_are_always_single.contains( domain ) ) {
 511                             domains_which_are_always_single.remove( domain );
 512                             domains_which_are_sometimes_single_sometimes_not.add( domain );
 513                         }
 514                         else {
 515                             domains_which_never_single.add( domain );
 516                         }
 517                     }
 518                 }
 519             }
 520         }
 521         try {
 522             writer.write( genome );
 523             writer.write( "\t" );
 524             if ( stats.getN() >= 1 ) {
 525                 writer.write( stats.arithmeticMean() + "" );
 526                 writer.write( "\t" );
 527                 if ( stats.getN() >= 2 ) {
 528                     writer.write( stats.sampleStandardDeviation() + "" );
 529                 }
 530                 else {
 531                     writer.write( "" );
 532                 }
 533                 writer.write( "\t" );
 534                 writer.write( stats.median() + "" );
 535                 writer.write( "\t" );
 536                 writer.write( stats.getN() + "" );
 537                 writer.write( "\t" );
 538                 writer.write( stats.getMin() + "" );
 539                 writer.write( "\t" );
 540                 writer.write( stats.getMax() + "" );
 541             }
 542             else {
 543                 writer.write( "\t" );
 544                 writer.write( "\t" );
 545                 writer.write( "\t" );
 546                 writer.write( "0" );
 547                 writer.write( "\t" );
 548                 writer.write( "\t" );
 549             }
 550             writer.write( "\n" );
 551         }
 552         catch ( final IOException e ) {
 553             e.printStackTrace();
 554         }
 555     }
 556
 557     public static void executeDomainLengthAnalysis( final String[][] input_file_properties,
 558                                                     final int number_of_genomes,
 559                                                     final DomainLengthsTable domain_lengths_table,
 560                                                     final File outfile )
 561             throws IOException {
 562         final DecimalFormat df = new DecimalFormat( "#.00" );
 563         checkForOutputFileWriteability( outfile );
 564         final BufferedWriter out = new BufferedWriter( new FileWriter( outfile ) );
 565         out.write( "MEAN BASED STATISTICS PER SPECIES" );
 566         out.write( ForesterUtil.LINE_SEPARATOR );
 567         out.write( domain_lengths_table.createMeanBasedStatisticsPerSpeciesTable().toString() );
 568         out.write( ForesterUtil.LINE_SEPARATOR );
 569         out.write( ForesterUtil.LINE_SEPARATOR );
 570         final List<DomainLengths> domain_lengths_list = domain_lengths_table.getDomainLengthsList();
 571         out.write( "OUTLIER SPECIES PER DOMAIN (Z>=1.5)" );
 572         out.write( ForesterUtil.LINE_SEPARATOR );
 573         for( final DomainLengths domain_lengths : domain_lengths_list ) {
 574             final List<Species> species_list = domain_lengths.getMeanBasedOutlierSpecies( 1.5 );
 575             if ( species_list.size() > 0 ) {
 576                 out.write( domain_lengths.getDomainId() + "\t" );
 577                 for( final Species species : species_list ) {
 578                     out.write( species + "\t" );
 579                 }
 580                 out.write( ForesterUtil.LINE_SEPARATOR );
 581             }
 582         }
 583         out.write( ForesterUtil.LINE_SEPARATOR );
 584         out.write( ForesterUtil.LINE_SEPARATOR );
 585         out.write( "OUTLIER SPECIES (Z 1.0)" );
 586         out.write( ForesterUtil.LINE_SEPARATOR );
 587         final DescriptiveStatistics stats_for_all_species = domain_lengths_table
 588                 .calculateMeanBasedStatisticsForAllSpecies();
 589         out.write( stats_for_all_species.asSummary() );
 590         out.write( ForesterUtil.LINE_SEPARATOR );
 591         final AsciiHistogram histo = new AsciiHistogram( stats_for_all_species );
 592         out.write( histo.toStringBuffer( 40, '=', 60, 4 ).toString() );
 593         out.write( ForesterUtil.LINE_SEPARATOR );
 594         final double population_sd = stats_for_all_species.sampleStandardDeviation();
 595         final double population_mean = stats_for_all_species.arithmeticMean();
 596         for( final Species species : domain_lengths_table.getSpecies() ) {
 597             final double x = domain_lengths_table.calculateMeanBasedStatisticsForSpecies( species ).arithmeticMean();
 598             final double z = ( x - population_mean ) / population_sd;
 599             out.write( species + "\t" + z );
 600             out.write( ForesterUtil.LINE_SEPARATOR );
 601         }
 602         out.write( ForesterUtil.LINE_SEPARATOR );
 603         for( final Species species : domain_lengths_table.getSpecies() ) {
 604             final DescriptiveStatistics stats_for_species = domain_lengths_table
 605                     .calculateMeanBasedStatisticsForSpecies( species );
 606             final double x = stats_for_species.arithmeticMean();
 607             final double z = ( x - population_mean ) / population_sd;
 608             if ( ( z <= -1.0 ) || ( z >= 1.0 ) ) {
 609                 out.write( species + "\t" + df.format( z ) + "\t" + stats_for_species.asSummary() );
 610                 out.write( ForesterUtil.LINE_SEPARATOR );
 611             }
 612         }
 613         out.close();
 614         System.gc();
 615     }
 616
 617     /**
 618      * Warning: This side-effects 'all_bin_domain_combinations_encountered'!
 619      *
 620      *
 621      * @param output_file
 622      * @param all_bin_domain_combinations_changed
 623      * @param sum_of_all_domains_encountered
 624      * @param all_bin_domain_combinations_encountered
 625      * @param is_gains_analysis
 626      * @param protein_length_stats_by_dc
 627      * @throws IOException
 628      */
 629     public static void executeFitchGainsAnalysis( final File output_file,
 630                                                   final List<BinaryDomainCombination> all_bin_domain_combinations_changed,
 631                                                   final int sum_of_all_domains_encountered,
 632                                                   final SortedSet<BinaryDomainCombination> all_bin_domain_combinations_encountered,
 633                                                   final boolean is_gains_analysis )
 634             throws IOException {
 635         checkForOutputFileWriteability( output_file );
 636         final Writer out = ForesterUtil.createBufferedWriter( output_file );
 637         final SortedMap<Object, Integer> bdc_to_counts = ForesterUtil
 638                 .listToSortedCountsMap( all_bin_domain_combinations_changed );
 639         final SortedSet<String> all_domains_in_combination_changed_more_than_once = new TreeSet<String>();
 640         final SortedSet<String> all_domains_in_combination_changed_only_once = new TreeSet<String>();
 641         int above_one = 0;
 642         int one = 0;
 643         for( final Object bdc_object : bdc_to_counts.keySet() ) {
 644             final BinaryDomainCombination bdc = ( BinaryDomainCombination ) bdc_object;
 645             final int count = bdc_to_counts.get( bdc_object );
 646             if ( count < 1 ) {
 647                 ForesterUtil.unexpectedFatalError( surfacing.PRG_NAME, "count < 1 " );
 648             }
 649             out.write( bdc + "\t" + count + ForesterUtil.LINE_SEPARATOR );
 650             if ( count > 1 ) {
 651                 all_domains_in_combination_changed_more_than_once.add( bdc.getId0() );
 652                 all_domains_in_combination_changed_more_than_once.add( bdc.getId1() );
 653                 above_one++;
 654             }
 655             else if ( count == 1 ) {
 656                 all_domains_in_combination_changed_only_once.add( bdc.getId0() );
 657                 all_domains_in_combination_changed_only_once.add( bdc.getId1() );
 658                 one++;
 659             }
 660         }
 661         final int all = all_bin_domain_combinations_encountered.size();
 662         int never_lost = -1;
 663         if ( !is_gains_analysis ) {
 664             all_bin_domain_combinations_encountered.removeAll( all_bin_domain_combinations_changed );
 665             never_lost = all_bin_domain_combinations_encountered.size();
 666             for( final BinaryDomainCombination bdc : all_bin_domain_combinations_encountered ) {
 667                 out.write( bdc + "\t" + "0" + ForesterUtil.LINE_SEPARATOR );
 668             }
 669         }
 670         if ( is_gains_analysis ) {
 671             out.write( "Sum of all distinct domain combinations appearing once               : " + one
 672                     + ForesterUtil.LINE_SEPARATOR );
 673             out.write( "Sum of all distinct domain combinations appearing more than once     : " + above_one
 674                     + ForesterUtil.LINE_SEPARATOR );
 675             out.write( "Sum of all distinct domains in combinations apppearing only once     : "
 676                     + all_domains_in_combination_changed_only_once.size() + ForesterUtil.LINE_SEPARATOR );
 677             out.write( "Sum of all distinct domains in combinations apppearing more than once: "
 678                     + all_domains_in_combination_changed_more_than_once.size() + ForesterUtil.LINE_SEPARATOR );
 679         }
 680         else {
 681             out.write( "Sum of all distinct domain combinations never lost                   : " + never_lost
 682                     + ForesterUtil.LINE_SEPARATOR );
 683             out.write( "Sum of all distinct domain combinations lost once                    : " + one
 684                     + ForesterUtil.LINE_SEPARATOR );
 685             out.write( "Sum of all distinct domain combinations lost more than once          : " + above_one
 686                     + ForesterUtil.LINE_SEPARATOR );
 687             out.write( "Sum of all distinct domains in combinations lost only once           : "
 688                     + all_domains_in_combination_changed_only_once.size() + ForesterUtil.LINE_SEPARATOR );
 689             out.write( "Sum of all distinct domains in combinations lost more than once: "
 690                     + all_domains_in_combination_changed_more_than_once.size() + ForesterUtil.LINE_SEPARATOR );
 691         }
 692         out.write( "All binary combinations                                              : " + all
 693                 + ForesterUtil.LINE_SEPARATOR );
 694         out.write( "All domains                                                          : "
 695                 + sum_of_all_domains_encountered );
 696         out.close();
 697         ForesterUtil
 698                 .programMessage( surfacing.PRG_NAME,
 699                                  "Wrote fitch domain combination dynamics counts analysis to \"" + output_file + "\"" );
 700     }
 701
 702     /**
 703      *
 704      * @param all_binary_domains_combination_lost_fitch
 705      * @param use_last_in_fitch_parsimony
 706      * @param perform_dc_fich
 707      * @param consider_directedness_and_adjacency_for_bin_combinations
 708      * @param all_binary_domains_combination_gained if null ignored, otherwise this is to list all binary domain combinations
 709      * which were gained under unweighted (Fitch) parsimony.
 710      */
 711     public static void executeParsimonyAnalysis( final long random_number_seed_for_fitch_parsimony,
 712                                                  final boolean radomize_fitch_parsimony,
 713                                                  final String outfile_name,
 714                                                  final DomainParsimonyCalculator domain_parsimony,
 715                                                  final Phylogeny phylogeny,
 716                                                  final Map<String, List<GoId>> domain_id_to_go_ids_map,
 717                                                  final Map<GoId, GoTerm> go_id_to_term_map,
 718                                                  final GoNameSpace go_namespace_limit,
 719                                                  final String parameters_str,
 720                                                  final Map<String, Set<String>>[] domain_id_to_secondary_features_maps,
 721                                                  final SortedSet<String> positive_filter,
 722                                                  final boolean output_binary_domain_combinations_for_graphs,
 723                                                  final List<BinaryDomainCombination> all_binary_domains_combination_gained_fitch,
 724                                                  final List<BinaryDomainCombination> all_binary_domains_combination_lost_fitch,
 725                                                  final BinaryDomainCombination.DomainCombinationType dc_type,
 726                                                  final Map<String, DescriptiveStatistics> protein_length_stats_by_dc,
 727                                                  final Map<String, DescriptiveStatistics> domain_number_stats_by_dc,
 728                                                  final Map<String, DescriptiveStatistics> domain_length_stats_by_domain,
 729                                                  final Map<String, Integer> tax_code_to_id_map,
 730                                                  final boolean write_to_nexus,
 731                                                  final boolean use_last_in_fitch_parsimony,
 732                                                  final boolean perform_dc_fich ) {
 733         final String sep = ForesterUtil.LINE_SEPARATOR + "###################" + ForesterUtil.LINE_SEPARATOR;
 734         final String date_time = ForesterUtil.getCurrentDateTime();
 735         final SortedSet<String> all_pfams_encountered = new TreeSet<String>();
 736         final SortedSet<String> all_pfams_gained_as_domains = new TreeSet<String>();
 737         final SortedSet<String> all_pfams_lost_as_domains = new TreeSet<String>();
 738         final SortedSet<String> all_pfams_gained_as_dom_combinations = new TreeSet<String>();
 739         final SortedSet<String> all_pfams_lost_as_dom_combinations = new TreeSet<String>();
 740         if ( write_to_nexus ) {
 741             writeToNexus( outfile_name, domain_parsimony, phylogeny );
 742         }
 743         // DOLLO DOMAINS
 744         // -------------
 745         Phylogeny local_phylogeny_l = phylogeny.copy();
 746         if ( ( positive_filter != null ) && ( positive_filter.size() > 0 ) ) {
 747             domain_parsimony.executeDolloParsimonyOnDomainPresence( positive_filter );
 748         }
 749         else {
 750             domain_parsimony.executeDolloParsimonyOnDomainPresence();
 751         }
 752         SurfacingUtil.writeMatrixToFile( domain_parsimony.getGainLossMatrix(),
 753                                          outfile_name + surfacing.PARSIMONY_OUTPUT_GL_SUFFIX_DOLLO_DOMAINS,
 754                                          Format.FORESTER );
 755         SurfacingUtil.writeMatrixToFile( domain_parsimony.getGainLossCountsMatrix(),
 756                                          outfile_name + surfacing.PARSIMONY_OUTPUT_GL_COUNTS_SUFFIX_DOLLO_DOMAINS,
 757                                          Format.FORESTER );
 758         SurfacingUtil.writeBinaryStatesMatrixAsListToFile( domain_parsimony.getGainLossMatrix(),
 759                                                            CharacterStateMatrix.GainLossStates.GAIN,
 760                                                            outfile_name + surfacing.PARSIMONY_OUTPUT_DOLLO_GAINS_D,
 761                                                            sep,
 762                                                            ForesterUtil.LINE_SEPARATOR,
 763                                                            null );
 764         SurfacingUtil.writeBinaryStatesMatrixAsListToFile( domain_parsimony.getGainLossMatrix(),
 765                                                            CharacterStateMatrix.GainLossStates.LOSS,
 766                                                            outfile_name + surfacing.PARSIMONY_OUTPUT_DOLLO_LOSSES_D,
 767                                                            sep,
 768                                                            ForesterUtil.LINE_SEPARATOR,
 769                                                            null );
 770         SurfacingUtil.writeBinaryStatesMatrixAsListToFile( domain_parsimony.getGainLossMatrix(),
 771                                                            null,
 772                                                            outfile_name + surfacing.PARSIMONY_OUTPUT_DOLLO_PRESENT_D,
 773                                                            sep,
 774                                                            ForesterUtil.LINE_SEPARATOR,
 775                                                            null );
 776         //HTML:
 777         writeBinaryStatesMatrixToList( domain_id_to_go_ids_map,
 778                                        go_id_to_term_map,
 779                                        go_namespace_limit,
 780                                        false,
 781                                        domain_parsimony.getGainLossMatrix(),
 782                                        CharacterStateMatrix.GainLossStates.GAIN,
 783                                        outfile_name + surfacing.PARSIMONY_OUTPUT_DOLLO_GAINS_HTML_D,
 784                                        sep,
 785                                        ForesterUtil.LINE_SEPARATOR,
 786                                        "Dollo Parsimony | Gains | Domains",
 787                                        "+",
 788                                        domain_id_to_secondary_features_maps,
 789                                        all_pfams_encountered,
 790                                        all_pfams_gained_as_domains,
 791                                        "_dollo_gains_d",
 792                                        tax_code_to_id_map );
 793         writeBinaryStatesMatrixToList( domain_id_to_go_ids_map,
 794                                        go_id_to_term_map,
 795                                        go_namespace_limit,
 796                                        false,
 797                                        domain_parsimony.getGainLossMatrix(),
 798                                        CharacterStateMatrix.GainLossStates.LOSS,
 799                                        outfile_name + surfacing.PARSIMONY_OUTPUT_DOLLO_LOSSES_HTML_D,
 800                                        sep,
 801                                        ForesterUtil.LINE_SEPARATOR,
 802                                        "Dollo Parsimony | Losses | Domains",
 803                                        "-",
 804                                        domain_id_to_secondary_features_maps,
 805                                        all_pfams_encountered,
 806                                        all_pfams_lost_as_domains,
 807                                        "_dollo_losses_d",
 808                                        tax_code_to_id_map );
 809         writeBinaryStatesMatrixToList( domain_id_to_go_ids_map,
 810                                        go_id_to_term_map,
 811                                        go_namespace_limit,
 812                                        false,
 813                                        domain_parsimony.getGainLossMatrix(),
 814                                        null,
 815                                        outfile_name + surfacing.PARSIMONY_OUTPUT_DOLLO_PRESENT_HTML_D,
 816                                        sep,
 817                                        ForesterUtil.LINE_SEPARATOR,
 818                                        "Dollo Parsimony | Present | Domains",
 819                                        "",
 820                                        domain_id_to_secondary_features_maps,
 821                                        all_pfams_encountered,
 822                                        null,
 823                                        "_dollo_present_d",
 824                                        tax_code_to_id_map );
 825         preparePhylogeny( local_phylogeny_l,
 826                           domain_parsimony,
 827                           date_time,
 828                           "Dollo parsimony on domain presence/absence",
 829                           "dollo_on_domains_" + outfile_name,
 830                           parameters_str );
 831         SurfacingUtil.writePhylogenyToFile( local_phylogeny_l,
 832                                             outfile_name + surfacing.DOMAINS_PARSIMONY_TREE_OUTPUT_SUFFIX_DOLLO );
 833         try {
 834             writeAllDomainsChangedOnAllSubtrees( local_phylogeny_l, true, outfile_name, "_dollo_all_gains_d" );
 835             writeAllDomainsChangedOnAllSubtrees( local_phylogeny_l, false, outfile_name, "_dollo_all_losses_d" );
 836         }
 837         catch ( final IOException e ) {
 838             e.printStackTrace();
 839             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getLocalizedMessage() );
 840         }
 841         if ( perform_dc_fich && ( domain_parsimony.calculateNumberOfBinaryDomainCombination() > 0 ) ) {
 842             // FITCH DOMAIN COMBINATIONS
 843             // -------------------------
 844             local_phylogeny_l = phylogeny.copy();
 845             String randomization = "no";
 846             if ( radomize_fitch_parsimony ) {
 847                 domain_parsimony
 848                         .executeFitchParsimonyOnBinaryDomainCombintion( random_number_seed_for_fitch_parsimony );
 849                 randomization = "yes, seed = " + random_number_seed_for_fitch_parsimony;
 850             }
 851             else {
 852                 domain_parsimony.executeFitchParsimonyOnBinaryDomainCombintion( use_last_in_fitch_parsimony );
 853             }
 854             SurfacingUtil.writeMatrixToFile( domain_parsimony.getGainLossMatrix(),
 855                                              outfile_name
 856                                                      + surfacing.PARSIMONY_OUTPUT_GL_SUFFIX_FITCH_BINARY_COMBINATIONS,
 857                                              Format.FORESTER );
 858             SurfacingUtil.writeMatrixToFile( domain_parsimony.getGainLossCountsMatrix(),
 859                                              outfile_name
 860                                                      + surfacing.PARSIMONY_OUTPUT_GL_COUNTS_SUFFIX_FITCH_BINARY_COMBINATIONS,
 861                                              Format.FORESTER );
 862             SurfacingUtil.writeBinaryStatesMatrixAsListToFile( domain_parsimony.getGainLossMatrix(),
 863                                                                CharacterStateMatrix.GainLossStates.GAIN,
 864                                                                outfile_name + surfacing.PARSIMONY_OUTPUT_FITCH_GAINS_BC,
 865                                                                sep,
 866                                                                ForesterUtil.LINE_SEPARATOR,
 867                                                                null );
 868             SurfacingUtil.writeBinaryStatesMatrixAsListToFile( domain_parsimony.getGainLossMatrix(),
 869                                                                CharacterStateMatrix.GainLossStates.LOSS,
 870                                                                outfile_name
 871                                                                        + surfacing.PARSIMONY_OUTPUT_FITCH_LOSSES_BC,
 872                                                                sep,
 873                                                                ForesterUtil.LINE_SEPARATOR,
 874                                                                null );
 875             SurfacingUtil.writeBinaryStatesMatrixAsListToFile( domain_parsimony.getGainLossMatrix(),
 876                                                                null,
 877                                                                outfile_name
 878                                                                        + surfacing.PARSIMONY_OUTPUT_FITCH_PRESENT_BC,
 879                                                                sep,
 880                                                                ForesterUtil.LINE_SEPARATOR,
 881                                                                null );
 882             if ( all_binary_domains_combination_gained_fitch != null ) {
 883                 collectChangedDomainCombinationsFromBinaryStatesMatrixAsListToFile( domain_parsimony
 884                         .getGainLossMatrix(), dc_type, all_binary_domains_combination_gained_fitch, true );
 885             }
 886             if ( all_binary_domains_combination_lost_fitch != null ) {
 887                 collectChangedDomainCombinationsFromBinaryStatesMatrixAsListToFile( domain_parsimony
 888                         .getGainLossMatrix(), dc_type, all_binary_domains_combination_lost_fitch, false );
 889             }
 890             if ( output_binary_domain_combinations_for_graphs ) {
 891                 SurfacingUtil.writeBinaryStatesMatrixAsListToFileForBinaryCombinationsForGraphAnalysis( domain_parsimony
 892                         .getGainLossMatrix(),
 893                                                                                                         null,
 894                                                                                                         outfile_name
 895                                                                                                                 + surfacing.PARSIMONY_OUTPUT_FITCH_PRESENT_BC_OUTPUTFILE_SUFFIX_FOR_GRAPH_ANALYSIS,
 896                                                                                                         sep,
 897                                                                                                         ForesterUtil.LINE_SEPARATOR,
 898                                                                                                         BinaryDomainCombination.OutputFormat.DOT );
 899             }
 900             // HTML:
 901             writeBinaryStatesMatrixToList( domain_id_to_go_ids_map,
 902                                            go_id_to_term_map,
 903                                            go_namespace_limit,
 904                                            true,
 905                                            domain_parsimony.getGainLossMatrix(),
 906                                            CharacterStateMatrix.GainLossStates.GAIN,
 907                                            outfile_name + surfacing.PARSIMONY_OUTPUT_FITCH_GAINS_HTML_BC,
 908                                            sep,
 909                                            ForesterUtil.LINE_SEPARATOR,
 910                                            "Fitch Parsimony | Gains | Domain Combinations",
 911                                            "+",
 912                                            null,
 913                                            all_pfams_encountered,
 914                                            all_pfams_gained_as_dom_combinations,
 915                                            "_fitch_gains_dc",
 916                                            tax_code_to_id_map );
 917             writeBinaryStatesMatrixToList( domain_id_to_go_ids_map,
 918                                            go_id_to_term_map,
 919                                            go_namespace_limit,
 920                                            true,
 921                                            domain_parsimony.getGainLossMatrix(),
 922                                            CharacterStateMatrix.GainLossStates.LOSS,
 923                                            outfile_name + surfacing.PARSIMONY_OUTPUT_FITCH_LOSSES_HTML_BC,
 924                                            sep,
 925                                            ForesterUtil.LINE_SEPARATOR,
 926                                            "Fitch Parsimony | Losses | Domain Combinations",
 927                                            "-",
 928                                            null,
 929                                            all_pfams_encountered,
 930                                            all_pfams_lost_as_dom_combinations,
 931                                            "_fitch_losses_dc",
 932                                            tax_code_to_id_map );
 933             //            writeBinaryStatesMatrixToList( domain_id_to_go_ids_map,
 934             //                                           go_id_to_term_map,
 935             //                                           go_namespace_limit,
 936             //                                           true,
 937             //                                           domain_parsimony.getGainLossMatrix(),
 938             //                                           null,
 939             //                                           outfile_name + surfacing.PARSIMONY_OUTPUT_FITCH_PRESENT_HTML_BC,
 940             //                                           sep,
 941             //                                           ForesterUtil.LINE_SEPARATOR,
 942             //                                           "Fitch Parsimony | Present | Domain Combinations",
 943             //                                           "",
 944             //                                           null,
 945             //                                           all_pfams_encountered,
 946             //                                           null,
 947             //                                           "_fitch_present_dc",
 948             //                                           tax_code_to_id_map );
 949             writeAllEncounteredPfamsToFile( domain_id_to_go_ids_map,
 950                                             go_id_to_term_map,
 951                                             outfile_name,
 952                                             all_pfams_encountered );
 953             writePfamsToFile( outfile_name + surfacing.ALL_PFAMS_GAINED_AS_DOMAINS_SUFFIX,
 954                               all_pfams_gained_as_domains );
 955             writePfamsToFile( outfile_name + surfacing.ALL_PFAMS_LOST_AS_DOMAINS_SUFFIX, all_pfams_lost_as_domains );
 956             writePfamsToFile( outfile_name + surfacing.ALL_PFAMS_GAINED_AS_DC_SUFFIX,
 957                               all_pfams_gained_as_dom_combinations );
 958             writePfamsToFile( outfile_name + surfacing.ALL_PFAMS_LOST_AS_DC_SUFFIX,
 959                               all_pfams_lost_as_dom_combinations );
 960             preparePhylogeny( local_phylogeny_l,
 961                               domain_parsimony,
 962                               date_time,
 963                               "Fitch parsimony on binary domain combination presence/absence randomization: "
 964                                       + randomization,
 965                               "fitch_on_binary_domain_combinations_" + outfile_name,
 966                               parameters_str );
 967             SurfacingUtil
 968                     .writePhylogenyToFile( local_phylogeny_l,
 969                                            outfile_name
 970                                                    + surfacing.BINARY_DOMAIN_COMBINATIONS_PARSIMONY_TREE_OUTPUT_SUFFIX_FITCH );
 971             calculateIndependentDomainCombinationGains( local_phylogeny_l,
 972                                                         outfile_name
 973                                                                 + surfacing.INDEPENDENT_DC_GAINS_FITCH_PARS_COUNTS_OUTPUT_SUFFIX,
 974                                                         outfile_name
 975                                                                 + surfacing.INDEPENDENT_DC_GAINS_FITCH_PARS_DC_OUTPUT_SUFFIX,
 976                                                         outfile_name
 977                                                                 + surfacing.INDEPENDENT_DC_GAINS_FITCH_PARS_DC_FOR_GO_MAPPING_OUTPUT_SUFFIX,
 978                                                         outfile_name
 979                                                                 + surfacing.INDEPENDENT_DC_GAINS_FITCH_PARS_DC_FOR_GO_MAPPING_OUTPUT_UNIQUE_SUFFIX,
 980                                                         outfile_name + "_indep_dc_gains_fitch_lca_ranks.txt",
 981                                                         outfile_name + "_indep_dc_gains_fitch_lca_taxonomies.txt",
 982                                                         outfile_name + "_indep_dc_gains_fitch_protein_statistics.txt",
 983                                                         protein_length_stats_by_dc,
 984                                                         domain_number_stats_by_dc,
 985                                                         domain_length_stats_by_domain );
 986         }
 987     }
 988
 989     public static void executeParsimonyAnalysisForSecondaryFeatures( final String outfile_name,
 990                                                                      final DomainParsimonyCalculator secondary_features_parsimony,
 991                                                                      final Phylogeny phylogeny,
 992                                                                      final String parameters_str,
 993                                                                      final Map<Species, MappingResults> mapping_results_map,
 994                                                                      final boolean use_last_in_fitch_parsimony ) {
 995         final String sep = ForesterUtil.LINE_SEPARATOR + "###################" + ForesterUtil.LINE_SEPARATOR;
 996         final String date_time = ForesterUtil.getCurrentDateTime();
 997         System.out.println();
 998         writeToNexus( outfile_name + surfacing.NEXUS_SECONDARY_FEATURES,
 999                       secondary_features_parsimony.createMatrixOfSecondaryFeaturePresenceOrAbsence( null ),
1000                       phylogeny );
1001         Phylogeny local_phylogeny_copy = phylogeny.copy();
1002         secondary_features_parsimony.executeDolloParsimonyOnSecondaryFeatures( mapping_results_map );
1003         SurfacingUtil.writeMatrixToFile( secondary_features_parsimony.getGainLossMatrix(),
1004                                          outfile_name + surfacing.PARSIMONY_OUTPUT_GL_SUFFIX_DOLLO_SECONDARY_FEATURES,
1005                                          Format.FORESTER );
1006         SurfacingUtil.writeMatrixToFile( secondary_features_parsimony.getGainLossCountsMatrix(),
1007                                          outfile_name
1008                                                  + surfacing.PARSIMONY_OUTPUT_GL_COUNTS_SUFFIX_DOLLO_SECONDARY_FEATURES,
1009                                          Format.FORESTER );
1010         SurfacingUtil.writeBinaryStatesMatrixAsListToFile( secondary_features_parsimony.getGainLossMatrix(),
1011                                                            CharacterStateMatrix.GainLossStates.GAIN,
1012                                                            outfile_name
1013                                                                    + surfacing.PARSIMONY_OUTPUT_DOLLO_GAINS_SECONDARY_FEATURES,
1014                                                            sep,
1015                                                            ForesterUtil.LINE_SEPARATOR,
1016                                                            null );
1017         SurfacingUtil.writeBinaryStatesMatrixAsListToFile( secondary_features_parsimony.getGainLossMatrix(),
1018                                                            CharacterStateMatrix.GainLossStates.LOSS,
1019                                                            outfile_name
1020                                                                    + surfacing.PARSIMONY_OUTPUT_DOLLO_LOSSES_SECONDARY_FEATURES,
1021                                                            sep,
1022                                                            ForesterUtil.LINE_SEPARATOR,
1023                                                            null );
1024         SurfacingUtil.writeBinaryStatesMatrixAsListToFile( secondary_features_parsimony.getGainLossMatrix(),
1025                                                            null,
1026                                                            outfile_name
1027                                                                    + surfacing.PARSIMONY_OUTPUT_DOLLO_PRESENT_SECONDARY_FEATURES,
1028                                                            sep,
1029                                                            ForesterUtil.LINE_SEPARATOR,
1030                                                            null );
1031         preparePhylogeny( local_phylogeny_copy,
1032                           secondary_features_parsimony,
1033                           date_time,
1034                           "Dollo parsimony on secondary feature presence/absence",
1035                           "dollo_on_secondary_features_" + outfile_name,
1036                           parameters_str );
1037         SurfacingUtil
1038                 .writePhylogenyToFile( local_phylogeny_copy,
1039                                        outfile_name + surfacing.SECONDARY_FEATURES_PARSIMONY_TREE_OUTPUT_SUFFIX_DOLLO );
1040         // FITCH DOMAIN COMBINATIONS
1041         // -------------------------
1042         local_phylogeny_copy = phylogeny.copy();
1043         final String randomization = "no";
1044         secondary_features_parsimony
1045                 .executeFitchParsimonyOnBinaryDomainCombintionOnSecondaryFeatures( use_last_in_fitch_parsimony );
1046         preparePhylogeny( local_phylogeny_copy,
1047                           secondary_features_parsimony,
1048                           date_time,
1049                           "Fitch parsimony on secondary binary domain combination presence/absence randomization: "
1050                                   + randomization,
1051                           "fitch_on_binary_domain_combinations_" + outfile_name,
1052                           parameters_str );
1053         SurfacingUtil
1054                 .writePhylogenyToFile( local_phylogeny_copy,
1055                                        outfile_name
1056                                                + surfacing.BINARY_DOMAIN_COMBINATIONS_PARSIMONY_TREE_OUTPUT_SUFFIX_FITCH_MAPPED );
1057         calculateIndependentDomainCombinationGains( local_phylogeny_copy,
1058                                                     outfile_name
1059                                                             + surfacing.INDEPENDENT_DC_GAINS_FITCH_PARS_COUNTS_MAPPED_OUTPUT_SUFFIX,
1060                                                     outfile_name
1061                                                             + surfacing.INDEPENDENT_DC_GAINS_FITCH_PARS_DC_MAPPED_OUTPUT_SUFFIX,
1062                                                     outfile_name
1063                                                             + surfacing.INDEPENDENT_DC_GAINS_FITCH_PARS_DC_FOR_GO_MAPPING_MAPPED_OUTPUT_SUFFIX,
1064                                                     outfile_name
1065                                                             + surfacing.INDEPENDENT_DC_GAINS_FITCH_PARS_DC_FOR_GO_MAPPING_MAPPED_OUTPUT_UNIQUE_SUFFIX,
1066                                                     outfile_name + "_MAPPED_indep_dc_gains_fitch_lca_ranks.txt",
1067                                                     outfile_name + "_MAPPED_indep_dc_gains_fitch_lca_taxonomies.txt",
1068                                                     null,
1069                                                     null,
1070                                                     null,
1071                                                     null );
1072     }
1073
1074     public static void executePlusMinusAnalysis( final File output_file,
1075                                                  final List<String> plus_minus_analysis_high_copy_base,
1076                                                  final List<String> plus_minus_analysis_high_copy_target,
1077                                                  final List<String> plus_minus_analysis_low_copy,
1078                                                  final List<GenomeWideCombinableDomains> gwcd_list,
1079                                                  final SortedMap<Species, List<Protein>> protein_lists_per_species,
1080                                                  final Map<String, List<GoId>> domain_id_to_go_ids_map,
1081                                                  final Map<GoId, GoTerm> go_id_to_term_map,
1082                                                  final List<Object> plus_minus_analysis_numbers ) {
1083         final Set<String> all_spec = new HashSet<String>();
1084         for( final GenomeWideCombinableDomains gwcd : gwcd_list ) {
1085             all_spec.add( gwcd.getSpecies().getSpeciesId() );
1086         }
1087         final File html_out_dom = new File( output_file + surfacing.PLUS_MINUS_DOM_SUFFIX_HTML );
1088         final File plain_out_dom = new File( output_file + surfacing.PLUS_MINUS_DOM_SUFFIX );
1089         final File html_out_dc = new File( output_file + surfacing.PLUS_MINUS_DC_SUFFIX_HTML );
1090         final File all_domains_go_ids_out_dom = new File( output_file + surfacing.PLUS_MINUS_ALL_GO_IDS_DOM_SUFFIX );
1091         final File passing_domains_go_ids_out_dom = new File( output_file
1092                 + surfacing.PLUS_MINUS_PASSING_GO_IDS_DOM_SUFFIX );
1093         final File proteins_file_base = new File( output_file + "" );
1094         final int min_diff = ( ( Integer ) plus_minus_analysis_numbers.get( 0 ) ).intValue();
1095         final double factor = ( ( Double ) plus_minus_analysis_numbers.get( 1 ) ).doubleValue();
1096         try {
1097             DomainCountsDifferenceUtil.calculateCopyNumberDifferences( gwcd_list,
1098                                                                        protein_lists_per_species,
1099                                                                        plus_minus_analysis_high_copy_base,
1100                                                                        plus_minus_analysis_high_copy_target,
1101                                                                        plus_minus_analysis_low_copy,
1102                                                                        min_diff,
1103                                                                        factor,
1104                                                                        plain_out_dom,
1105                                                                        html_out_dom,
1106                                                                        html_out_dc,
1107                                                                        domain_id_to_go_ids_map,
1108                                                                        go_id_to_term_map,
1109                                                                        all_domains_go_ids_out_dom,
1110                                                                        passing_domains_go_ids_out_dom,
1111                                                                        proteins_file_base );
1112         }
1113         catch ( final IOException e ) {
1114             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getLocalizedMessage() );
1115         }
1116         ForesterUtil.programMessage( surfacing.PRG_NAME,
1117                                      "Wrote plus minus domain analysis results to \"" + html_out_dom + "\"" );
1118         ForesterUtil.programMessage( surfacing.PRG_NAME,
1119                                      "Wrote plus minus domain analysis results to \"" + plain_out_dom + "\"" );
1120         ForesterUtil.programMessage( surfacing.PRG_NAME,
1121                                      "Wrote plus minus domain analysis results to \"" + html_out_dc + "\"" );
1122         ForesterUtil.programMessage( surfacing.PRG_NAME,
1123                                      "Wrote plus minus domain analysis based passing GO ids to \""
1124                                              + passing_domains_go_ids_out_dom + "\"" );
1125         ForesterUtil.programMessage( surfacing.PRG_NAME,
1126                                      "Wrote plus minus domain analysis based all GO ids to \""
1127                                              + all_domains_go_ids_out_dom + "\"" );
1128     }
1129
1130     public static void extractProteinNames( final List<Protein> proteins,
1131                                             final List<String> query_domain_ids_nc_order,
1132                                             final Writer out,
1133                                             final String separator,
1134                                             final String limit_to_species )
1135             throws IOException {
1136         for( final Protein protein : proteins ) {
1137             if ( ForesterUtil.isEmpty( limit_to_species )
1138                     || protein.getSpecies().getSpeciesId().equalsIgnoreCase( limit_to_species ) ) {
1139                 if ( protein.contains( query_domain_ids_nc_order, true ) ) {
1140                     out.write( protein.getSpecies().getSpeciesId() );
1141                     out.write( separator );
1142                     out.write( protein.getProteinId().getId() );
1143                     out.write( separator );
1144                     out.write( "[" );
1145                     final Set<String> visited_domain_ids = new HashSet<String>();
1146                     boolean first = true;
1147                     for( final Domain domain : protein.getProteinDomains() ) {
1148                         if ( !visited_domain_ids.contains( domain.getDomainId() ) ) {
1149                             visited_domain_ids.add( domain.getDomainId() );
1150                             if ( first ) {
1151                                 first = false;
1152                             }
1153                             else {
1154                                 out.write( " " );
1155                             }
1156                             out.write( domain.getDomainId() );
1157                             out.write( " {" );
1158                             out.write( "" + domain.getTotalCount() );
1159                             out.write( "}" );
1160                         }
1161                     }
1162                     out.write( "]" );
1163                     out.write( separator );
1164                     if ( !( ForesterUtil.isEmpty( protein.getDescription() )
1165                             || protein.getDescription().equals( SurfacingConstants.NONE ) ) ) {
1166                         out.write( protein.getDescription() );
1167                     }
1168                     out.write( separator );
1169                     if ( !( ForesterUtil.isEmpty( protein.getAccession() )
1170                             || protein.getAccession().equals( SurfacingConstants.NONE ) ) ) {
1171                         out.write( protein.getAccession() );
1172                     }
1173                     out.write( SurfacingConstants.NL );
1174                 }
1175             }
1176         }
1177         out.flush();
1178     }
1179
1180     public static void extractProteinNames( final SortedMap<Species, List<Protein>> protein_lists_per_species,
1181                                             final String domain_id,
1182                                             final Writer out,
1183                                             final String separator,
1184                                             final String limit_to_species,
1185                                             final double domain_e_cutoff )
1186             throws IOException {
1187         //System.out.println( "Per domain E-value: " + domain_e_cutoff );
1188         for( final Species species : protein_lists_per_species.keySet() ) {
1189             //System.out.println( species + ":" );
1190             for( final Protein protein : protein_lists_per_species.get( species ) ) {
1191                 if ( ForesterUtil.isEmpty( limit_to_species )
1192                         || protein.getSpecies().getSpeciesId().equalsIgnoreCase( limit_to_species ) ) {
1193                     final List<Domain> domains = protein.getProteinDomains( domain_id );
1194                     if ( domains.size() > 0 ) {
1195                         out.write( protein.getSpecies().getSpeciesId() );
1196                         out.write( separator );
1197                         out.write( protein.getProteinId().getId() );
1198                         out.write( separator );
1199                         out.write( domain_id.toString() );
1200                         out.write( separator );
1201                         int prev_to = -1;
1202                         for( final Domain domain : domains ) {
1203                             if ( ( domain_e_cutoff < 0 ) || ( domain.getPerDomainEvalue() <= domain_e_cutoff ) ) {
1204                                 out.write( "/" );
1205                                 out.write( domain.getFrom() + "-" + domain.getTo() );
1206                                 if ( prev_to >= 0 ) {
1207                                     final int l = domain.getFrom() - prev_to;
1208                                     // System.out.println( l );
1209                                 }
1210                                 prev_to = domain.getTo();
1211                             }
1212                         }
1213                         out.write( "/" );
1214                         out.write( separator );
1215                         final List<Domain> domain_list = new ArrayList<Domain>();
1216                         for( final Domain domain : protein.getProteinDomains() ) {
1217                             if ( ( domain_e_cutoff < 0 ) || ( domain.getPerDomainEvalue() <= domain_e_cutoff ) ) {
1218                                 domain_list.add( domain );
1219                             }
1220                         }
1221                         final Domain domain_ary[] = new Domain[ domain_list.size() ];
1222                         for( int i = 0; i < domain_list.size(); ++i ) {
1223                             domain_ary[ i ] = domain_list.get( i );
1224                         }
1225                         Arrays.sort( domain_ary, new DomainComparator( true ) );
1226                         out.write( "{" );
1227                         boolean first = true;
1228                         for( final Domain domain : domain_ary ) {
1229                             if ( first ) {
1230                                 first = false;
1231                             }
1232                             else {
1233                                 out.write( "," );
1234                             }
1235                             out.write( domain.getDomainId().toString() );
1236                             out.write( ":" + domain.getFrom() + "-" + domain.getTo() );
1237                             out.write( ":" + domain.getPerDomainEvalue() );
1238                         }
1239                         out.write( "}" );
1240                         if ( !( ForesterUtil.isEmpty( protein.getDescription() )
1241                                 || protein.getDescription().equals( SurfacingConstants.NONE ) ) ) {
1242                             out.write( protein.getDescription() );
1243                         }
1244                         out.write( separator );
1245                         if ( !( ForesterUtil.isEmpty( protein.getAccession() )
1246                                 || protein.getAccession().equals( SurfacingConstants.NONE ) ) ) {
1247                             out.write( protein.getAccession() );
1248                         }
1249                         out.write( SurfacingConstants.NL );
1250                     }
1251                 }
1252             }
1253         }
1254         out.flush();
1255     }
1256
1257     public static SortedSet<String> getAllDomainIds( final List<GenomeWideCombinableDomains> gwcd_list ) {
1258         final SortedSet<String> all_domains_ids = new TreeSet<String>();
1259         for( final GenomeWideCombinableDomains gwcd : gwcd_list ) {
1260             final Set<String> all_domains = gwcd.getAllDomainIds();
1261             //    for( final Domain domain : all_domains ) {
1262             all_domains_ids.addAll( all_domains );
1263             //    }
1264         }
1265         return all_domains_ids;
1266     }
1267
1268     public static SortedMap<String, Integer> getDomainCounts( final List<Protein> protein_domain_collections ) {
1269         final SortedMap<String, Integer> map = new TreeMap<String, Integer>();
1270         for( final Protein protein_domain_collection : protein_domain_collections ) {
1271             for( final Object name : protein_domain_collection.getProteinDomains() ) {
1272                 final BasicDomain protein_domain = ( BasicDomain ) name;
1273                 final String id = protein_domain.getDomainId();
1274                 if ( map.containsKey( id ) ) {
1275                     map.put( id, map.get( id ) + 1 );
1276                 }
1277                 else {
1278                     map.put( id, 1 );
1279                 }
1280             }
1281         }
1282         return map;
1283     }
1284
1285     public static int getNumberOfNodesLackingName( final Phylogeny p, final StringBuilder names ) {
1286         final PhylogenyNodeIterator it = p.iteratorPostorder();
1287         int c = 0;
1288         while ( it.hasNext() ) {
1289             final PhylogenyNode n = it.next();
1290             if ( ForesterUtil.isEmpty( n.getName() )
1291                     && ( !n.getNodeData().isHasTaxonomy()
1292                             || ForesterUtil.isEmpty( n.getNodeData().getTaxonomy().getScientificName() ) )
1293                     && ( !n.getNodeData().isHasTaxonomy()
1294                             || ForesterUtil.isEmpty( n.getNodeData().getTaxonomy().getCommonName() ) ) ) {
1295                 if ( n.getParent() != null ) {
1296                     names.append( " " );
1297                     names.append( n.getParent().getName() );
1298                 }
1299                 final List l = n.getAllExternalDescendants();
1300                 for( final Object object : l ) {
1301                     System.out.println( l.toString() );
1302                 }
1303                 ++c;
1304             }
1305         }
1306         return c;
1307     }
1308
1309     public static void log( final String msg, final Writer w ) {
1310         try {
1311             w.write( msg );
1312             w.write( ForesterUtil.LINE_SEPARATOR );
1313         }
1314         catch ( final IOException e ) {
1315             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getLocalizedMessage() );
1316         }
1317     }
1318
1319     public static Phylogeny[] obtainAndPreProcessIntrees( final File[] intree_files,
1320                                                           final int number_of_genomes,
1321                                                           final String[][] input_file_properties ) {
1322         final Phylogeny[] intrees = new Phylogeny[ intree_files.length ];
1323         int i = 0;
1324         for( final File intree_file : intree_files ) {
1325             Phylogeny intree = null;
1326             final String error = ForesterUtil.isReadableFile( intree_file );
1327             if ( !ForesterUtil.isEmpty( error ) ) {
1328                 ForesterUtil.fatalError( surfacing.PRG_NAME,
1329                                          "cannot read input tree file [" + intree_file + "]: " + error );
1330             }
1331             try {
1332                 final Phylogeny[] p_array = ParserBasedPhylogenyFactory.getInstance()
1333                         .create( intree_file, ParserUtils.createParserDependingOnFileType( intree_file, true ) );
1334                 if ( p_array.length < 1 ) {
1335                     ForesterUtil.fatalError( surfacing.PRG_NAME,
1336                                              "file [" + intree_file
1337                                                      + "] does not contain any phylogeny in phyloXML format" );
1338                 }
1339                 else if ( p_array.length > 1 ) {
1340                     ForesterUtil.fatalError( surfacing.PRG_NAME,
1341                                              "file [" + intree_file
1342                                                      + "] contains more than one phylogeny in phyloXML format" );
1343                 }
1344                 intree = p_array[ 0 ];
1345             }
1346             catch ( final Exception e ) {
1347                 ForesterUtil.fatalError( surfacing.PRG_NAME,
1348                                          "failed to read input tree from file [" + intree_file + "]: " + error );
1349             }
1350             if ( ( intree == null ) || intree.isEmpty() ) {
1351                 ForesterUtil.fatalError( surfacing.PRG_NAME, "input tree [" + intree_file + "] is empty" );
1352             }
1353             if ( !intree.isRooted() ) {
1354                 ForesterUtil.fatalError( surfacing.PRG_NAME, "input tree [" + intree_file + "] is not rooted" );
1355             }
1356             final StringBuilder parent_names = new StringBuilder();
1357             final int nodes_lacking_name = getNumberOfNodesLackingName( intree, parent_names );
1358             if ( nodes_lacking_name > 0 ) {
1359                 ForesterUtil.fatalError( surfacing.PRG_NAME,
1360                                          "input tree [" + intree_file + "] has " + nodes_lacking_name
1361                                                  + " node(s) lacking a name [parent names:" + parent_names + "]" );
1362             }
1363             preparePhylogenyForParsimonyAnalyses( intree, input_file_properties );
1364             if ( !intree.isCompletelyBinary() ) {
1365                 ForesterUtil.printWarningMessage( surfacing.PRG_NAME,
1366                                                   "input tree [" + intree_file + "] is not completely binary" );
1367             }
1368             intrees[ i++ ] = intree;
1369         }
1370         return intrees;
1371     }
1372
1373     public static Phylogeny obtainFirstIntree( final File intree_file ) {
1374         Phylogeny intree = null;
1375         final String error = ForesterUtil.isReadableFile( intree_file );
1376         if ( !ForesterUtil.isEmpty( error ) ) {
1377             ForesterUtil.fatalError( surfacing.PRG_NAME,
1378                                      "cannot read input tree file [" + intree_file + "]: " + error );
1379         }
1380         try {
1381             final Phylogeny[] phys = ParserBasedPhylogenyFactory.getInstance()
1382                     .create( intree_file, ParserUtils.createParserDependingOnFileType( intree_file, true ) );
1383             if ( phys.length < 1 ) {
1384                 ForesterUtil
1385                         .fatalError( surfacing.PRG_NAME,
1386                                      "file [" + intree_file + "] does not contain any phylogeny in phyloXML format" );
1387             }
1388             else if ( phys.length > 1 ) {
1389                 ForesterUtil
1390                         .fatalError( surfacing.PRG_NAME,
1391                                      "file [" + intree_file + "] contains more than one phylogeny in phyloXML format" );
1392             }
1393             intree = phys[ 0 ];
1394         }
1395         catch ( final Exception e ) {
1396             ForesterUtil.fatalError( surfacing.PRG_NAME,
1397                                      "failed to read input tree from file [" + intree_file + "]: " + error );
1398         }
1399         if ( ( intree == null ) || intree.isEmpty() ) {
1400             ForesterUtil.fatalError( surfacing.PRG_NAME, "input tree [" + intree_file + "] is empty" );
1401         }
1402         if ( !intree.isRooted() ) {
1403             ForesterUtil.fatalError( surfacing.PRG_NAME, "input tree [" + intree_file + "] is not rooted" );
1404         }
1405         return intree;
1406     }
1407
1408     public static String obtainHexColorStringDependingOnTaxonomyGroup( final String tax_code, final Phylogeny phy )
1409             throws IllegalArgumentException {
1410         if ( !_TAXCODE_HEXCOLORSTRING_MAP.containsKey( tax_code ) ) {
1411             if ( ( phy != null ) && !phy.isEmpty() ) {
1412                 //                final List<PhylogenyNode> nodes = phy.getNodesViaTaxonomyCode( tax_code );
1413                 //                Color c = null;
1414                 //                if ( ( nodes == null ) || nodes.isEmpty() ) {
1415                 //                    throw new IllegalArgumentException( "code " + tax_code + " is not found" );
1416                 //                }
1417                 //                if ( nodes.size() != 1 ) {
1418                 //                    throw new IllegalArgumentException( "code " + tax_code + " is not unique" );
1419                 //                }
1420                 //                PhylogenyNode n = nodes.get( 0 );
1421                 //                while ( n != null ) {
1422                 //                    if ( n.getNodeData().isHasTaxonomy()
1423                 //                            && !ForesterUtil.isEmpty( n.getNodeData().getTaxonomy().getScientificName() ) ) {
1424                 //                        c = ForesterUtil.obtainColorDependingOnTaxonomyGroup( n.getNodeData().getTaxonomy()
1425                 //                                .getScientificName(), tax_code );
1426                 //                    }
1427                 //                    if ( ( c == null ) && !ForesterUtil.isEmpty( n.getName() ) ) {
1428                 //                        c = ForesterUtil.obtainColorDependingOnTaxonomyGroup( n.getName(), tax_code );
1429                 //                    }
1430                 //                    if ( c != null ) {
1431                 //                        break;
1432                 //                    }
1433                 //                    n = n.getParent();
1434                 //                }
1435                 final String group = obtainTaxonomyGroup( tax_code, phy );
1436                 final Color c = ForesterUtil.obtainColorDependingOnTaxonomyGroup( group );
1437                 if ( c == null ) {
1438                     throw new IllegalArgumentException( "no color found for taxonomy group \"" + group
1439                             + "\" for code \"" + tax_code + "\"" );
1440                 }
1441                 final String hex = String.format( "#%02x%02x%02x", c.getRed(), c.getGreen(), c.getBlue() );
1442                 _TAXCODE_HEXCOLORSTRING_MAP.put( tax_code, hex );
1443             }
1444             else {
1445                 throw new IllegalArgumentException( "unable to obtain color for code " + tax_code
1446                         + " (tree is null or empty and code is not in map)" );
1447             }
1448         }
1449         return _TAXCODE_HEXCOLORSTRING_MAP.get( tax_code );
1450     }
1451
1452     public static String obtainTaxonomyGroup( final String tax_code, final Phylogeny species_tree )
1453             throws IllegalArgumentException {
1454         if ( !_TAXCODE_TAXGROUP_MAP.containsKey( tax_code ) ) {
1455             if ( ( species_tree != null ) && !species_tree.isEmpty() ) {
1456                 final List<PhylogenyNode> nodes = species_tree.getNodesViaTaxonomyCode( tax_code );
1457                 if ( ( nodes == null ) || nodes.isEmpty() ) {
1458                     throw new IllegalArgumentException( "code " + tax_code + " is not found" );
1459                 }
1460                 if ( nodes.size() != 1 ) {
1461                     throw new IllegalArgumentException( "code " + tax_code + " is not unique" );
1462                 }
1463                 PhylogenyNode n = nodes.get( 0 );
1464                 String group = null;
1465                 while ( n != null ) {
1466                     if ( n.getNodeData().isHasTaxonomy()
1467                             && !ForesterUtil.isEmpty( n.getNodeData().getTaxonomy().getScientificName() ) ) {
1468                         group = ForesterUtil
1469                                 .obtainNormalizedTaxonomyGroup( n.getNodeData().getTaxonomy().getScientificName() );
1470                     }
1471                     if ( ForesterUtil.isEmpty( group ) && !ForesterUtil.isEmpty( n.getName() ) ) {
1472                         group = ForesterUtil.obtainNormalizedTaxonomyGroup( n.getName() );
1473                     }
1474                     if ( !ForesterUtil.isEmpty( group ) ) {
1475                         break;
1476                     }
1477                     n = n.getParent();
1478                 }
1479                 if ( ForesterUtil.isEmpty( group ) ) {
1480                     group = TaxonomyGroups.OTHER;
1481                 }
1482                 _TAXCODE_TAXGROUP_MAP.put( tax_code, group );
1483             }
1484             else {
1485                 throw new IllegalArgumentException( "unable to obtain group for code " + tax_code
1486                         + " (tree is null or empty and code is not in map)" );
1487             }
1488         }
1489         return _TAXCODE_TAXGROUP_MAP.get( tax_code );
1490     }
1491
1492     public static void performDomainArchitectureAnalysis( final SortedMap<String, Set<String>> domain_architecutures,
1493                                                           final SortedMap<String, Integer> domain_architecuture_counts,
1494                                                           final int min_count,
1495                                                           final File da_counts_outfile,
1496                                                           final File unique_da_outfile ) {
1497         checkForOutputFileWriteability( da_counts_outfile );
1498         checkForOutputFileWriteability( unique_da_outfile );
1499         try {
1500             final BufferedWriter da_counts_out = new BufferedWriter( new FileWriter( da_counts_outfile ) );
1501             final BufferedWriter unique_da_out = new BufferedWriter( new FileWriter( unique_da_outfile ) );
1502             final Iterator<Entry<String, Integer>> it = domain_architecuture_counts.entrySet().iterator();
1503             while ( it.hasNext() ) {
1504                 final Map.Entry<String, Integer> e = it.next();
1505                 final String da = e.getKey();
1506                 final int count = e.getValue();
1507                 if ( count >= min_count ) {
1508                     da_counts_out.write( da );
1509                     da_counts_out.write( "\t" );
1510                     da_counts_out.write( String.valueOf( count ) );
1511                     da_counts_out.write( ForesterUtil.LINE_SEPARATOR );
1512                 }
1513                 if ( count == 1 ) {
1514                     final Iterator<Entry<String, Set<String>>> it2 = domain_architecutures.entrySet().iterator();
1515                     while ( it2.hasNext() ) {
1516                         final Map.Entry<String, Set<String>> e2 = it2.next();
1517                         final String genome = e2.getKey();
1518                         final Set<String> das = e2.getValue();
1519                         if ( das.contains( da ) ) {
1520                             unique_da_out.write( genome );
1521                             unique_da_out.write( "\t" );
1522                             unique_da_out.write( da );
1523                             unique_da_out.write( ForesterUtil.LINE_SEPARATOR );
1524                         }
1525                     }
1526                 }
1527             }
1528             unique_da_out.close();
1529             da_counts_out.close();
1530         }
1531         catch ( final IOException e ) {
1532             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getMessage() );
1533         }
1534         ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote distance matrices to \"" + da_counts_outfile + "\"" );
1535         ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote distance matrices to \"" + unique_da_outfile + "\"" );
1536         //
1537     }
1538
1539     public static void preparePhylogeny( final Phylogeny p,
1540                                          final DomainParsimonyCalculator domain_parsimony,
1541                                          final String date_time,
1542                                          final String method,
1543                                          final String name,
1544                                          final String parameters_str ) {
1545         domain_parsimony.decoratePhylogenyWithDomains( p );
1546         final StringBuilder desc = new StringBuilder();
1547         desc.append( "[Method: " + method + "] [Date: " + date_time + "] " );
1548         desc.append( "[Cost: " + domain_parsimony.getCost() + "] " );
1549         desc.append( "[Gains: " + domain_parsimony.getTotalGains() + "] " );
1550         desc.append( "[Losses: " + domain_parsimony.getTotalLosses() + "] " );
1551         desc.append( "[Unchanged: " + domain_parsimony.getTotalUnchanged() + "] " );
1552         desc.append( "[Parameters: " + parameters_str + "]" );
1553         p.setName( name );
1554         p.setDescription( desc.toString() );
1555         p.setConfidence( new Confidence( domain_parsimony.getCost(), "parsimony" ) );
1556         p.setRerootable( false );
1557         p.setRooted( true );
1558     }
1559
1560     public static void preparePhylogenyForParsimonyAnalyses( final Phylogeny intree,
1561                                                              final String[][] input_file_properties ) {
1562         final String[] genomes = new String[ input_file_properties.length ];
1563         for( int i = 0; i < input_file_properties.length; ++i ) {
1564             if ( intree.getNodes( input_file_properties[ i ][ 1 ] ).size() > 1 ) {
1565                 ForesterUtil.fatalError( surfacing.PRG_NAME,
1566                                          "node named [" + input_file_properties[ i ][ 1 ]
1567                                                  + "] is not unique in input tree " + intree.getName() );
1568             }
1569             genomes[ i ] = input_file_properties[ i ][ 1 ];
1570         }
1571         //
1572         final PhylogenyNodeIterator it = intree.iteratorPostorder();
1573         while ( it.hasNext() ) {
1574             final PhylogenyNode n = it.next();
1575             if ( ForesterUtil.isEmpty( n.getName() ) ) {
1576                 if ( n.getNodeData().isHasTaxonomy()
1577                         && !ForesterUtil.isEmpty( n.getNodeData().getTaxonomy().getTaxonomyCode() ) ) {
1578                     n.setName( n.getNodeData().getTaxonomy().getTaxonomyCode() );
1579                 }
1580                 else if ( n.getNodeData().isHasTaxonomy()
1581                         && !ForesterUtil.isEmpty( n.getNodeData().getTaxonomy().getScientificName() ) ) {
1582                     n.setName( n.getNodeData().getTaxonomy().getScientificName() );
1583                 }
1584                 else if ( n.getNodeData().isHasTaxonomy()
1585                         && !ForesterUtil.isEmpty( n.getNodeData().getTaxonomy().getCommonName() ) ) {
1586                     n.setName( n.getNodeData().getTaxonomy().getCommonName() );
1587                 }
1588                 else {
1589                     ForesterUtil
1590                             .fatalError( surfacing.PRG_NAME,
1591                                          "node with no name, scientific name, common name, or taxonomy code present" );
1592                 }
1593             }
1594         }
1595         final List<String> igns = PhylogenyMethods.deleteExternalNodesPositiveSelection( genomes, intree );
1596         if ( igns.size() > 0 ) {
1597             System.out.println( "Not using the following " + igns.size() + " nodes:" );
1598             for( int i = 0; i < igns.size(); ++i ) {
1599                 System.out.println( " " + i + ": " + igns.get( i ) );
1600             }
1601             System.out.println( "--" );
1602         }
1603         //Test for node names:
1604         final SortedSet<String> not_found = new TreeSet<String>();
1605         final SortedSet<String> not_unique = new TreeSet<String>();
1606         for( final String[] input_file_propertie : input_file_properties ) {
1607             final String name = input_file_propertie[ 1 ];
1608             final List<PhylogenyNode> nodes = intree.getNodes( name );
1609             if ( ( nodes == null ) || ( nodes.size() < 1 ) ) {
1610                 not_found.add( name );
1611             }
1612             if ( nodes.size() > 1 ) {
1613                 not_unique.add( name );
1614             }
1615         }
1616         if ( not_found.size() > 0 ) {
1617             ForesterUtil.fatalError( surfacing.PRG_NAME,
1618                                      "the following " + not_found.size()
1619                                              + " node(s) are not present in the input tree: " + not_found );
1620         }
1621         if ( not_unique.size() > 0 ) {
1622             ForesterUtil.fatalError( surfacing.PRG_NAME,
1623                                      "the following " + not_unique.size()
1624                                              + " node(s) are not unique in the input tree: " + not_unique );
1625         }
1626     }
1627
1628     public static void printOutPercentageOfMultidomainProteins( final SortedMap<Integer, Integer> all_genomes_domains_per_potein_histo,
1629                                                                 final Writer log_writer ) {
1630         int sum = 0;
1631         for( final Entry<Integer, Integer> entry : all_genomes_domains_per_potein_histo.entrySet() ) {
1632             sum += entry.getValue();
1633         }
1634         final double percentage = ( 100.0 * ( sum - all_genomes_domains_per_potein_histo.get( 1 ) ) ) / sum;
1635         ForesterUtil.programMessage( surfacing.PRG_NAME, "Percentage of multidomain proteins: " + percentage + "%" );
1636         log( "Percentage of multidomain proteins:            : " + percentage + "%", log_writer );
1637     }
1638
1639     public static void processFilter( final File filter_file, final SortedSet<String> filter ) {
1640         SortedSet<String> filter_str = null;
1641         try {
1642             filter_str = ForesterUtil.file2set( filter_file );
1643         }
1644         catch ( final IOException e ) {
1645             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getMessage() );
1646         }
1647         if ( filter_str != null ) {
1648             for( final String string : filter_str ) {
1649                 filter.add( string );
1650             }
1651         }
1652         if ( surfacing.VERBOSE ) {
1653             System.out.println( "Filter:" );
1654             for( final String domainId : filter ) {
1655                 System.out.println( domainId );
1656             }
1657         }
1658     }
1659
1660     public static String[][] processInputGenomesFile( final File input_genomes ) {
1661         String[][] input_file_properties = null;
1662         try {
1663             input_file_properties = ForesterUtil.file22dArray( input_genomes );
1664         }
1665         catch ( final IOException e ) {
1666             ForesterUtil.fatalError( surfacing.PRG_NAME,
1667                                      "genomes files is to be in the following format \"<hmmpfam output file> <species>\": "
1668                                              + e.getLocalizedMessage() );
1669         }
1670         final Set<String> specs = new HashSet<String>();
1671         final Set<String> paths = new HashSet<String>();
1672         for( int i = 0; i < input_file_properties.length; ++i ) {
1673             if ( !PhyloXmlUtil.TAXOMONY_CODE_PATTERN.matcher( input_file_properties[ i ][ 1 ] ).matches() ) {
1674                 ForesterUtil.fatalError( surfacing.PRG_NAME,
1675                                          "illegal format for species code: " + input_file_properties[ i ][ 1 ] );
1676             }
1677             if ( specs.contains( input_file_properties[ i ][ 1 ] ) ) {
1678                 ForesterUtil.fatalError( surfacing.PRG_NAME,
1679                                          "species code " + input_file_properties[ i ][ 1 ] + " is not unique" );
1680             }
1681             specs.add( input_file_properties[ i ][ 1 ] );
1682             if ( paths.contains( input_file_properties[ i ][ 0 ] ) ) {
1683                 ForesterUtil.fatalError( surfacing.PRG_NAME,
1684                                          "path " + input_file_properties[ i ][ 0 ] + " is not unique" );
1685             }
1686             paths.add( input_file_properties[ i ][ 0 ] );
1687             final String error = ForesterUtil.isReadableFile( new File( input_file_properties[ i ][ 0 ] ) );
1688             if ( !ForesterUtil.isEmpty( error ) ) {
1689                 ForesterUtil.fatalError( surfacing.PRG_NAME, error );
1690             }
1691         }
1692         return input_file_properties;
1693     }
1694
1695     public static void processPlusMinusAnalysisOption( final CommandLineArguments cla,
1696                                                        final List<String> high_copy_base,
1697                                                        final List<String> high_copy_target,
1698                                                        final List<String> low_copy,
1699                                                        final List<Object> numbers ) {
1700         if ( cla.isOptionSet( surfacing.PLUS_MINUS_ANALYSIS_OPTION ) ) {
1701             if ( !cla.isOptionValueSet( surfacing.PLUS_MINUS_ANALYSIS_OPTION ) ) {
1702                 ForesterUtil.fatalError( surfacing.PRG_NAME,
1703                                          "no value for 'plus-minus' file: -" + surfacing.PLUS_MINUS_ANALYSIS_OPTION
1704                                                  + "=<file>" );
1705             }
1706             final File plus_minus_file = new File( cla.getOptionValue( surfacing.PLUS_MINUS_ANALYSIS_OPTION ) );
1707             final String msg = ForesterUtil.isReadableFile( plus_minus_file );
1708             if ( !ForesterUtil.isEmpty( msg ) ) {
1709                 ForesterUtil.fatalError( surfacing.PRG_NAME, "can not read from \"" + plus_minus_file + "\": " + msg );
1710             }
1711             processPlusMinusFile( plus_minus_file, high_copy_base, high_copy_target, low_copy, numbers );
1712         }
1713     }
1714
1715     // First numbers is minimal difference, second is factor.
1716     public static void processPlusMinusFile( final File plus_minus_file,
1717                                              final List<String> high_copy_base,
1718                                              final List<String> high_copy_target,
1719                                              final List<String> low_copy,
1720                                              final List<Object> numbers ) {
1721         Set<String> species_set = null;
1722         int min_diff = surfacing.PLUS_MINUS_ANALYSIS_MIN_DIFF_DEFAULT;
1723         double factor = surfacing.PLUS_MINUS_ANALYSIS_FACTOR_DEFAULT;
1724         try {
1725             species_set = ForesterUtil.file2set( plus_minus_file );
1726         }
1727         catch ( final IOException e ) {
1728             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getMessage() );
1729         }
1730         if ( species_set != null ) {
1731             for( final String species : species_set ) {
1732                 final String species_trimmed = species.substring( 1 );
1733                 if ( species.startsWith( "+" ) ) {
1734                     if ( low_copy.contains( species_trimmed ) ) {
1735                         ForesterUtil.fatalError( surfacing.PRG_NAME,
1736                                                  "species/genome names can not appear with both '+' and '-' suffix, as appears the case for: \""
1737                                                          + species_trimmed + "\"" );
1738                     }
1739                     high_copy_base.add( species_trimmed );
1740                 }
1741                 else if ( species.startsWith( "*" ) ) {
1742                     if ( low_copy.contains( species_trimmed ) ) {
1743                         ForesterUtil.fatalError( surfacing.PRG_NAME,
1744                                                  "species/genome names can not appear with both '*' and '-' suffix, as appears the case for: \""
1745                                                          + species_trimmed + "\"" );
1746                     }
1747                     high_copy_target.add( species_trimmed );
1748                 }
1749                 else if ( species.startsWith( "-" ) ) {
1750                     if ( high_copy_base.contains( species_trimmed ) || high_copy_target.contains( species_trimmed ) ) {
1751                         ForesterUtil.fatalError( surfacing.PRG_NAME,
1752                                                  "species/genome names can not appear with both '+' or '*' and '-' suffix, as appears the case for: \""
1753                                                          + species_trimmed + "\"" );
1754                     }
1755                     low_copy.add( species_trimmed );
1756                 }
1757                 else if ( species.startsWith( "$D" ) ) {
1758                     try {
1759                         min_diff = Integer.parseInt( species.substring( 3 ) );
1760                     }
1761                     catch ( final NumberFormatException e ) {
1762                         ForesterUtil.fatalError( surfacing.PRG_NAME,
1763                                                  "could not parse integer value for minimal difference from: \""
1764                                                          + species.substring( 3 ) + "\"" );
1765                     }
1766                 }
1767                 else if ( species.startsWith( "$F" ) ) {
1768                     try {
1769                         factor = Double.parseDouble( species.substring( 3 ) );
1770                     }
1771                     catch ( final NumberFormatException e ) {
1772                         ForesterUtil.fatalError( surfacing.PRG_NAME,
1773                                                  "could not parse double value for factor from: \""
1774                                                          + species.substring( 3 ) + "\"" );
1775                     }
1776                 }
1777                 else if ( species.startsWith( "#" ) ) {
1778                     // Comment, ignore.
1779                 }
1780                 else {
1781                     ForesterUtil.fatalError( surfacing.PRG_NAME,
1782                                              "species/genome names in 'plus minus' file must begin with '*' (high copy target genome), '+' (high copy base genomes), '-' (low copy genomes), '$D=<integer>' minimal Difference (default is 1), '$F=<double>' factor (default is 1.0), double), or '#' (ignore) suffix, encountered: \""
1783                                                      + species + "\"" );
1784                 }
1785                 numbers.add( new Integer( min_diff + "" ) );
1786                 numbers.add( new Double( factor + "" ) );
1787             }
1788         }
1789         else {
1790             ForesterUtil.fatalError( surfacing.PRG_NAME, "'plus minus' file [" + plus_minus_file + "] appears empty" );
1791         }
1792     }
1793
1794     /*
1795      * species | protein id | n-terminal domain | c-terminal domain | n-terminal domain per domain E-value | c-terminal domain per domain E-value
1796      *
1797      *
1798      */
1799     static public StringBuffer proteinToDomainCombinations( final Protein protein,
1800                                                             final String protein_id,
1801                                                             final String separator ) {
1802         final StringBuffer sb = new StringBuffer();
1803         if ( protein.getSpecies() == null ) {
1804             throw new IllegalArgumentException( "species must not be null" );
1805         }
1806         if ( ForesterUtil.isEmpty( protein.getSpecies().getSpeciesId() ) ) {
1807             throw new IllegalArgumentException( "species id must not be empty" );
1808         }
1809         final List<Domain> domains = protein.getProteinDomains();
1810         if ( domains.size() > 1 ) {
1811             final Map<String, Integer> counts = new HashMap<String, Integer>();
1812             for( final Domain domain : domains ) {
1813                 final String id = domain.getDomainId();
1814                 if ( counts.containsKey( id ) ) {
1815                     counts.put( id, counts.get( id ) + 1 );
1816                 }
1817                 else {
1818                     counts.put( id, 1 );
1819                 }
1820             }
1821             final Set<String> dcs = new HashSet<String>();
1822             for( int i = 1; i < domains.size(); ++i ) {
1823                 for( int j = 0; j < i; ++j ) {
1824                     Domain domain_n = domains.get( i );
1825                     Domain domain_c = domains.get( j );
1826                     if ( domain_n.getFrom() > domain_c.getFrom() ) {
1827                         domain_n = domains.get( j );
1828                         domain_c = domains.get( i );
1829                     }
1830                     final String dc = domain_n.getDomainId() + domain_c.getDomainId();
1831                     if ( !dcs.contains( dc ) ) {
1832                         dcs.add( dc );
1833                         sb.append( protein.getSpecies() );
1834                         sb.append( separator );
1835                         sb.append( protein_id );
1836                         sb.append( separator );
1837                         sb.append( domain_n.getDomainId() );
1838                         sb.append( separator );
1839                         sb.append( domain_c.getDomainId() );
1840                         sb.append( separator );
1841                         sb.append( domain_n.getPerDomainEvalue() );
1842                         sb.append( separator );
1843                         sb.append( domain_c.getPerDomainEvalue() );
1844                         sb.append( separator );
1845                         sb.append( counts.get( domain_n.getDomainId() ) );
1846                         sb.append( separator );
1847                         sb.append( counts.get( domain_c.getDomainId() ) );
1848                         sb.append( ForesterUtil.LINE_SEPARATOR );
1849                     }
1850                 }
1851             }
1852         }
1853         else if ( domains.size() == 1 ) {
1854             sb.append( protein.getSpecies() );
1855             sb.append( separator );
1856             sb.append( protein_id );
1857             sb.append( separator );
1858             sb.append( domains.get( 0 ).getDomainId() );
1859             sb.append( separator );
1860             sb.append( separator );
1861             sb.append( domains.get( 0 ).getPerDomainEvalue() );
1862             sb.append( separator );
1863             sb.append( separator );
1864             sb.append( 1 );
1865             sb.append( separator );
1866             sb.append( ForesterUtil.LINE_SEPARATOR );
1867         }
1868         else {
1869             sb.append( protein.getSpecies() );
1870             sb.append( separator );
1871             sb.append( protein_id );
1872             sb.append( separator );
1873             sb.append( separator );
1874             sb.append( separator );
1875             sb.append( separator );
1876             sb.append( separator );
1877             sb.append( separator );
1878             sb.append( ForesterUtil.LINE_SEPARATOR );
1879         }
1880         return sb;
1881     }
1882
1883     public static List<Domain> sortDomainsWithAscendingConfidenceValues( final Protein protein ) {
1884         final List<Domain> domains = new ArrayList<Domain>();
1885         for( final Domain d : protein.getProteinDomains() ) {
1886             domains.add( d );
1887         }
1888         Collections.sort( domains, SurfacingUtil.ASCENDING_CONFIDENCE_VALUE_ORDER );
1889         return domains;
1890     }
1891
1892     public static int storeDomainArchitectures( final String genome,
1893                                                 final SortedMap<String, Set<String>> domain_architecutures,
1894                                                 final List<Protein> protein_list,
1895                                                 final Map<String, Integer> distinct_domain_architecuture_counts ) {
1896         final Set<String> da = new HashSet<String>();
1897         domain_architecutures.put( genome, da );
1898         for( final Protein protein : protein_list ) {
1899             final String da_str = ( ( BasicProtein ) protein ).toDomainArchitectureString( "~", 3, "=" );
1900             if ( !da.contains( da_str ) ) {
1901                 if ( !distinct_domain_architecuture_counts.containsKey( da_str ) ) {
1902                     distinct_domain_architecuture_counts.put( da_str, 1 );
1903                 }
1904                 else {
1905                     distinct_domain_architecuture_counts.put( da_str,
1906                                                               distinct_domain_architecuture_counts.get( da_str ) + 1 );
1907                 }
1908                 da.add( da_str );
1909             }
1910         }
1911         return da.size();
1912     }
1913
1914     public static void writeAllDomainsChangedOnAllSubtrees( final Phylogeny p,
1915                                                             final boolean get_gains,
1916                                                             final String outdir,
1917                                                             final String suffix_for_filename )
1918             throws IOException {
1919         CharacterStateMatrix.GainLossStates state = CharacterStateMatrix.GainLossStates.GAIN;
1920         if ( !get_gains ) {
1921             state = CharacterStateMatrix.GainLossStates.LOSS;
1922         }
1923         final File base_dir = createBaseDirForPerNodeDomainFiles( surfacing.BASE_DIRECTORY_PER_SUBTREE_DOMAIN_GAIN_LOSS_FILES,
1924                                                                   false,
1925                                                                   state,
1926                                                                   outdir );
1927         for( final PhylogenyNodeIterator it = p.iteratorPostorder(); it.hasNext(); ) {
1928             final PhylogenyNode node = it.next();
1929             if ( !node.isExternal() ) {
1930                 final SortedSet<String> domains = collectAllDomainsChangedOnSubtree( node, get_gains );
1931                 if ( domains.size() > 0 ) {
1932                     final Writer writer = ForesterUtil.createBufferedWriter( base_dir + ForesterUtil.FILE_SEPARATOR
1933                             + node.getName() + suffix_for_filename );
1934                     for( final String domain : domains ) {
1935                         writer.write( domain );
1936                         writer.write( ForesterUtil.LINE_SEPARATOR );
1937                     }
1938                     writer.close();
1939                 }
1940             }
1941         }
1942     }
1943
1944     public static void writeBinaryDomainCombinationsFileForGraphAnalysis( final String[][] input_file_properties,
1945                                                                           final File output_dir,
1946                                                                           final GenomeWideCombinableDomains gwcd,
1947                                                                           final int i,
1948                                                                           final GenomeWideCombinableDomainsSortOrder dc_sort_order ) {
1949         File dc_outfile_dot = new File( input_file_properties[ i ][ 1 ]
1950                 + surfacing.DOMAIN_COMBINITONS_OUTPUTFILE_SUFFIX_FOR_GRAPH_ANALYSIS );
1951         if ( output_dir != null ) {
1952             dc_outfile_dot = new File( output_dir + ForesterUtil.FILE_SEPARATOR + dc_outfile_dot );
1953         }
1954         checkForOutputFileWriteability( dc_outfile_dot );
1955         final SortedSet<BinaryDomainCombination> binary_combinations = createSetOfAllBinaryDomainCombinationsPerGenome( gwcd );
1956         try {
1957             final BufferedWriter out_dot = new BufferedWriter( new FileWriter( dc_outfile_dot ) );
1958             for( final BinaryDomainCombination bdc : binary_combinations ) {
1959                 out_dot.write( bdc.toGraphDescribingLanguage( BinaryDomainCombination.OutputFormat.DOT, null, null )
1960                         .toString() );
1961                 out_dot.write( SurfacingConstants.NL );
1962             }
1963             out_dot.close();
1964         }
1965         catch ( final IOException e ) {
1966             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getMessage() );
1967         }
1968         if ( input_file_properties[ i ].length == 3 ) {
1969             ForesterUtil
1970                     .programMessage( surfacing.PRG_NAME,
1971                                      "Wrote binary domain combination for \"" + input_file_properties[ i ][ 0 ] + "\" ("
1972                                              + input_file_properties[ i ][ 1 ] + ", " + input_file_properties[ i ][ 2 ]
1973                                              + ") to: \"" + dc_outfile_dot + "\"" );
1974         }
1975         else {
1976             ForesterUtil.programMessage( surfacing.PRG_NAME,
1977                                          "Wrote binary domain combination for \"" + input_file_properties[ i ][ 0 ]
1978                                                  + "\" (" + input_file_properties[ i ][ 1 ] + ") to: \""
1979                                                  + dc_outfile_dot + "\"" );
1980         }
1981     }
1982
1983     public static void writeBinaryStatesMatrixAsListToFile( final CharacterStateMatrix<CharacterStateMatrix.GainLossStates> matrix,
1984                                                             final CharacterStateMatrix.GainLossStates state,
1985                                                             final String filename,
1986                                                             final String indentifier_characters_separator,
1987                                                             final String character_separator,
1988                                                             final Map<String, String> descriptions ) {
1989         final File outfile = new File( filename );
1990         checkForOutputFileWriteability( outfile );
1991         final SortedSet<String> sorted_ids = new TreeSet<String>();
1992         for( int i = 0; i < matrix.getNumberOfIdentifiers(); ++i ) {
1993             sorted_ids.add( matrix.getIdentifier( i ) );
1994         }
1995         try {
1996             final BufferedWriter out = new BufferedWriter( new FileWriter( outfile ) );
1997             for( final String id : sorted_ids ) {
1998                 out.write( indentifier_characters_separator );
1999                 out.write( "#" + id );
2000                 out.write( indentifier_characters_separator );
2001                 for( int c = 0; c < matrix.getNumberOfCharacters(); ++c ) {
2002                     // Not nice:
2003                     // using null to indicate either UNCHANGED_PRESENT or GAIN.
2004                     if ( ( matrix.getState( id, c ) == state ) || ( ( state == null )
2005                             && ( ( matrix.getState( id, c ) == CharacterStateMatrix.GainLossStates.GAIN ) || ( matrix
2006                                     .getState( id,
2007                                                c ) == CharacterStateMatrix.GainLossStates.UNCHANGED_PRESENT ) ) ) ) {
2008                         out.write( matrix.getCharacter( c ) );
2009                         if ( ( descriptions != null ) && !descriptions.isEmpty()
2010                                 && descriptions.containsKey( matrix.getCharacter( c ) ) ) {
2011                             out.write( "\t" );
2012                             out.write( descriptions.get( matrix.getCharacter( c ) ) );
2013                         }
2014                         out.write( character_separator );
2015                     }
2016                 }
2017             }
2018             out.flush();
2019             out.close();
2020         }
2021         catch ( final IOException e ) {
2022             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getMessage() );
2023         }
2024         ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote characters list: \"" + filename + "\"" );
2025     }
2026
2027     public static void writeBinaryStatesMatrixAsListToFileForBinaryCombinationsForGraphAnalysis( final CharacterStateMatrix<CharacterStateMatrix.GainLossStates> matrix,
2028                                                                                                  final CharacterStateMatrix.GainLossStates state,
2029                                                                                                  final String filename,
2030                                                                                                  final String indentifier_characters_separator,
2031                                                                                                  final String character_separator,
2032                                                                                                  final BinaryDomainCombination.OutputFormat bc_output_format ) {
2033         final File outfile = new File( filename );
2034         checkForOutputFileWriteability( outfile );
2035         final SortedSet<String> sorted_ids = new TreeSet<String>();
2036         for( int i = 0; i < matrix.getNumberOfIdentifiers(); ++i ) {
2037             sorted_ids.add( matrix.getIdentifier( i ) );
2038         }
2039         try {
2040             final BufferedWriter out = new BufferedWriter( new FileWriter( outfile ) );
2041             for( final String id : sorted_ids ) {
2042                 out.write( indentifier_characters_separator );
2043                 out.write( "#" + id );
2044                 out.write( indentifier_characters_separator );
2045                 for( int c = 0; c < matrix.getNumberOfCharacters(); ++c ) {
2046                     // Not nice:
2047                     // using null to indicate either UNCHANGED_PRESENT or GAIN.
2048                     if ( ( matrix.getState( id, c ) == state ) || ( ( state == null )
2049                             && ( ( matrix.getState( id, c ) == CharacterStateMatrix.GainLossStates.GAIN ) || ( matrix
2050                                     .getState( id,
2051                                                c ) == CharacterStateMatrix.GainLossStates.UNCHANGED_PRESENT ) ) ) ) {
2052                         BinaryDomainCombination bdc = null;
2053                         try {
2054                             bdc = BasicBinaryDomainCombination.obtainInstance( matrix.getCharacter( c ) );
2055                         }
2056                         catch ( final Exception e ) {
2057                             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getLocalizedMessage() );
2058                         }
2059                         out.write( bdc.toGraphDescribingLanguage( bc_output_format, null, null ).toString() );
2060                         out.write( character_separator );
2061                     }
2062                 }
2063             }
2064             out.flush();
2065             out.close();
2066         }
2067         catch ( final IOException e ) {
2068             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getMessage() );
2069         }
2070         ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote characters list: \"" + filename + "\"" );
2071     }
2072
2073     public static void writeBinaryStatesMatrixToList( final Map<String, List<GoId>> domain_id_to_go_ids_map,
2074                                                       final Map<GoId, GoTerm> go_id_to_term_map,
2075                                                       final GoNameSpace go_namespace_limit,
2076                                                       final boolean domain_combinations,
2077                                                       final CharacterStateMatrix<CharacterStateMatrix.GainLossStates> matrix,
2078                                                       final CharacterStateMatrix.GainLossStates state,
2079                                                       final String filename,
2080                                                       final String indentifier_characters_separator,
2081                                                       final String character_separator,
2082                                                       final String title_for_html,
2083                                                       final String prefix_for_html,
2084                                                       final Map<String, Set<String>>[] domain_id_to_secondary_features_maps,
2085                                                       final SortedSet<String> all_pfams_encountered,
2086                                                       final SortedSet<String> pfams_gained_or_lost,
2087                                                       final String suffix_for_per_node_events_file,
2088                                                       final Map<String, Integer> tax_code_to_id_map ) {
2089         if ( ( go_namespace_limit != null ) && ( ( go_id_to_term_map == null ) || ( go_id_to_term_map.size() < 1 ) ) ) {
2090             throw new IllegalArgumentException( "attempt to use GO namespace limit without a GO-id to term map" );
2091         }
2092         else if ( ( ( domain_id_to_go_ids_map == null ) || ( domain_id_to_go_ids_map.size() < 1 ) ) ) {
2093             throw new IllegalArgumentException( "attempt to output detailed HTML without a Pfam to GO map" );
2094         }
2095         else if ( ( ( go_id_to_term_map == null ) || ( go_id_to_term_map.size() < 1 ) ) ) {
2096             throw new IllegalArgumentException( "attempt to output detailed HTML without a GO-id to term map" );
2097         }
2098         final File outfile = new File( filename );
2099         checkForOutputFileWriteability( outfile );
2100         final SortedSet<String> sorted_ids = new TreeSet<String>();
2101         for( int i = 0; i < matrix.getNumberOfIdentifiers(); ++i ) {
2102             sorted_ids.add( matrix.getIdentifier( i ) );
2103         }
2104         try {
2105             final Writer out = new BufferedWriter( new FileWriter( outfile ) );
2106             final File per_node_go_mapped_domain_gain_loss_files_base_dir = createBaseDirForPerNodeDomainFiles( surfacing.BASE_DIRECTORY_PER_NODE_DOMAIN_GAIN_LOSS_FILES,
2107                                                                                                                 domain_combinations,
2108                                                                                                                 state,
2109                                                                                                                 filename );
2110             Writer per_node_go_mapped_domain_gain_loss_outfile_writer = null;
2111             File per_node_go_mapped_domain_gain_loss_outfile = null;
2112             int per_node_counter = 0;
2113             out.write( "<html>" );
2114             out.write( SurfacingConstants.NL );
2115             writeHtmlHead( out, title_for_html );
2116             out.write( SurfacingConstants.NL );
2117             out.write( "<body>" );
2118             out.write( SurfacingConstants.NL );
2119             out.write( "<h1>" );
2120             out.write( SurfacingConstants.NL );
2121             out.write( title_for_html );
2122             out.write( SurfacingConstants.NL );
2123             out.write( "</h1>" );
2124             out.write( SurfacingConstants.NL );
2125             out.write( "<table>" );
2126             out.write( SurfacingConstants.NL );
2127             for( final String id : sorted_ids ) {
2128                 final Matcher matcher = PATTERN_SP_STYLE_TAXONOMY.matcher( id );
2129                 if ( matcher.matches() ) {
2130                     continue;
2131                 }
2132                 out.write( "<tr>" );
2133                 out.write( "<td>" );
2134                 out.write( "<a href=\"#" + id + "\">" + id + "</a>" );
2135                 out.write( "</td>" );
2136                 out.write( "</tr>" );
2137                 out.write( SurfacingConstants.NL );
2138             }
2139             out.write( "</table>" );
2140             out.write( SurfacingConstants.NL );
2141             for( final String id : sorted_ids ) {
2142                 final Matcher matcher = PATTERN_SP_STYLE_TAXONOMY.matcher( id );
2143                 if ( matcher.matches() ) {
2144                     continue;
2145                 }
2146                 out.write( SurfacingConstants.NL );
2147                 out.write( "<h2>" );
2148                 out.write( "<a name=\"" + id + "\">" + id + "</a>" );
2149                 writeTaxonomyLinks( out, id, tax_code_to_id_map );
2150                 out.write( "</h2>" );
2151                 out.write( SurfacingConstants.NL );
2152                 out.write( "<table>" );
2153                 out.write( SurfacingConstants.NL );
2154                 out.write( "<tr>" );
2155                 out.write( "<td><b>" );
2156                 out.write( "Pfam domain(s)" );
2157                 out.write( "</b></td><td><b>" );
2158                 out.write( "GO term acc" );
2159                 out.write( "</b></td><td><b>" );
2160                 out.write( "GO term" );
2161                 out.write( "</b></td><td><b>" );
2162                 out.write( "GO namespace" );
2163                 out.write( "</b></td>" );
2164                 out.write( "</tr>" );
2165                 out.write( SurfacingConstants.NL );
2166                 out.write( "</tr>" );
2167                 out.write( SurfacingConstants.NL );
2168                 per_node_counter = 0;
2169                 if ( matrix.getNumberOfCharacters() > 0 ) {
2170                     per_node_go_mapped_domain_gain_loss_outfile = new File( per_node_go_mapped_domain_gain_loss_files_base_dir
2171                             + ForesterUtil.FILE_SEPARATOR + id + suffix_for_per_node_events_file );
2172                     SurfacingUtil.checkForOutputFileWriteability( per_node_go_mapped_domain_gain_loss_outfile );
2173                     per_node_go_mapped_domain_gain_loss_outfile_writer = ForesterUtil
2174                             .createBufferedWriter( per_node_go_mapped_domain_gain_loss_outfile );
2175                 }
2176                 else {
2177                     per_node_go_mapped_domain_gain_loss_outfile = null;
2178                     per_node_go_mapped_domain_gain_loss_outfile_writer = null;
2179                 }
2180                 for( int c = 0; c < matrix.getNumberOfCharacters(); ++c ) {
2181                     // Not nice:
2182                     // using null to indicate either UNCHANGED_PRESENT or GAIN.
2183                     if ( ( matrix.getState( id, c ) == state ) || ( ( state == null )
2184                             && ( ( matrix.getState( id, c ) == CharacterStateMatrix.GainLossStates.UNCHANGED_PRESENT )
2185                                     || ( matrix.getState( id, c ) == CharacterStateMatrix.GainLossStates.GAIN ) ) ) ) {
2186                         final String character = matrix.getCharacter( c );
2187                         String domain_0 = "";
2188                         String domain_1 = "";
2189                         if ( character.indexOf( BinaryDomainCombination.SEPARATOR ) > 0 ) {
2190                             final String[] s = character.split( BinaryDomainCombination.SEPARATOR );
2191                             if ( s.length != 2 ) {
2192                                 throw new AssertionError( "this should not have happened: unexpected format for domain combination: ["
2193                                         + character + "]" );
2194                             }
2195                             domain_0 = s[ 0 ];
2196                             domain_1 = s[ 1 ];
2197                         }
2198                         else {
2199                             domain_0 = character;
2200                         }
2201                         writeDomainData( domain_id_to_go_ids_map,
2202                                          go_id_to_term_map,
2203                                          go_namespace_limit,
2204                                          out,
2205                                          domain_0,
2206                                          domain_1,
2207                                          prefix_for_html,
2208                                          character_separator,
2209                                          domain_id_to_secondary_features_maps,
2210                                          null );
2211                         all_pfams_encountered.add( domain_0 );
2212                         if ( pfams_gained_or_lost != null ) {
2213                             pfams_gained_or_lost.add( domain_0 );
2214                         }
2215                         if ( !ForesterUtil.isEmpty( domain_1 ) ) {
2216                             all_pfams_encountered.add( domain_1 );
2217                             if ( pfams_gained_or_lost != null ) {
2218                                 pfams_gained_or_lost.add( domain_1 );
2219                             }
2220                         }
2221                         if ( per_node_go_mapped_domain_gain_loss_outfile_writer != null ) {
2222                             writeDomainsToIndividualFilePerTreeNode( per_node_go_mapped_domain_gain_loss_outfile_writer,
2223                                                                      domain_0,
2224                                                                      domain_1 );
2225                             per_node_counter++;
2226                         }
2227                     }
2228                 }
2229                 if ( per_node_go_mapped_domain_gain_loss_outfile_writer != null ) {
2230                     per_node_go_mapped_domain_gain_loss_outfile_writer.close();
2231                     if ( per_node_counter < 1 ) {
2232                         per_node_go_mapped_domain_gain_loss_outfile.delete();
2233                     }
2234                     per_node_counter = 0;
2235                 }
2236                 out.write( "</table>" );
2237                 out.write( SurfacingConstants.NL );
2238                 out.write( "<hr>" );
2239                 out.write( SurfacingConstants.NL );
2240             } // for( final String id : sorted_ids ) {
2241             out.write( "</body>" );
2242             out.write( SurfacingConstants.NL );
2243             out.write( "</html>" );
2244             out.write( SurfacingConstants.NL );
2245             out.flush();
2246             out.close();
2247         }
2248         catch ( final IOException e ) {
2249             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getMessage() );
2250         }
2251         ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote characters detailed HTML list: \"" + filename + "\"" );
2252     }
2253
2254     public static void writeDomainCombinationsCountsFile( final String[][] input_file_properties,
2255                                                           final File output_dir,
2256                                                           final Writer per_genome_domain_promiscuity_statistics_writer,
2257                                                           final GenomeWideCombinableDomains gwcd,
2258                                                           final int i,
2259                                                           final GenomeWideCombinableDomains.GenomeWideCombinableDomainsSortOrder dc_sort_order ) {
2260         File dc_outfile = new File( input_file_properties[ i ][ 1 ]
2261                 + surfacing.DOMAIN_COMBINITON_COUNTS_OUTPUTFILE_SUFFIX );
2262         if ( output_dir != null ) {
2263             dc_outfile = new File( output_dir + ForesterUtil.FILE_SEPARATOR + dc_outfile );
2264         }
2265         checkForOutputFileWriteability( dc_outfile );
2266         try {
2267             final BufferedWriter out = new BufferedWriter( new FileWriter( dc_outfile ) );
2268             out.write( gwcd.toStringBuilder( dc_sort_order ).toString() );
2269             out.close();
2270         }
2271         catch ( final IOException e ) {
2272             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getMessage() );
2273         }
2274         final DescriptiveStatistics stats = gwcd.getPerGenomeDomainPromiscuityStatistics();
2275         try {
2276             per_genome_domain_promiscuity_statistics_writer.write( input_file_properties[ i ][ 1 ] + "\t" );
2277             per_genome_domain_promiscuity_statistics_writer
2278                     .write( FORMATTER_3.format( stats.arithmeticMean() ) + "\t" );
2279             if ( stats.getN() < 2 ) {
2280                 per_genome_domain_promiscuity_statistics_writer.write( "n/a" + "\t" );
2281             }
2282             else {
2283                 per_genome_domain_promiscuity_statistics_writer
2284                         .write( FORMATTER_3.format( stats.sampleStandardDeviation() ) + "\t" );
2285             }
2286             per_genome_domain_promiscuity_statistics_writer.write( FORMATTER_3.format( stats.median() ) + "\t" );
2287             per_genome_domain_promiscuity_statistics_writer.write( ( int ) stats.getMin() + "\t" );
2288             per_genome_domain_promiscuity_statistics_writer.write( ( int ) stats.getMax() + "\t" );
2289             per_genome_domain_promiscuity_statistics_writer.write( stats.getN() + "\t" );
2290             final SortedSet<String> mpds = gwcd.getMostPromiscuosDomain();
2291             for( final String mpd : mpds ) {
2292                 per_genome_domain_promiscuity_statistics_writer.write( mpd + " " );
2293             }
2294             per_genome_domain_promiscuity_statistics_writer.write( ForesterUtil.LINE_SEPARATOR );
2295         }
2296         catch ( final IOException e ) {
2297             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getMessage() );
2298         }
2299         if ( input_file_properties[ i ].length == 3 ) {
2300             ForesterUtil.programMessage( surfacing.PRG_NAME,
2301                                          "Wrote domain combination counts for \"" + input_file_properties[ i ][ 0 ]
2302                                                  + "\" (" + input_file_properties[ i ][ 1 ] + ", "
2303                                                  + input_file_properties[ i ][ 2 ] + ") to: \"" + dc_outfile + "\"" );
2304         }
2305         else {
2306             ForesterUtil.programMessage( surfacing.PRG_NAME,
2307                                          "Wrote domain combination counts for \"" + input_file_properties[ i ][ 0 ]
2308                                                  + "\" (" + input_file_properties[ i ][ 1 ] + ") to: \"" + dc_outfile
2309                                                  + "\"" );
2310         }
2311     }
2312
2313     public static void writeDomainSimilaritiesToFile( final StringBuilder html_desc,
2314                                                       final StringBuilder html_title,
2315                                                       final Writer simple_tab_writer,
2316                                                       final Writer single_writer,
2317                                                       Map<Character, Writer> split_writers,
2318                                                       final SortedSet<DomainSimilarity> similarities,
2319                                                       final boolean treat_as_binary,
2320                                                       final List<Species> species_order,
2321                                                       final DomainSimilarity.PRINT_OPTION print_option,
2322                                                       final DomainSimilarity.DomainSimilarityScoring scoring,
2323                                                       final boolean verbose,
2324                                                       final Map<String, Integer> tax_code_to_id_map,
2325                                                       final Phylogeny phy,
2326                                                       final Set<String> pos_filter_doms )
2327             throws IOException {
2328         if ( ( single_writer != null ) && ( ( split_writers == null ) || split_writers.isEmpty() ) ) {
2329             split_writers = new HashMap<Character, Writer>();
2330             split_writers.put( '_', single_writer );
2331         }
2332         switch ( print_option ) {
2333             case SIMPLE_TAB_DELIMITED:
2334                 break;
2335             case HTML:
2336                 for( final Character key : split_writers.keySet() ) {
2337                     final Writer w = split_writers.get( key );
2338                     w.write( "<html>" );
2339                     w.write( SurfacingConstants.NL );
2340                     if ( key != '_' ) {
2341                         writeHtmlHead( w, "DC analysis (" + html_title + ") " + key.toString().toUpperCase() );
2342                     }
2343                     else {
2344                         writeHtmlHead( w, "DC analysis (" + html_title + ")" );
2345                     }
2346                     w.write( SurfacingConstants.NL );
2347                     w.write( "<body>" );
2348                     w.write( SurfacingConstants.NL );
2349                     w.write( html_desc.toString() );
2350                     w.write( SurfacingConstants.NL );
2351                     w.write( "<hr>" );
2352                     w.write( SurfacingConstants.NL );
2353                     w.write( "<br>" );
2354                     w.write( SurfacingConstants.NL );
2355                     w.write( "<table>" );
2356                     w.write( SurfacingConstants.NL );
2357                     w.write( "<tr><td><b>Domains:</b></td></tr>" );
2358                     w.write( SurfacingConstants.NL );
2359                 }
2360                 break;
2361         }
2362         //
2363         for( final DomainSimilarity similarity : similarities ) {
2364             if ( ( species_order != null ) && !species_order.isEmpty() ) {
2365                 ( similarity ).setSpeciesOrder( species_order );
2366             }
2367             if ( single_writer != null ) {
2368                 if ( !ForesterUtil.isEmpty( pos_filter_doms )
2369                         && pos_filter_doms.contains( similarity.getDomainId() ) ) {
2370                     single_writer.write( "<tr><td><b><a href=\"#" + similarity.getDomainId()
2371                             + "\"><span style=\"color:#00ff00\">" + similarity.getDomainId()
2372                             + "</span></a></b></td></tr>" );
2373                 }
2374                 else {
2375                     single_writer.write( "<tr><td><b><a href=\"#" + similarity.getDomainId() + "\">"
2376                             + similarity.getDomainId() + "</a></b></td></tr>" );
2377                 }
2378                 single_writer.write( SurfacingConstants.NL );
2379             }
2380             else {
2381                 Writer local_writer = split_writers
2382                         .get( ( similarity.getDomainId().charAt( 0 ) + "" ).toLowerCase().charAt( 0 ) );
2383                 if ( local_writer == null ) {
2384                     local_writer = split_writers.get( '0' );
2385                 }
2386                 if ( !ForesterUtil.isEmpty( pos_filter_doms )
2387                         && pos_filter_doms.contains( similarity.getDomainId() ) ) {
2388                     local_writer.write( "<tr><td><b><a href=\"#" + similarity.getDomainId()
2389                             + "\"><span style=\"color:#00ff00\">" + similarity.getDomainId()
2390                             + "</span></a></b></td></tr>" );
2391                 }
2392                 else {
2393                     local_writer.write( "<tr><td><b><a href=\"#" + similarity.getDomainId() + "\">"
2394                             + similarity.getDomainId() + "</a></b></td></tr>" );
2395                 }
2396                 local_writer.write( SurfacingConstants.NL );
2397             }
2398         }
2399         for( final Writer w : split_writers.values() ) {
2400             w.write( "</table>" );
2401             w.write( SurfacingConstants.NL );
2402             w.write( "<hr>" );
2403             w.write( SurfacingConstants.NL );
2404             //
2405             w.write( "<table>" );
2406             w.write( SurfacingConstants.NL );
2407             w.write( "<tr><td><b>" );
2408             w.write( "Species group colors:" );
2409             w.write( "</b></td></tr>" );
2410             w.write( SurfacingConstants.NL );
2411             writeColorLabels( "Deuterostomia", TaxonomyColors.DEUTEROSTOMIA_COLOR, w );
2412             writeColorLabels( "Protostomia", TaxonomyColors.PROTOSTOMIA_COLOR, w );
2413             writeColorLabels( "Cnidaria", TaxonomyColors.CNIDARIA_COLOR, w );
2414             writeColorLabels( "Placozoa", TaxonomyColors.PLACOZOA_COLOR, w );
2415             writeColorLabels( "Ctenophora (comb jellies)", TaxonomyColors.CTENOPHORA_COLOR, w );
2416             writeColorLabels( "Porifera (sponges)", TaxonomyColors.PORIFERA_COLOR, w );
2417             writeColorLabels( "Choanoflagellida", TaxonomyColors.CHOANOFLAGELLIDA, w );
2418             writeColorLabels( "Ichthyosporea & Filasterea", TaxonomyColors.ICHTHYOSPOREA_AND_FILASTEREA, w );
2419             writeColorLabels( "Dikarya (Ascomycota & Basidiomycota, so-called \"higher fungi\")",
2420                               TaxonomyColors.DIKARYA_COLOR,
2421                               w );
2422             writeColorLabels( "other Fungi", TaxonomyColors.OTHER_FUNGI_COLOR, w );
2423             writeColorLabels( "Nucleariidae and Fonticula group",
2424                               TaxonomyColors.NUCLEARIIDAE_AND_FONTICULA_GROUP_COLOR,
2425                               w );
2426             writeColorLabels( "Amoebozoa", TaxonomyColors.AMOEBOZOA_COLOR, w );
2427             writeColorLabels( "Embryophyta (plants)", TaxonomyColors.EMBRYOPHYTA_COLOR, w );
2428             writeColorLabels( "Chlorophyta (green algae)", TaxonomyColors.CHLOROPHYTA_COLOR, w );
2429             writeColorLabels( "Rhodophyta (red algae)", TaxonomyColors.RHODOPHYTA_COLOR, w );
2430             writeColorLabels( "Glaucocystophyce (Glaucophyta)", TaxonomyColors.GLAUCOPHYTA_COLOR, w );
2431             writeColorLabels( "Hacrobia (Cryptophyta & Haptophyceae & Centroheliozoa)",
2432                               TaxonomyColors.HACROBIA_COLOR,
2433                               w );
2434             writeColorLabels( "Stramenopiles (Chromophyta, heterokonts)", TaxonomyColors.STRAMENOPILES_COLOR, w );
2435             writeColorLabels( "Alveolata", TaxonomyColors.ALVEOLATA_COLOR, w );
2436             writeColorLabels( "Rhizaria", TaxonomyColors.RHIZARIA_COLOR, w );
2437             writeColorLabels( "Excavata", TaxonomyColors.EXCAVATA_COLOR, w );
2438             writeColorLabels( "Apusozoa", TaxonomyColors.APUSOZOA_COLOR, w );
2439             writeColorLabels( "Archaea", TaxonomyColors.ARCHAEA_COLOR, w );
2440             writeColorLabels( "Bacteria", TaxonomyColors.BACTERIA_COLOR, w );
2441             w.write( "</table>" );
2442             w.write( SurfacingConstants.NL );
2443             //
2444             w.write( "<hr>" );
2445             w.write( SurfacingConstants.NL );
2446             w.write( "<table>" );
2447             w.write( SurfacingConstants.NL );
2448         }
2449         //
2450         for( final DomainSimilarity similarity : similarities ) {
2451             if ( ( species_order != null ) && !species_order.isEmpty() ) {
2452                 ( similarity ).setSpeciesOrder( species_order );
2453             }
2454             if ( simple_tab_writer != null ) {
2455                 simple_tab_writer.write( similarity
2456                         .toStringBuffer( PRINT_OPTION.SIMPLE_TAB_DELIMITED, tax_code_to_id_map, null ).toString() );
2457             }
2458             if ( single_writer != null ) {
2459                 single_writer.write( similarity.toStringBuffer( print_option, tax_code_to_id_map, phy ).toString() );
2460                 single_writer.write( SurfacingConstants.NL );
2461             }
2462             else {
2463                 Writer local_writer = split_writers
2464                         .get( ( similarity.getDomainId().charAt( 0 ) + "" ).toLowerCase().charAt( 0 ) );
2465                 if ( local_writer == null ) {
2466                     local_writer = split_writers.get( '0' );
2467                 }
2468                 local_writer.write( similarity.toStringBuffer( print_option, tax_code_to_id_map, phy ).toString() );
2469                 local_writer.write( SurfacingConstants.NL );
2470             }
2471         }
2472         switch ( print_option ) {
2473             case HTML:
2474                 for( final Writer w : split_writers.values() ) {
2475                     w.write( SurfacingConstants.NL );
2476                     w.write( "</table>" );
2477                     w.write( SurfacingConstants.NL );
2478                     w.write( "</font>" );
2479                     w.write( SurfacingConstants.NL );
2480                     w.write( "</body>" );
2481                     w.write( SurfacingConstants.NL );
2482                     w.write( "</html>" );
2483                     w.write( SurfacingConstants.NL );
2484                 }
2485                 break;
2486             default:
2487                 break;
2488         }
2489         for( final Writer w : split_writers.values() ) {
2490             w.close();
2491         }
2492     }
2493
2494     public static void writeHtmlHead( final Writer w, final String title ) throws IOException {
2495         w.write( SurfacingConstants.NL );
2496         w.write( "<head>" );
2497         w.write( "<title>" );
2498         w.write( title );
2499         w.write( "</title>" );
2500         w.write( SurfacingConstants.NL );
2501         w.write( "<style>" );
2502         w.write( SurfacingConstants.NL );
2503         w.write( "a:visited { color : #000066; text-decoration : none; }" );
2504         w.write( SurfacingConstants.NL );
2505         w.write( "a:link { color : #000066; text-decoration : none; }" );
2506         w.write( SurfacingConstants.NL );
2507         w.write( "a:active { color : ##000066; text-decoration : none; }" );
2508         w.write( SurfacingConstants.NL );
2509         w.write( "a:hover { color : #FFFFFF; background-color : #000000; text-decoration : none; }" );
2510         w.write( SurfacingConstants.NL );
2511         //
2512         w.write( "a.pl:visited { color : #505050; text-decoration : none; font-size: 7px;}" );
2513         w.write( SurfacingConstants.NL );
2514         w.write( "a.pl:link { color : #505050; text-decoration : none; font-size: 7px;}" );
2515         w.write( SurfacingConstants.NL );
2516         w.write( "a.pl:active { color : #505050; text-decoration : none; font-size: 7px;}" );
2517         w.write( SurfacingConstants.NL );
2518         w.write( "a.pl:hover { color : #FFFFFF; background-color : #000000; text-decoration : none; font-size: 7px;}" );
2519         w.write( SurfacingConstants.NL );
2520         //
2521         w.write( "a.ps:visited { color : #707070; text-decoration : none; font-size: 7px;}" );
2522         w.write( SurfacingConstants.NL );
2523         w.write( "a.ps:link { color : #707070; text-decoration : none; font-size: 7px;}" );
2524         w.write( SurfacingConstants.NL );
2525         w.write( "a.ps:active { color : #707070; text-decoration : none; font-size: 7px;}" );
2526         w.write( SurfacingConstants.NL );
2527         w.write( "a.ps:hover { color : #FFFFFF; background-color : #000000; text-decoration : none; font-size: 7px;}" );
2528         w.write( SurfacingConstants.NL );
2529         //
2530         w.write( "td { text-align: left; vertical-align: top; font-family: Verdana, Arial, Helvetica; font-size: 8pt}" );
2531         w.write( SurfacingConstants.NL );
2532         w.write( "h1 { color : #0000FF; font-family: Verdana, Arial, Helvetica; font-size: 18pt; font-weight: bold }" );
2533         w.write( SurfacingConstants.NL );
2534         w.write( "h2 { color : #0000FF; font-family: Verdana, Arial, Helvetica; font-size: 16pt; font-weight: bold }" );
2535         w.write( SurfacingConstants.NL );
2536         w.write( "</style>" );
2537         w.write( SurfacingConstants.NL );
2538         w.write( "</head>" );
2539         w.write( SurfacingConstants.NL );
2540     }
2541
2542     public static void writeMatrixToFile( final CharacterStateMatrix<?> matrix,
2543                                           final String filename,
2544                                           final Format format ) {
2545         final File outfile = new File( filename );
2546         checkForOutputFileWriteability( outfile );
2547         try {
2548             final BufferedWriter out = new BufferedWriter( new FileWriter( outfile ) );
2549             matrix.toWriter( out, format );
2550             out.flush();
2551             out.close();
2552         }
2553         catch ( final IOException e ) {
2554             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getMessage() );
2555         }
2556         ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote matrix: \"" + filename + "\"" );
2557     }
2558
2559     public static void writeMatrixToFile( final File matrix_outfile, final List<DistanceMatrix> matrices ) {
2560         checkForOutputFileWriteability( matrix_outfile );
2561         try {
2562             final BufferedWriter out = new BufferedWriter( new FileWriter( matrix_outfile ) );
2563             for( final DistanceMatrix distance_matrix : matrices ) {
2564                 out.write( distance_matrix.toStringBuffer( DistanceMatrix.Format.PHYLIP ).toString() );
2565                 out.write( ForesterUtil.LINE_SEPARATOR );
2566                 out.flush();
2567             }
2568             out.close();
2569         }
2570         catch ( final IOException e ) {
2571             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getMessage() );
2572         }
2573         ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote distance matrices to \"" + matrix_outfile + "\"" );
2574     }
2575
2576     public static void writePhylogenyToFile( final Phylogeny phylogeny, final String filename ) {
2577         final PhylogenyWriter writer = new PhylogenyWriter();
2578         try {
2579             writer.toPhyloXML( new File( filename ), phylogeny, 1 );
2580         }
2581         catch ( final IOException e ) {
2582             ForesterUtil.printWarningMessage( surfacing.PRG_NAME,
2583                                               "failed to write phylogeny to \"" + filename + "\": " + e );
2584         }
2585         ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote phylogeny to \"" + filename + "\"" );
2586     }
2587
2588     public static void writePresentToNexus( final File output_file,
2589                                             final File positive_filter_file,
2590                                             final SortedSet<String> filter,
2591                                             final List<GenomeWideCombinableDomains> gwcd_list ) {
2592         try {
2593             writeMatrixToFile( DomainParsimonyCalculator
2594                     .createMatrixOfDomainPresenceOrAbsence( gwcd_list, positive_filter_file == null ? null : filter ),
2595                                output_file + surfacing.DOMAINS_PRESENT_NEXUS,
2596                                Format.NEXUS_BINARY );
2597             writeMatrixToFile( DomainParsimonyCalculator
2598                     .createMatrixOfBinaryDomainCombinationPresenceOrAbsence( gwcd_list ),
2599                                output_file + surfacing.BDC_PRESENT_NEXUS,
2600                                Format.NEXUS_BINARY );
2601         }
2602         catch ( final Exception e ) {
2603             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getLocalizedMessage() );
2604         }
2605     }
2606
2607     public static void writeProteinListsForAllSpecies( final File output_dir,
2608                                                        final SortedMap<Species, List<Protein>> protein_lists_per_species,
2609                                                        final List<GenomeWideCombinableDomains> gwcd_list,
2610                                                        final double domain_e_cutoff,
2611                                                        final Set<String> pos_filter_doms ) {
2612         final SortedSet<String> all_domains = new TreeSet<String>();
2613         for( final GenomeWideCombinableDomains gwcd : gwcd_list ) {
2614             all_domains.addAll( gwcd.getAllDomainIds() );
2615         }
2616         for( final String domain : all_domains ) {
2617             if ( !ForesterUtil.isEmpty( pos_filter_doms ) && !pos_filter_doms.contains( domain ) ) {
2618                 continue;
2619             }
2620             final File out = new File( output_dir + ForesterUtil.FILE_SEPARATOR + domain
2621                     + surfacing.SEQ_EXTRACT_SUFFIX );
2622             checkForOutputFileWriteability( out );
2623             try {
2624                 final Writer proteins_file_writer = new BufferedWriter( new FileWriter( out ) );
2625                 extractProteinNames( protein_lists_per_species,
2626                                      domain,
2627                                      proteins_file_writer,
2628                                      "\t",
2629                                      surfacing.LIMIT_SPEC_FOR_PROT_EX,
2630                                      domain_e_cutoff );
2631                 proteins_file_writer.close();
2632             }
2633             catch ( final IOException e ) {
2634                 ForesterUtil.fatalError( surfacing.PRG_NAME, e.getLocalizedMessage() );
2635             }
2636             ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote proteins list to \"" + out + "\"" );
2637         }
2638     }
2639
2640     public static void writeTaxonomyLinks( final Writer writer,
2641                                            final String species,
2642                                            final Map<String, Integer> tax_code_to_id_map )
2643             throws IOException {
2644         if ( ( species.length() > 1 ) && ( species.indexOf( '_' ) < 1 ) ) {
2645             writer.write( " [" );
2646             if ( ( tax_code_to_id_map != null ) && tax_code_to_id_map.containsKey( species ) ) {
2647                 writer.write( "<a href=\"" + SurfacingConstants.UNIPROT_TAXONOMY_ID_LINK
2648                         + tax_code_to_id_map.get( species ) + "\" target=\"taxonomy_window\">uniprot</a>" );
2649             }
2650             else {
2651                 writer.write( "<a href=\"" + SurfacingConstants.EOL_LINK + species
2652                         + "\" target=\"taxonomy_window\">eol</a>" );
2653                 writer.write( "|" );
2654                 writer.write( "<a href=\"" + SurfacingConstants.GOOGLE_SCHOLAR_SEARCH + species
2655                         + "\" target=\"taxonomy_window\">scholar</a>" );
2656                 writer.write( "|" );
2657                 writer.write( "<a href=\"" + SurfacingConstants.GOOGLE_WEB_SEARCH_LINK + species
2658                         + "\" target=\"taxonomy_window\">google</a>" );
2659             }
2660             writer.write( "]" );
2661         }
2662     }
2663
2664     private final static void addToCountMap( final Map<String, Integer> map, final String s ) {
2665         if ( map.containsKey( s ) ) {
2666             map.put( s, map.get( s ) + 1 );
2667         }
2668         else {
2669             map.put( s, 1 );
2670         }
2671     }
2672
2673     private static void calculateIndependentDomainCombinationGains( final Phylogeny local_phylogeny_l,
2674                                                                     final String outfilename_for_counts,
2675                                                                     final String outfilename_for_dc,
2676                                                                     final String outfilename_for_dc_for_go_mapping,
2677                                                                     final String outfilename_for_dc_for_go_mapping_unique,
2678                                                                     final String outfilename_for_rank_counts,
2679                                                                     final String outfilename_for_ancestor_species_counts,
2680                                                                     final String outfilename_for_protein_stats,
2681                                                                     final Map<String, DescriptiveStatistics> protein_length_stats_by_dc,
2682                                                                     final Map<String, DescriptiveStatistics> domain_number_stats_by_dc,
2683                                                                     final Map<String, DescriptiveStatistics> domain_length_stats_by_domain ) {
2684         try {
2685             //
2686             //            if ( protein_length_stats_by_dc != null ) {
2687             //                for( final Entry<?, DescriptiveStatistics> entry : protein_length_stats_by_dc.entrySet() ) {
2688             //                    System.out.print( entry.getKey().toString() );
2689             //                    System.out.print( ": " );
2690             //                    double[] a = entry.getValue().getDataAsDoubleArray();
2691             //                    for( int i = 0; i < a.length; i++ ) {
2692             //                        System.out.print( a[ i ] + " " );
2693             //                    }
2694             //                    System.out.println();
2695             //                }
2696             //            }
2697             //            if ( domain_number_stats_by_dc != null ) {
2698             //                for( final Entry<?, DescriptiveStatistics> entry : domain_number_stats_by_dc.entrySet() ) {
2699             //                    System.out.print( entry.getKey().toString() );
2700             //                    System.out.print( ": " );
2701             //                    double[] a = entry.getValue().getDataAsDoubleArray();
2702             //                    for( int i = 0; i < a.length; i++ ) {
2703             //                        System.out.print( a[ i ] + " " );
2704             //                    }
2705             //                    System.out.println();
2706             //                }
2707             //            }
2708             //
2709             final BufferedWriter out_counts = new BufferedWriter( new FileWriter( outfilename_for_counts ) );
2710             final BufferedWriter out_dc = new BufferedWriter( new FileWriter( outfilename_for_dc ) );
2711             final BufferedWriter out_dc_for_go_mapping = new BufferedWriter( new FileWriter( outfilename_for_dc_for_go_mapping ) );
2712             final BufferedWriter out_dc_for_go_mapping_unique = new BufferedWriter( new FileWriter( outfilename_for_dc_for_go_mapping_unique ) );
2713             final SortedMap<String, Integer> dc_gain_counts = new TreeMap<String, Integer>();
2714             for( final PhylogenyNodeIterator it = local_phylogeny_l.iteratorPostorder(); it.hasNext(); ) {
2715                 final PhylogenyNode n = it.next();
2716                 final Set<String> gained_dc = n.getNodeData().getBinaryCharacters().getGainedCharacters();
2717                 for( final String dc : gained_dc ) {
2718                     if ( dc_gain_counts.containsKey( dc ) ) {
2719                         dc_gain_counts.put( dc, dc_gain_counts.get( dc ) + 1 );
2720                     }
2721                     else {
2722                         dc_gain_counts.put( dc, 1 );
2723                     }
2724                 }
2725             }
2726             final SortedMap<Integer, Integer> histogram = new TreeMap<Integer, Integer>();
2727             final SortedMap<Integer, StringBuilder> domain_lists = new TreeMap<Integer, StringBuilder>();
2728             final SortedMap<Integer, DescriptiveStatistics> dc_reapp_counts_to_protein_length_stats = new TreeMap<Integer, DescriptiveStatistics>();
2729             final SortedMap<Integer, DescriptiveStatistics> dc_reapp_counts_to_domain_number_stats = new TreeMap<Integer, DescriptiveStatistics>();
2730             final SortedMap<Integer, DescriptiveStatistics> dc_reapp_counts_to_domain_lengths_stats = new TreeMap<Integer, DescriptiveStatistics>();
2731             final SortedMap<Integer, PriorityQueue<String>> domain_lists_go = new TreeMap<Integer, PriorityQueue<String>>();
2732             final SortedMap<Integer, SortedSet<String>> domain_lists_go_unique = new TreeMap<Integer, SortedSet<String>>();
2733             final Set<String> dcs = dc_gain_counts.keySet();
2734             final SortedSet<String> more_than_once = new TreeSet<String>();
2735             DescriptiveStatistics gained_once_lengths_stats = new BasicDescriptiveStatistics();
2736             DescriptiveStatistics gained_once_domain_count_stats = new BasicDescriptiveStatistics();
2737             DescriptiveStatistics gained_multiple_times_lengths_stats = new BasicDescriptiveStatistics();
2738             final DescriptiveStatistics gained_multiple_times_domain_count_stats = new BasicDescriptiveStatistics();
2739             long gained_multiple_times_domain_length_sum = 0;
2740             long gained_once_domain_length_sum = 0;
2741             long gained_multiple_times_domain_length_count = 0;
2742             long gained_once_domain_length_count = 0;
2743             for( final String dc : dcs ) {
2744                 final int count = dc_gain_counts.get( dc );
2745                 if ( histogram.containsKey( count ) ) {
2746                     histogram.put( count, histogram.get( count ) + 1 );
2747                     domain_lists.get( count ).append( ", " + dc );
2748                     domain_lists_go.get( count ).addAll( splitDomainCombination( dc ) );
2749                     domain_lists_go_unique.get( count ).addAll( splitDomainCombination( dc ) );
2750                 }
2751                 else {
2752                     histogram.put( count, 1 );
2753                     domain_lists.put( count, new StringBuilder( dc ) );
2754                     final PriorityQueue<String> q = new PriorityQueue<String>();
2755                     q.addAll( splitDomainCombination( dc ) );
2756                     domain_lists_go.put( count, q );
2757                     final SortedSet<String> set = new TreeSet<String>();
2758                     set.addAll( splitDomainCombination( dc ) );
2759                     domain_lists_go_unique.put( count, set );
2760                 }
2761                 if ( protein_length_stats_by_dc != null ) {
2762                     if ( !dc_reapp_counts_to_protein_length_stats.containsKey( count ) ) {
2763                         dc_reapp_counts_to_protein_length_stats.put( count, new BasicDescriptiveStatistics() );
2764                     }
2765                     dc_reapp_counts_to_protein_length_stats.get( count )
2766                             .addValue( protein_length_stats_by_dc.get( dc ).arithmeticMean() );
2767                 }
2768                 if ( domain_number_stats_by_dc != null ) {
2769                     if ( !dc_reapp_counts_to_domain_number_stats.containsKey( count ) ) {
2770                         dc_reapp_counts_to_domain_number_stats.put( count, new BasicDescriptiveStatistics() );
2771                     }
2772                     dc_reapp_counts_to_domain_number_stats.get( count )
2773                             .addValue( domain_number_stats_by_dc.get( dc ).arithmeticMean() );
2774                 }
2775                 if ( domain_length_stats_by_domain != null ) {
2776                     if ( !dc_reapp_counts_to_domain_lengths_stats.containsKey( count ) ) {
2777                         dc_reapp_counts_to_domain_lengths_stats.put( count, new BasicDescriptiveStatistics() );
2778                     }
2779                     final String[] ds = dc.split( "=" );
2780                     dc_reapp_counts_to_domain_lengths_stats.get( count )
2781                             .addValue( domain_length_stats_by_domain.get( ds[ 0 ] ).arithmeticMean() );
2782                     dc_reapp_counts_to_domain_lengths_stats.get( count )
2783                             .addValue( domain_length_stats_by_domain.get( ds[ 1 ] ).arithmeticMean() );
2784                 }
2785                 if ( count > 1 ) {
2786                     more_than_once.add( dc );
2787                     if ( protein_length_stats_by_dc != null ) {
2788                         final DescriptiveStatistics s = protein_length_stats_by_dc.get( dc );
2789                         for( final double element : s.getData() ) {
2790                             gained_multiple_times_lengths_stats.addValue( element );
2791                         }
2792                     }
2793                     if ( domain_number_stats_by_dc != null ) {
2794                         final DescriptiveStatistics s = domain_number_stats_by_dc.get( dc );
2795                         for( final double element : s.getData() ) {
2796                             gained_multiple_times_domain_count_stats.addValue( element );
2797                         }
2798                     }
2799                     if ( domain_length_stats_by_domain != null ) {
2800                         final String[] ds = dc.split( "=" );
2801                         final DescriptiveStatistics s0 = domain_length_stats_by_domain.get( ds[ 0 ] );
2802                         final DescriptiveStatistics s1 = domain_length_stats_by_domain.get( ds[ 1 ] );
2803                         for( final double element : s0.getData() ) {
2804                             gained_multiple_times_domain_length_sum += element;
2805                             ++gained_multiple_times_domain_length_count;
2806                         }
2807                         for( final double element : s1.getData() ) {
2808                             gained_multiple_times_domain_length_sum += element;
2809                             ++gained_multiple_times_domain_length_count;
2810                         }
2811                     }
2812                 }
2813                 else {
2814                     if ( protein_length_stats_by_dc != null ) {
2815                         final DescriptiveStatistics s = protein_length_stats_by_dc.get( dc );
2816                         for( final double element : s.getData() ) {
2817                             gained_once_lengths_stats.addValue( element );
2818                         }
2819                     }
2820                     if ( domain_number_stats_by_dc != null ) {
2821                         final DescriptiveStatistics s = domain_number_stats_by_dc.get( dc );
2822                         for( final double element : s.getData() ) {
2823                             gained_once_domain_count_stats.addValue( element );
2824                         }
2825                     }
2826                     if ( domain_length_stats_by_domain != null ) {
2827                         final String[] ds = dc.split( "=" );
2828                         final DescriptiveStatistics s0 = domain_length_stats_by_domain.get( ds[ 0 ] );
2829                         final DescriptiveStatistics s1 = domain_length_stats_by_domain.get( ds[ 1 ] );
2830                         for( final double element : s0.getData() ) {
2831                             gained_once_domain_length_sum += element;
2832                             ++gained_once_domain_length_count;
2833                         }
2834                         for( final double element : s1.getData() ) {
2835                             gained_once_domain_length_sum += element;
2836                             ++gained_once_domain_length_count;
2837                         }
2838                     }
2839                 }
2840             }
2841             final Set<Integer> histogram_keys = histogram.keySet();
2842             for( final Integer histogram_key : histogram_keys ) {
2843                 final int count = histogram.get( histogram_key );
2844                 final StringBuilder dc = domain_lists.get( histogram_key );
2845                 out_counts.write( histogram_key + "\t" + count + ForesterUtil.LINE_SEPARATOR );
2846                 out_dc.write( histogram_key + "\t" + dc + ForesterUtil.LINE_SEPARATOR );
2847                 out_dc_for_go_mapping.write( "#" + histogram_key + ForesterUtil.LINE_SEPARATOR );
2848                 final Object[] sorted = domain_lists_go.get( histogram_key ).toArray();
2849                 Arrays.sort( sorted );
2850                 for( final Object domain : sorted ) {
2851                     out_dc_for_go_mapping.write( domain + ForesterUtil.LINE_SEPARATOR );
2852                 }
2853                 out_dc_for_go_mapping_unique.write( "#" + histogram_key + ForesterUtil.LINE_SEPARATOR );
2854                 for( final String domain : domain_lists_go_unique.get( histogram_key ) ) {
2855                     out_dc_for_go_mapping_unique.write( domain + ForesterUtil.LINE_SEPARATOR );
2856                 }
2857             }
2858             out_counts.close();
2859             out_dc.close();
2860             out_dc_for_go_mapping.close();
2861             out_dc_for_go_mapping_unique.close();
2862             final SortedMap<String, Integer> lca_rank_counts = new TreeMap<String, Integer>();
2863             final SortedMap<String, Integer> lca_ancestor_species_counts = new TreeMap<String, Integer>();
2864             for( final String dc : more_than_once ) {
2865                 final List<PhylogenyNode> nodes = new ArrayList<PhylogenyNode>();
2866                 for( final PhylogenyNodeIterator it = local_phylogeny_l.iteratorExternalForward(); it.hasNext(); ) {
2867                     final PhylogenyNode n = it.next();
2868                     if ( n.getNodeData().getBinaryCharacters().getGainedCharacters().contains( dc ) ) {
2869                         nodes.add( n );
2870                     }
2871                 }
2872                 for( int i = 0; i < ( nodes.size() - 1 ); ++i ) {
2873                     for( int j = i + 1; j < nodes.size(); ++j ) {
2874                         final PhylogenyNode lca = PhylogenyMethods.calculateLCA( nodes.get( i ), nodes.get( j ) );
2875                         String rank = "unknown";
2876                         if ( lca.getNodeData().isHasTaxonomy()
2877                                 && !ForesterUtil.isEmpty( lca.getNodeData().getTaxonomy().getRank() ) ) {
2878                             rank = lca.getNodeData().getTaxonomy().getRank();
2879                         }
2880                         addToCountMap( lca_rank_counts, rank );
2881                         String lca_species;
2882                         if ( lca.getNodeData().isHasTaxonomy()
2883                                 && !ForesterUtil.isEmpty( lca.getNodeData().getTaxonomy().getScientificName() ) ) {
2884                             lca_species = lca.getNodeData().getTaxonomy().getScientificName();
2885                         }
2886                         else if ( lca.getNodeData().isHasTaxonomy()
2887                                 && !ForesterUtil.isEmpty( lca.getNodeData().getTaxonomy().getCommonName() ) ) {
2888                             lca_species = lca.getNodeData().getTaxonomy().getCommonName();
2889                         }
2890                         else {
2891                             lca_species = lca.getName();
2892                         }
2893                         addToCountMap( lca_ancestor_species_counts, lca_species );
2894                     }
2895                 }
2896             }
2897             final BufferedWriter out_for_rank_counts = new BufferedWriter( new FileWriter( outfilename_for_rank_counts ) );
2898             final BufferedWriter out_for_ancestor_species_counts = new BufferedWriter( new FileWriter( outfilename_for_ancestor_species_counts ) );
2899             ForesterUtil.map2writer( out_for_rank_counts, lca_rank_counts, "\t", ForesterUtil.LINE_SEPARATOR );
2900             ForesterUtil.map2writer( out_for_ancestor_species_counts,
2901                                      lca_ancestor_species_counts,
2902                                      "\t",
2903                                      ForesterUtil.LINE_SEPARATOR );
2904             out_for_rank_counts.close();
2905             out_for_ancestor_species_counts.close();
2906             if ( !ForesterUtil.isEmpty( outfilename_for_protein_stats ) && ( ( domain_length_stats_by_domain != null )
2907                     || ( protein_length_stats_by_dc != null ) || ( domain_number_stats_by_dc != null ) ) ) {
2908                 final BufferedWriter w = new BufferedWriter( new FileWriter( outfilename_for_protein_stats ) );
2909                 w.write( "Domain Lengths: " );
2910                 w.write( "\n" );
2911                 if ( domain_length_stats_by_domain != null ) {
2912                     for( final Entry<Integer, DescriptiveStatistics> entry : dc_reapp_counts_to_domain_lengths_stats
2913                             .entrySet() ) {
2914                         w.write( entry.getKey().toString() );
2915                         w.write( "\t" + entry.getValue().arithmeticMean() );
2916                         w.write( "\t" + entry.getValue().median() );
2917                         w.write( "\n" );
2918                     }
2919                 }
2920                 w.flush();
2921                 w.write( "\n" );
2922                 w.write( "\n" );
2923                 w.write( "Protein Lengths: " );
2924                 w.write( "\n" );
2925                 if ( protein_length_stats_by_dc != null ) {
2926                     for( final Entry<Integer, DescriptiveStatistics> entry : dc_reapp_counts_to_protein_length_stats
2927                             .entrySet() ) {
2928                         w.write( entry.getKey().toString() );
2929                         w.write( "\t" + entry.getValue().arithmeticMean() );
2930                         w.write( "\t" + entry.getValue().median() );
2931                         w.write( "\n" );
2932                     }
2933                 }
2934                 w.flush();
2935                 w.write( "\n" );
2936                 w.write( "\n" );
2937                 w.write( "Number of domains: " );
2938                 w.write( "\n" );
2939                 if ( domain_number_stats_by_dc != null ) {
2940                     for( final Entry<Integer, DescriptiveStatistics> entry : dc_reapp_counts_to_domain_number_stats
2941                             .entrySet() ) {
2942                         w.write( entry.getKey().toString() );
2943                         w.write( "\t" + entry.getValue().arithmeticMean() );
2944                         w.write( "\t" + entry.getValue().median() );
2945                         w.write( "\n" );
2946                     }
2947                 }
2948                 w.flush();
2949                 w.write( "\n" );
2950                 w.write( "\n" );
2951                 w.write( "Gained once, domain lengths:" );
2952                 w.write( "\n" );
2953                 w.write( "N: " + gained_once_domain_length_count );
2954                 w.write( "\n" );
2955                 w.write( "Avg: " + ( ( double ) gained_once_domain_length_sum / gained_once_domain_length_count ) );
2956                 w.write( "\n" );
2957                 w.write( "\n" );
2958                 w.write( "Gained multiple times, domain lengths:" );
2959                 w.write( "\n" );
2960                 w.write( "N: " + gained_multiple_times_domain_length_count );
2961                 w.write( "\n" );
2962                 w.write( "Avg: " + ( ( double ) gained_multiple_times_domain_length_sum
2963                         / gained_multiple_times_domain_length_count ) );
2964                 w.write( "\n" );
2965                 w.write( "\n" );
2966                 w.write( "\n" );
2967                 w.write( "\n" );
2968                 w.write( "Gained once, protein lengths:" );
2969                 w.write( "\n" );
2970                 w.write( gained_once_lengths_stats.toString() );
2971                 gained_once_lengths_stats = null;
2972                 w.write( "\n" );
2973                 w.write( "\n" );
2974                 w.write( "Gained once, domain counts:" );
2975                 w.write( "\n" );
2976                 w.write( gained_once_domain_count_stats.toString() );
2977                 gained_once_domain_count_stats = null;
2978                 w.write( "\n" );
2979                 w.write( "\n" );
2980                 w.write( "Gained multiple times, protein lengths:" );
2981                 w.write( "\n" );
2982                 w.write( gained_multiple_times_lengths_stats.toString() );
2983                 gained_multiple_times_lengths_stats = null;
2984                 w.write( "\n" );
2985                 w.write( "\n" );
2986                 w.write( "Gained multiple times, domain counts:" );
2987                 w.write( "\n" );
2988                 w.write( gained_multiple_times_domain_count_stats.toString() );
2989                 w.flush();
2990                 w.close();
2991             }
2992         }
2993         catch ( final IOException e ) {
2994             ForesterUtil.printWarningMessage( surfacing.PRG_NAME, "Failure to write: " + e );
2995         }
2996         ForesterUtil.programMessage( surfacing.PRG_NAME,
2997                                      "Wrote independent domain combination gains fitch counts to ["
2998                                              + outfilename_for_counts + "]" );
2999         ForesterUtil.programMessage( surfacing.PRG_NAME,
3000                                      "Wrote independent domain combination gains fitch lists to [" + outfilename_for_dc
3001                                              + "]" );
3002         ForesterUtil.programMessage( surfacing.PRG_NAME,
3003                                      "Wrote independent domain combination gains fitch lists to (for GO mapping) ["
3004                                              + outfilename_for_dc_for_go_mapping + "]" );
3005         ForesterUtil.programMessage( surfacing.PRG_NAME,
3006                                      "Wrote independent domain combination gains fitch lists to (for GO mapping, unique) ["
3007                                              + outfilename_for_dc_for_go_mapping_unique + "]" );
3008     }
3009
3010     private static SortedSet<String> collectAllDomainsChangedOnSubtree( final PhylogenyNode subtree_root,
3011                                                                         final boolean get_gains ) {
3012         final SortedSet<String> domains = new TreeSet<String>();
3013         for( final PhylogenyNode descendant : PhylogenyMethods.getAllDescendants( subtree_root ) ) {
3014             final BinaryCharacters chars = descendant.getNodeData().getBinaryCharacters();
3015             if ( get_gains ) {
3016                 domains.addAll( chars.getGainedCharacters() );
3017             }
3018             else {
3019                 domains.addAll( chars.getLostCharacters() );
3020             }
3021         }
3022         return domains;
3023     }
3024
3025     private static File createBaseDirForPerNodeDomainFiles( final String base_dir,
3026                                                             final boolean domain_combinations,
3027                                                             final CharacterStateMatrix.GainLossStates state,
3028                                                             final String outfile ) {
3029         File per_node_go_mapped_domain_gain_loss_files_base_dir = new File( new File( outfile ).getParent()
3030                 + ForesterUtil.FILE_SEPARATOR + base_dir );
3031         if ( !per_node_go_mapped_domain_gain_loss_files_base_dir.exists() ) {
3032             per_node_go_mapped_domain_gain_loss_files_base_dir.mkdir();
3033         }
3034         if ( domain_combinations ) {
3035             per_node_go_mapped_domain_gain_loss_files_base_dir = new File( per_node_go_mapped_domain_gain_loss_files_base_dir
3036                     + ForesterUtil.FILE_SEPARATOR + "DC" );
3037         }
3038         else {
3039             per_node_go_mapped_domain_gain_loss_files_base_dir = new File( per_node_go_mapped_domain_gain_loss_files_base_dir
3040                     + ForesterUtil.FILE_SEPARATOR + "DOMAINS" );
3041         }
3042         if ( !per_node_go_mapped_domain_gain_loss_files_base_dir.exists() ) {
3043             per_node_go_mapped_domain_gain_loss_files_base_dir.mkdir();
3044         }
3045         if ( state == GainLossStates.GAIN ) {
3046             per_node_go_mapped_domain_gain_loss_files_base_dir = new File( per_node_go_mapped_domain_gain_loss_files_base_dir
3047                     + ForesterUtil.FILE_SEPARATOR + "GAINS" );
3048         }
3049         else if ( state == GainLossStates.LOSS ) {
3050             per_node_go_mapped_domain_gain_loss_files_base_dir = new File( per_node_go_mapped_domain_gain_loss_files_base_dir
3051                     + ForesterUtil.FILE_SEPARATOR + "LOSSES" );
3052         }
3053         else {
3054             per_node_go_mapped_domain_gain_loss_files_base_dir = new File( per_node_go_mapped_domain_gain_loss_files_base_dir
3055                     + ForesterUtil.FILE_SEPARATOR + "PRESENT" );
3056         }
3057         if ( !per_node_go_mapped_domain_gain_loss_files_base_dir.exists() ) {
3058             per_node_go_mapped_domain_gain_loss_files_base_dir.mkdir();
3059         }
3060         return per_node_go_mapped_domain_gain_loss_files_base_dir;
3061     }
3062
3063     private static SortedSet<BinaryDomainCombination> createSetOfAllBinaryDomainCombinationsPerGenome( final GenomeWideCombinableDomains gwcd ) {
3064         final SortedMap<String, CombinableDomains> cds = gwcd.getAllCombinableDomainsIds();
3065         final SortedSet<BinaryDomainCombination> binary_combinations = new TreeSet<BinaryDomainCombination>();
3066         for( final String domain_id : cds.keySet() ) {
3067             final CombinableDomains cd = cds.get( domain_id );
3068             binary_combinations.addAll( cd.toBinaryDomainCombinations() );
3069         }
3070         return binary_combinations;
3071     }
3072
3073     private static void printSomeStats( final DescriptiveStatistics stats, final AsciiHistogram histo, final Writer w )
3074             throws IOException {
3075         w.write( "<hr>" );
3076         w.write( "<br>" );
3077         w.write( SurfacingConstants.NL );
3078         w.write( "<tt><pre>" );
3079         w.write( SurfacingConstants.NL );
3080         if ( histo != null ) {
3081             w.write( histo.toStringBuffer( 20, '|', 40, 5 ).toString() );
3082             w.write( SurfacingConstants.NL );
3083         }
3084         w.write( "</pre></tt>" );
3085         w.write( SurfacingConstants.NL );
3086         w.write( "<table>" );
3087         w.write( SurfacingConstants.NL );
3088         w.write( "<tr><td>N: </td><td>" + stats.getN() + "</td></tr>" );
3089         w.write( SurfacingConstants.NL );
3090         w.write( "<tr><td>Min: </td><td>" + stats.getMin() + "</td></tr>" );
3091         w.write( SurfacingConstants.NL );
3092         w.write( "<tr><td>Max: </td><td>" + stats.getMax() + "</td></tr>" );
3093         w.write( SurfacingConstants.NL );
3094         w.write( "<tr><td>Mean: </td><td>" + stats.arithmeticMean() + "</td></tr>" );
3095         w.write( SurfacingConstants.NL );
3096         if ( stats.getN() > 1 ) {
3097             w.write( "<tr><td>SD: </td><td>" + stats.sampleStandardDeviation() + "</td></tr>" );
3098         }
3099         else {
3100             w.write( "<tr><td>SD: </td><td>n/a</td></tr>" );
3101         }
3102         w.write( SurfacingConstants.NL );
3103         w.write( "</table>" );
3104         w.write( SurfacingConstants.NL );
3105         w.write( "<br>" );
3106         w.write( SurfacingConstants.NL );
3107     }
3108
3109     private static List<String> splitDomainCombination( final String dc ) {
3110         final String[] s = dc.split( "=" );
3111         if ( s.length != 2 ) {
3112             ForesterUtil.printErrorMessage( surfacing.PRG_NAME,
3113                                             "Stringyfied domain combination has illegal format: " + dc );
3114             System.exit( -1 );
3115         }
3116         final List<String> l = new ArrayList<String>( 2 );
3117         l.add( s[ 0 ] );
3118         l.add( s[ 1 ] );
3119         return l;
3120     }
3121
3122     private static void writeAllEncounteredPfamsToFile( final Map<String, List<GoId>> domain_id_to_go_ids_map,
3123                                                         final Map<GoId, GoTerm> go_id_to_term_map,
3124                                                         final String outfile_name,
3125                                                         final SortedSet<String> all_pfams_encountered ) {
3126         final File all_pfams_encountered_file = new File( outfile_name + surfacing.ALL_PFAMS_ENCOUNTERED_SUFFIX );
3127         final File all_pfams_encountered_with_go_annotation_file = new File( outfile_name
3128                 + surfacing.ALL_PFAMS_ENCOUNTERED_WITH_GO_ANNOTATION_SUFFIX );
3129         final File encountered_pfams_summary_file = new File( outfile_name
3130                 + surfacing.ENCOUNTERED_PFAMS_SUMMARY_SUFFIX );
3131         int biological_process_counter = 0;
3132         int cellular_component_counter = 0;
3133         int molecular_function_counter = 0;
3134         int pfams_with_mappings_counter = 0;
3135         int pfams_without_mappings_counter = 0;
3136         int pfams_without_mappings_to_bp_or_mf_counter = 0;
3137         int pfams_with_mappings_to_bp_or_mf_counter = 0;
3138         try {
3139             final Writer all_pfams_encountered_writer = new BufferedWriter( new FileWriter( all_pfams_encountered_file ) );
3140             final Writer all_pfams_encountered_with_go_annotation_writer = new BufferedWriter( new FileWriter( all_pfams_encountered_with_go_annotation_file ) );
3141             final Writer summary_writer = new BufferedWriter( new FileWriter( encountered_pfams_summary_file ) );
3142             summary_writer.write( "# Pfam to GO mapping summary" );
3143             summary_writer.write( ForesterUtil.LINE_SEPARATOR );
3144             summary_writer.write( "# Actual summary is at the end of this file." );
3145             summary_writer.write( ForesterUtil.LINE_SEPARATOR );
3146             summary_writer.write( "# Encountered Pfams without a GO mapping:" );
3147             summary_writer.write( ForesterUtil.LINE_SEPARATOR );
3148             for( final String pfam : all_pfams_encountered ) {
3149                 all_pfams_encountered_writer.write( pfam );
3150                 all_pfams_encountered_writer.write( ForesterUtil.LINE_SEPARATOR );
3151                 final String domain_id = new String( pfam );
3152                 if ( domain_id_to_go_ids_map.containsKey( domain_id ) ) {
3153                     ++pfams_with_mappings_counter;
3154                     all_pfams_encountered_with_go_annotation_writer.write( pfam );
3155                     all_pfams_encountered_with_go_annotation_writer.write( ForesterUtil.LINE_SEPARATOR );
3156                     final List<GoId> go_ids = domain_id_to_go_ids_map.get( domain_id );
3157                     boolean maps_to_bp = false;
3158                     boolean maps_to_cc = false;
3159                     boolean maps_to_mf = false;
3160                     for( final GoId go_id : go_ids ) {
3161                         final GoTerm go_term = go_id_to_term_map.get( go_id );
3162                         if ( go_term.getGoNameSpace().isBiologicalProcess() ) {
3163                             maps_to_bp = true;
3164                         }
3165                         else if ( go_term.getGoNameSpace().isCellularComponent() ) {
3166                             maps_to_cc = true;
3167                         }
3168                         else if ( go_term.getGoNameSpace().isMolecularFunction() ) {
3169                             maps_to_mf = true;
3170                         }
3171                     }
3172                     if ( maps_to_bp ) {
3173                         ++biological_process_counter;
3174                     }
3175                     if ( maps_to_cc ) {
3176                         ++cellular_component_counter;
3177                     }
3178                     if ( maps_to_mf ) {
3179                         ++molecular_function_counter;
3180                     }
3181                     if ( maps_to_bp || maps_to_mf ) {
3182                         ++pfams_with_mappings_to_bp_or_mf_counter;
3183                     }
3184                     else {
3185                         ++pfams_without_mappings_to_bp_or_mf_counter;
3186                     }
3187                 }
3188                 else {
3189                     ++pfams_without_mappings_to_bp_or_mf_counter;
3190                     ++pfams_without_mappings_counter;
3191                     summary_writer.write( pfam );
3192                     summary_writer.write( ForesterUtil.LINE_SEPARATOR );
3193                 }
3194             }
3195             all_pfams_encountered_writer.close();
3196             all_pfams_encountered_with_go_annotation_writer.close();
3197             ForesterUtil.programMessage( surfacing.PRG_NAME,
3198                                          "Wrote all [" + all_pfams_encountered.size() + "] encountered Pfams to: \""
3199                                                  + all_pfams_encountered_file + "\"" );
3200             ForesterUtil.programMessage( surfacing.PRG_NAME,
3201                                          "Wrote all [" + pfams_with_mappings_counter
3202                                                  + "] encountered Pfams with GO mappings to: \""
3203                                                  + all_pfams_encountered_with_go_annotation_file + "\"" );
3204             ForesterUtil.programMessage( surfacing.PRG_NAME,
3205                                          "Wrote summary (including all [" + pfams_without_mappings_counter
3206                                                  + "] encountered Pfams without GO mappings) to: \""
3207                                                  + encountered_pfams_summary_file + "\"" );
3208             ForesterUtil.programMessage( surfacing.PRG_NAME,
3209                                          "Sum of Pfams encountered                : " + all_pfams_encountered.size() );
3210             ForesterUtil.programMessage( surfacing.PRG_NAME,
3211                                          "Pfams without a mapping                 : " + pfams_without_mappings_counter
3212                                                  + " [" + ( ( 100 * pfams_without_mappings_counter )
3213                                                          / all_pfams_encountered.size() )
3214                                                  + "%]" );
3215             ForesterUtil.programMessage( surfacing.PRG_NAME,
3216                                          "Pfams without mapping to proc. or func. : "
3217                                                  + pfams_without_mappings_to_bp_or_mf_counter + " ["
3218                                                  + ( ( 100 * pfams_without_mappings_to_bp_or_mf_counter )
3219                                                          / all_pfams_encountered.size() )
3220                                                  + "%]" );
3221             ForesterUtil
3222                     .programMessage( surfacing.PRG_NAME,
3223                                      "Pfams with a mapping                    : " + pfams_with_mappings_counter + " ["
3224                                              + ( ( 100 * pfams_with_mappings_counter ) / all_pfams_encountered.size() )
3225                                              + "%]" );
3226             ForesterUtil.programMessage( surfacing.PRG_NAME,
3227                                          "Pfams with a mapping to proc. or func.  : "
3228                                                  + pfams_with_mappings_to_bp_or_mf_counter + " ["
3229                                                  + ( ( 100 * pfams_with_mappings_to_bp_or_mf_counter )
3230                                                          / all_pfams_encountered.size() )
3231                                                  + "%]" );
3232             ForesterUtil
3233                     .programMessage( surfacing.PRG_NAME,
3234                                      "Pfams with mapping to biological process: " + biological_process_counter + " ["
3235                                              + ( ( 100 * biological_process_counter ) / all_pfams_encountered.size() )
3236                                              + "%]" );
3237             ForesterUtil
3238                     .programMessage( surfacing.PRG_NAME,
3239                                      "Pfams with mapping to molecular function: " + molecular_function_counter + " ["
3240                                              + ( ( 100 * molecular_function_counter ) / all_pfams_encountered.size() )
3241                                              + "%]" );
3242             ForesterUtil
3243                     .programMessage( surfacing.PRG_NAME,
3244                                      "Pfams with mapping to cellular component: " + cellular_component_counter + " ["
3245                                              + ( ( 100 * cellular_component_counter ) / all_pfams_encountered.size() )
3246                                              + "%]" );
3247             summary_writer.write( ForesterUtil.LINE_SEPARATOR );
3248             summary_writer.write( "# Sum of Pfams encountered                : " + all_pfams_encountered.size() );
3249             summary_writer.write( ForesterUtil.LINE_SEPARATOR );
3250             summary_writer.write( "# Pfams without a mapping                 : " + pfams_without_mappings_counter + " ["
3251                     + ( ( 100 * pfams_without_mappings_counter ) / all_pfams_encountered.size() ) + "%]" );
3252             summary_writer.write( ForesterUtil.LINE_SEPARATOR );
3253             summary_writer.write( "# Pfams without mapping to proc. or func. : "
3254                     + pfams_without_mappings_to_bp_or_mf_counter + " ["
3255                     + ( ( 100 * pfams_without_mappings_to_bp_or_mf_counter ) / all_pfams_encountered.size() ) + "%]" );
3256             summary_writer.write( ForesterUtil.LINE_SEPARATOR );
3257             summary_writer.write( "# Pfams with a mapping                    : " + pfams_with_mappings_counter + " ["
3258                     + ( ( 100 * pfams_with_mappings_counter ) / all_pfams_encountered.size() ) + "%]" );
3259             summary_writer.write( ForesterUtil.LINE_SEPARATOR );
3260             summary_writer.write( "# Pfams with a mapping to proc. or func.  : "
3261                     + pfams_with_mappings_to_bp_or_mf_counter + " ["
3262                     + ( ( 100 * pfams_with_mappings_to_bp_or_mf_counter ) / all_pfams_encountered.size() ) + "%]" );
3263             summary_writer.write( ForesterUtil.LINE_SEPARATOR );
3264             summary_writer.write( "# Pfams with mapping to biological process: " + biological_process_counter + " ["
3265                     + ( ( 100 * biological_process_counter ) / all_pfams_encountered.size() ) + "%]" );
3266             summary_writer.write( ForesterUtil.LINE_SEPARATOR );
3267             summary_writer.write( "# Pfams with mapping to molecular function: " + molecular_function_counter + " ["
3268                     + ( ( 100 * molecular_function_counter ) / all_pfams_encountered.size() ) + "%]" );
3269             summary_writer.write( ForesterUtil.LINE_SEPARATOR );
3270             summary_writer.write( "# Pfams with mapping to cellular component: " + cellular_component_counter + " ["
3271                     + ( ( 100 * cellular_component_counter ) / all_pfams_encountered.size() ) + "%]" );
3272             summary_writer.write( ForesterUtil.LINE_SEPARATOR );
3273             summary_writer.close();
3274         }
3275         catch ( final IOException e ) {
3276             ForesterUtil.printWarningMessage( surfacing.PRG_NAME, "Failure to write: " + e );
3277         }
3278     }
3279
3280     private final static void writeColorLabels( final String l, final Color c, final Writer w ) throws IOException {
3281         w.write( "<tr><td><b><span style=\"color:" );
3282         w.write( String.format( "#%02x%02x%02x", c.getRed(), c.getGreen(), c.getBlue() ) );
3283         w.write( "\">" );
3284         w.write( l );
3285         w.write( "</span></b></td></tr>" );
3286         w.write( SurfacingConstants.NL );
3287     }
3288
3289     private static void writeDomainData( final Map<String, List<GoId>> domain_id_to_go_ids_map,
3290                                          final Map<GoId, GoTerm> go_id_to_term_map,
3291                                          final GoNameSpace go_namespace_limit,
3292                                          final Writer out,
3293                                          final String domain_0,
3294                                          final String domain_1,
3295                                          final String prefix_for_html,
3296                                          final String character_separator_for_non_html_output,
3297                                          final Map<String, Set<String>>[] domain_id_to_secondary_features_maps,
3298                                          final Set<GoId> all_go_ids )
3299             throws IOException {
3300         boolean any_go_annotation_present = false;
3301         boolean first_has_no_go = false;
3302         int domain_count = 2; // To distinguish between domains and binary domain combinations.
3303         if ( ForesterUtil.isEmpty( domain_1 ) ) {
3304             domain_count = 1;
3305         }
3306         // The following has a difficult to understand logic.
3307         for( int d = 0; d < domain_count; ++d ) {
3308             List<GoId> go_ids = null;
3309             boolean go_annotation_present = false;
3310             if ( d == 0 ) {
3311                 if ( domain_id_to_go_ids_map.containsKey( domain_0 ) ) {
3312                     go_annotation_present = true;
3313                     any_go_annotation_present = true;
3314                     go_ids = domain_id_to_go_ids_map.get( domain_0 );
3315                 }
3316                 else {
3317                     first_has_no_go = true;
3318                 }
3319             }
3320             else {
3321                 if ( domain_id_to_go_ids_map.containsKey( domain_1 ) ) {
3322                     go_annotation_present = true;
3323                     any_go_annotation_present = true;
3324                     go_ids = domain_id_to_go_ids_map.get( domain_1 );
3325                 }
3326             }
3327             if ( go_annotation_present ) {
3328                 boolean first = ( ( d == 0 ) || ( ( d == 1 ) && first_has_no_go ) );
3329                 for( final GoId go_id : go_ids ) {
3330                     out.write( "<tr>" );
3331                     if ( first ) {
3332                         first = false;
3333                         writeDomainIdsToHtml( out,
3334                                               domain_0,
3335                                               domain_1,
3336                                               prefix_for_html,
3337                                               domain_id_to_secondary_features_maps );
3338                     }
3339                     else {
3340                         out.write( "<td></td>" );
3341                     }
3342                     if ( !go_id_to_term_map.containsKey( go_id ) ) {
3343                         throw new IllegalArgumentException( "GO-id [" + go_id + "] not found in GO-id to GO-term map" );
3344                     }
3345                     final GoTerm go_term = go_id_to_term_map.get( go_id );
3346                     if ( ( go_namespace_limit == null ) || go_namespace_limit.equals( go_term.getGoNameSpace() ) ) {
3347                         // final String top = GoUtils.getPenultimateGoTerm( go_term, go_id_to_term_map ).getName();
3348                         final String go_id_str = go_id.getId();
3349                         out.write( "<td>" );
3350                         out.write( "<a href=\"" + SurfacingConstants.AMIGO_LINK + go_id_str
3351                                 + "\" target=\"amigo_window\">" + go_id_str + "</a>" );
3352                         out.write( "</td><td>" );
3353                         out.write( go_term.getName() );
3354                         if ( domain_count == 2 ) {
3355                             out.write( " (" + d + ")" );
3356                         }
3357                         out.write( "</td><td>" );
3358                         // out.write( top );
3359                         // out.write( "</td><td>" );
3360                         out.write( "[" );
3361                         out.write( go_term.getGoNameSpace().toShortString() );
3362                         out.write( "]" );
3363                         out.write( "</td>" );
3364                         if ( all_go_ids != null ) {
3365                             all_go_ids.add( go_id );
3366                         }
3367                     }
3368                     else {
3369                         out.write( "<td>" );
3370                         out.write( "</td><td>" );
3371                         out.write( "</td><td>" );
3372                         out.write( "</td><td>" );
3373                         out.write( "</td>" );
3374                     }
3375                     out.write( "</tr>" );
3376                     out.write( SurfacingConstants.NL );
3377                 }
3378             }
3379         } //  for( int d = 0; d < domain_count; ++d )
3380         if ( !any_go_annotation_present ) {
3381             out.write( "<tr>" );
3382             writeDomainIdsToHtml( out, domain_0, domain_1, prefix_for_html, domain_id_to_secondary_features_maps );
3383             out.write( "<td>" );
3384             out.write( "</td><td>" );
3385             out.write( "</td><td>" );
3386             out.write( "</td><td>" );
3387             out.write( "</td>" );
3388             out.write( "</tr>" );
3389             out.write( SurfacingConstants.NL );
3390         }
3391     }
3392
3393     private static void writeDomainIdsToHtml( final Writer out,
3394                                               final String domain_0,
3395                                               final String domain_1,
3396                                               final String prefix_for_detailed_html,
3397                                               final Map<String, Set<String>>[] domain_id_to_secondary_features_maps )
3398             throws IOException {
3399         out.write( "<td>" );
3400         if ( !ForesterUtil.isEmpty( prefix_for_detailed_html ) ) {
3401             out.write( prefix_for_detailed_html );
3402             out.write( " " );
3403         }
3404         out.write( "<a href=\"" + SurfacingConstants.PFAM_FAMILY_ID_LINK + domain_0 + "\">" + domain_0 + "</a>" );
3405         out.write( "</td>" );
3406     }
3407
3408     private static void writeDomainsToIndividualFilePerTreeNode( final Writer individual_files_writer,
3409                                                                  final String domain_0,
3410                                                                  final String domain_1 )
3411             throws IOException {
3412         individual_files_writer.write( domain_0 );
3413         individual_files_writer.write( ForesterUtil.LINE_SEPARATOR );
3414         if ( !ForesterUtil.isEmpty( domain_1 ) ) {
3415             individual_files_writer.write( domain_1 );
3416             individual_files_writer.write( ForesterUtil.LINE_SEPARATOR );
3417         }
3418     }
3419
3420     private static void writePfamsToFile( final String outfile_name, final SortedSet<String> pfams ) {
3421         try {
3422             final Writer writer = new BufferedWriter( new FileWriter( new File( outfile_name ) ) );
3423             for( final String pfam : pfams ) {
3424                 writer.write( pfam );
3425                 writer.write( ForesterUtil.LINE_SEPARATOR );
3426             }
3427             writer.close();
3428             ForesterUtil.programMessage( surfacing.PRG_NAME,
3429                                          "Wrote " + pfams.size() + " pfams to [" + outfile_name + "]" );
3430         }
3431         catch ( final IOException e ) {
3432             ForesterUtil.printWarningMessage( surfacing.PRG_NAME, "Failure to write: " + e );
3433         }
3434     }
3435
3436     private static void writeToNexus( final String outfile_name,
3437                                       final CharacterStateMatrix<BinaryStates> matrix,
3438                                       final Phylogeny phylogeny ) {
3439         if ( !( matrix instanceof BasicCharacterStateMatrix ) ) {
3440             throw new IllegalArgumentException( "can only write matrices of type [" + BasicCharacterStateMatrix.class
3441                     + "] to nexus" );
3442         }
3443         final BasicCharacterStateMatrix<BinaryStates> my_matrix = ( org.forester.evoinference.matrix.character.BasicCharacterStateMatrix<BinaryStates> ) matrix;
3444         final List<Phylogeny> phylogenies = new ArrayList<Phylogeny>( 1 );
3445         phylogenies.add( phylogeny );
3446         try {
3447             final BufferedWriter w = new BufferedWriter( new FileWriter( outfile_name ) );
3448             w.write( NexusConstants.NEXUS );
3449             w.write( ForesterUtil.LINE_SEPARATOR );
3450             my_matrix.writeNexusTaxaBlock( w );
3451             my_matrix.writeNexusBinaryChractersBlock( w );
3452             PhylogenyWriter.writeNexusTreesBlock( w, phylogenies, NH_CONVERSION_SUPPORT_VALUE_STYLE.NONE );
3453             w.flush();
3454             w.close();
3455             ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote Nexus file: \"" + outfile_name + "\"" );
3456         }
3457         catch ( final IOException e ) {
3458             ForesterUtil.fatalError( surfacing.PRG_NAME, e.getMessage() );
3459         }
3460     }
3461
3462     private static void writeToNexus( final String outfile_name,
3463                                       final DomainParsimonyCalculator domain_parsimony,
3464                                       final Phylogeny phylogeny ) {
3465         writeToNexus( outfile_name + surfacing.NEXUS_EXTERNAL_DOMAINS,
3466                       domain_parsimony.createMatrixOfDomainPresenceOrAbsence(),
3467                       phylogeny );
3468         writeToNexus( outfile_name + surfacing.NEXUS_EXTERNAL_DOMAIN_COMBINATIONS,
3469                       domain_parsimony.createMatrixOfBinaryDomainCombinationPresenceOrAbsence(),
3470                       phylogeny );
3471     }
3472
3473     final static class DomainComparator implements Comparator<Domain> {
3474
3475         final private boolean _ascending;
3476
3477         public DomainComparator( final boolean ascending ) {
3478             _ascending = ascending;
3479         }
3480
3481         @Override
3482         public final int compare( final Domain d0, final Domain d1 ) {
3483             if ( d0.getFrom() < d1.getFrom() ) {
3484                 return _ascending ? -1 : 1;
3485             }
3486             else if ( d0.getFrom() > d1.getFrom() ) {
3487                 return _ascending ? 1 : -1;
3488             }
3489             return 0;
3490         }
3491     }
3492 }