removed unneeded import
[jalview.git] / forester / java / src / org / forester / surfacing / MinimalDomainomeCalculator.java
index 9abb02a..0c9ca73 100644 (file)
@@ -19,6 +19,7 @@ import java.util.TreeSet;
 
 import org.forester.application.surfacing;
 import org.forester.phylogeny.Phylogeny;
+import org.forester.phylogeny.PhylogenyMethods;
 import org.forester.phylogeny.PhylogenyNode;
 import org.forester.phylogeny.iterators.PhylogenyNodeIterator;
 import org.forester.protein.Domain;
@@ -30,12 +31,14 @@ import org.forester.util.ForesterUtil;
 
 public final class MinimalDomainomeCalculator {
 
-    static final public void calcOme( final boolean use_domain_architectures,
-                                      final Phylogeny tre,
-                                      final SortedMap<Species, List<Protein>> protein_lists_per_species,
-                                      final String separator,
-                                      final double ie_cutoff,
-                                      final String outfile_base )
+    public final static void calc( final boolean use_domain_architectures,
+                                   final Phylogeny tre,
+                                   final int target_level,
+                                   final SortedMap<Species, List<Protein>> protein_lists_per_species,
+                                   final String separator,
+                                   final double ie_cutoff,
+                                   final String outfile_base,
+                                   final boolean write_protein_files )
             throws IOException {
         final SortedMap<String, SortedSet<String>> species_to_features_map = new TreeMap<String, SortedSet<String>>();
         if ( protein_lists_per_species == null || tre == null ) {
@@ -59,8 +62,20 @@ public final class MinimalDomainomeCalculator {
         final BufferedWriter out_table = new BufferedWriter( new FileWriter( outfile_table ) );
         out.write( "SPECIES\tCOMMON NAME\tCODE\tRANK\t#EXT NODES\tEXT NODE CODES\t#" + x + "\t" + x + "" );
         out.write( ForesterUtil.LINE_SEPARATOR );
+        ///////////
+        //////////
+        SortedMap<String, List<Protein>> protein_lists_per_quasi_species = null;
+        if ( target_level >= 1 ) {
+            protein_lists_per_quasi_species = makeProteinListsPerQuasiSpecies( tre,
+                                                                               target_level,
+                                                                               protein_lists_per_species );
+           
+        }
+        /////////
+        ///////////
         for( final PhylogenyNodeIterator iter = tre.iteratorPostorder(); iter.hasNext(); ) {
             final PhylogenyNode node = iter.next();
+            final int node_level = PhylogenyMethods.calculateLevel( node );
             final String species_name = node.getNodeData().isHasTaxonomy()
                     ? node.getNodeData().getTaxonomy().getScientificName() : node.getName();
             final String common = node.getNodeData().isHasTaxonomy() ? node.getNodeData().getTaxonomy().getCommonName()
@@ -68,67 +83,126 @@ public final class MinimalDomainomeCalculator {
             final String tcode = node.getNodeData().isHasTaxonomy() ? node.getNodeData().getTaxonomy().getTaxonomyCode()
                     : "";
             final String rank = node.getNodeData().isHasTaxonomy() ? node.getNodeData().getTaxonomy().getRank() : "";
-            out.write( species_name );
-            if ( !ForesterUtil.isEmpty( common ) ) {
-                out.write( "\t" + common );
-            }
-            else {
-                out.write( "\t" );
-            }
-            if ( !ForesterUtil.isEmpty( tcode ) ) {
-                out.write( "\t" + tcode );
-            }
-            else {
-                out.write( "\t" );
-            }
-            if ( !ForesterUtil.isEmpty( rank ) ) {
-                out.write( "\t" + rank );
-            }
-            else {
-                out.write( "\t" );
-            }
             final List<PhylogenyNode> external_descs = node.getAllExternalDescendants();
-            if ( node.isInternal() ) {
-                out.write( "\t" + external_descs.size() + "\t" );
-            }
-            else {
-                out.write( "\t\t" );
+            if ( ( target_level < 1 ) || ( node_level >= target_level ) ) {
+                out.write( species_name );
+                if ( !ForesterUtil.isEmpty( common ) ) {
+                    out.write( "\t" + common );
+                }
+                else {
+                    out.write( "\t" );
+                }
+                if ( !ForesterUtil.isEmpty( tcode ) ) {
+                    out.write( "\t" + tcode );
+                }
+                else {
+                    out.write( "\t" );
+                }
+                if ( !ForesterUtil.isEmpty( rank ) ) {
+                    out.write( "\t" + rank );
+                }
+                else {
+                    out.write( "\t" );
+                }
+                if ( node.isInternal() ) {
+                    out.write( "\t" + external_descs.size() + "\t" );
+                }
+                else {
+                    out.write( "\t\t" );
+                }
             }
             final List<Set<String>> features_per_genome_list = new ArrayList<Set<String>>();
             boolean first = true;
-            for( final PhylogenyNode external_desc : external_descs ) {
-                final String code = external_desc.getNodeData().getTaxonomy().getTaxonomyCode();
-                if ( node.isInternal() ) {
-                    if ( first ) {
-                        first = false;
-                    }
-                    else {
-                        out.write( ", " );
+            if ( target_level >= 1 ) {
+                ////////////
+                ////////////
+                if ( node_level >= target_level ) {
+                    final List<PhylogenyNode> given_level_descs = PhylogenyMethods
+                            .getAllDescendantsOfGivenLevel( node, target_level );
+                    for( final PhylogenyNode given_level_desc : given_level_descs ) {
+                        final String spec_name = given_level_desc.getNodeData().isHasTaxonomy()
+                                ? given_level_desc.getNodeData().getTaxonomy().getScientificName()
+                                : given_level_desc.getName();
+                        if ( node.isInternal() ) {
+                            if ( first ) {
+                                first = false;
+                            }
+                            else {
+                                out.write( ", " );
+                            }
+                            out.write( "sp_n=" + spec_name );
+                        }
+                        final List<Protein> proteins_per_species = protein_lists_per_quasi_species.get( spec_name );
+                        if ( proteins_per_species != null ) {
+                            final SortedSet<String> features_per_genome = new TreeSet<String>();
+                            for( final Protein protein : proteins_per_species ) {
+                                if ( use_domain_architectures ) {
+                                    final String da = protein.toDomainArchitectureString( separator, ie_cutoff );
+                                    features_per_genome.add( da );
+                                }
+                                else {
+                                    List<Domain> domains = protein.getProteinDomains();
+                                    for( final Domain domain : domains ) {
+                                        if ( ( ie_cutoff <= -1 ) || ( domain.getPerDomainEvalue() <= ie_cutoff ) ) {
+                                            features_per_genome.add( domain.getDomainId() );
+                                        }
+                                    }
+                                }
+                            }
+                            System.out.println( ">>>>>>>>>>>>>> features_per_genome.size()=" + features_per_genome.size() );
+                            if ( features_per_genome.size() > 0 ) {
+                                features_per_genome_list.add( features_per_genome );
+                            }
+                            else {
+                                System.out.println( "error!" );
+                                System.exit( -1 );
+                            }
+                        }
+                        else {
+                            System.out.println( "error!" );
+                            System.exit( -1 );
+                        }
                     }
-                    out.write( code );
                 }
-                final List<Protein> proteins_per_species = protein_lists_per_species.get( new BasicSpecies( code ) );
-                if ( proteins_per_species != null ) {
-                    final SortedSet<String> features_per_genome = new TreeSet<String>();
-                    for( final Protein protein : proteins_per_species ) {
-                        if ( use_domain_architectures ) {
-                            final String da = protein.toDomainArchitectureString( separator, ie_cutoff );
-                            features_per_genome.add( da );
+                ///////////
+                ///////////
+            }
+            else {
+                for( final PhylogenyNode external_desc : external_descs ) {
+                    final String code = external_desc.getNodeData().getTaxonomy().getTaxonomyCode();
+                    if ( node.isInternal() ) {
+                        if ( first ) {
+                            first = false;
                         }
                         else {
-                            List<Domain> domains = protein.getProteinDomains();
-                            for( final Domain domain : domains ) {
-                                if ( ( ie_cutoff <= -1 ) || ( domain.getPerDomainEvalue() <= ie_cutoff ) ) {
-                                    features_per_genome.add( domain.getDomainId() );
+                            out.write( ", " );
+                        }
+                        out.write( code );
+                    }
+                    final List<Protein> proteins_per_species = protein_lists_per_species
+                            .get( new BasicSpecies( code ) );
+                    if ( proteins_per_species != null ) {
+                        final SortedSet<String> features_per_genome = new TreeSet<String>();
+                        for( final Protein protein : proteins_per_species ) {
+                            if ( use_domain_architectures ) {
+                                final String da = protein.toDomainArchitectureString( separator, ie_cutoff );
+                                features_per_genome.add( da );
+                            }
+                            else {
+                                List<Domain> domains = protein.getProteinDomains();
+                                for( final Domain domain : domains ) {
+                                    if ( ( ie_cutoff <= -1 ) || ( domain.getPerDomainEvalue() <= ie_cutoff ) ) {
+                                        features_per_genome.add( domain.getDomainId() );
+                                    }
                                 }
                             }
                         }
+                        if ( features_per_genome.size() > 0 ) {
+                            features_per_genome_list.add( features_per_genome );
+                        }
                     }
-                    if ( features_per_genome.size() > 0 ) {
-                        features_per_genome_list.add( features_per_genome );
-                    }
-                }
-            }
+                } // for( final PhylogenyNode external_desc : external_descs )
+            } // else
             if ( features_per_genome_list.size() > 0 ) {
                 SortedSet<String> intersection = calcIntersection( features_per_genome_list );
                 out.write( "\t" + intersection.size() + "\t" );
@@ -192,26 +266,80 @@ public final class MinimalDomainomeCalculator {
         out_table.close();
         ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote minimal DAome data to           : " + outfile );
         ForesterUtil.programMessage( surfacing.PRG_NAME, "Wrote minimal DAome data to (as table): " + outfile_table );
-        for( String f : all_features ) {
+        if ( write_protein_files ) {
+            final String protdirname;
             final String a;
+            final String b;
             if ( use_domain_architectures ) {
-                a = "DA_";
+                a = "_DA";
+                b = "domain architectures (DAs)";
+                protdirname = "_DAS";
             }
             else {
-                a = "domain_";
+                a = "_domain";
+                b = "domains";
+                protdirname = "_DOMAINS";
+            }
+            final File prot_dir = new File( outfile_base + protdirname );
+            final boolean success = prot_dir.mkdir();
+            if ( !success ) {
+                throw new IOException( "failed to create dir " + prot_dir );
+            }
+            int total = 0;
+            final String dir = outfile_base + protdirname + "/";
+            for( final String feat : all_features ) {
+                final File extract_outfile = new File( dir + feat + a + surfacing.SEQ_EXTRACT_SUFFIX );
+                SurfacingUtil.checkForOutputFileWriteability( extract_outfile );
+                final Writer proteins_file_writer = new BufferedWriter( new FileWriter( extract_outfile ) );
+                final int counter = extractProteinFeatures( use_domain_architectures,
+                                                            protein_lists_per_species,
+                                                            feat,
+                                                            proteins_file_writer,
+                                                            ie_cutoff,
+                                                            separator );
+                if ( counter < 1 ) {
+                    ForesterUtil.printWarningMessage( "surfacing", feat + " not present (in " + b + " extraction)" );
+                }
+                total += counter;
+                proteins_file_writer.close();
+            }
+            ForesterUtil.programMessage( "surfacing",
+                                         "Wrote " + total + " individual " + b + " from a total of "
+                                                 + all_features.size() + " into: " + dir );
+        }
+    }
+
+    private final static SortedMap<String, List<Protein>> makeProteinListsPerQuasiSpecies( final Phylogeny tre,
+                                                                                           final int level,
+                                                                                           final SortedMap<Species, List<Protein>> protein_lists_per_species ) {
+        final SortedMap<String, List<Protein>> protein_lists_per_quasi_species = new TreeMap<String, List<Protein>>();
+        System.out.println( "---------------------------------" );
+        System.out.println( "level=" + level );
+        for( final PhylogenyNodeIterator iter = tre.iteratorPostorder(); iter.hasNext(); ) {
+            final PhylogenyNode node = iter.next();
+            final int node_level = PhylogenyMethods.calculateLevel( node );
+            if ( node_level == level ) {
+                System.out.println( "level=" + level );
+                final List<PhylogenyNode> external_descs = node.getAllExternalDescendants();
+                final List<Protein> protein_list_per_quasi_species = new ArrayList<Protein>();
+                for( final PhylogenyNode external_desc : external_descs ) {
+                    final String code = external_desc.getNodeData().getTaxonomy().getTaxonomyCode();
+                    final List<Protein> proteins_per_species = protein_lists_per_species
+                            .get( new BasicSpecies( code ) );
+                    //System.out.println( code );
+                    for( Protein protein : proteins_per_species ) {
+                        protein_list_per_quasi_species.add( protein );
+                    }
+                }
+                final String species_name = node.getNodeData().isHasTaxonomy()
+                        ? node.getNodeData().getTaxonomy().getScientificName() : node.getName();
+                System.out.println( "species_name=" + species_name );
+                protein_lists_per_quasi_species.put( species_name, protein_list_per_quasi_species );
+                System.out.println( ">>>>" + protein_list_per_quasi_species.size() );
             }
-            final File prot_dir = new File( outfile_base + "_prot" );
-            prot_dir.mkdir();
-            final File outt = new File( outfile_base + "_prot/" + a + f + surfacing.SEQ_EXTRACT_SUFFIX );
-            final Writer proteins_file_writer = new BufferedWriter( new FileWriter( outt ) );
-            extractProteinFeatures( use_domain_architectures,
-                                    protein_lists_per_species,
-                                    f,
-                                    proteins_file_writer,
-                                    ie_cutoff,
-                                    separator );
-            proteins_file_writer.close();
         }
+      
+        return protein_lists_per_quasi_species;
     }
 
     private final static SortedSet<String> calcIntersection( final List<Set<String>> features_per_genome_list ) {
@@ -226,13 +354,14 @@ public final class MinimalDomainomeCalculator {
         return my_first;
     }
 
-    public static void extractProteinFeatures( final boolean use_domain_architectures,
-                                               final SortedMap<Species, List<Protein>> protein_lists_per_species,
-                                               final String domain_id,
-                                               final Writer out,
-                                               final double ie_cutoff,
-                                               final String domain_separator )
+    private final static int extractProteinFeatures( final boolean use_domain_architectures,
+                                                     final SortedMap<Species, List<Protein>> protein_lists_per_species,
+                                                     final String domain_id,
+                                                     final Writer out,
+                                                     final double ie_cutoff,
+                                                     final String domain_separator )
             throws IOException {
+        int counter = 0;
         final String separator_for_output = "\t";
         for( final Species species : protein_lists_per_species.keySet() ) {
             final List<Protein> proteins_per_species = protein_lists_per_species.get( species );
@@ -261,6 +390,7 @@ public final class MinimalDomainomeCalculator {
                         out.write( from + "-" + to );
                         out.write( "/" );
                         out.write( SurfacingConstants.NL );
+                        ++counter;
                     }
                 }
                 else {
@@ -315,11 +445,13 @@ public final class MinimalDomainomeCalculator {
                             out.write( protein.getAccession() );
                         }
                         out.write( SurfacingConstants.NL );
+                        ++counter;
                     }
                 }
             }
         }
         out.flush();
+        return counter;
     }
 
     public static void main( final String[] args ) {