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