forester/java/src/org/forester/msa/MsaMethods.java

   1 // $Id:
   2 // forester -- software libraries and applications
   3 // for genomics and evolutionary biology research.
   4 //
   5 // Copyright (C) 2010 Christian M Zmasek
   6 // Copyright (C) 2010 Sanford-Burnham Medical Research Institute
   7 // All rights reserved
   8 //
   9 // This library is free software; you can redistribute it and/or
  10 // modify it under the terms of the GNU Lesser General Public
  11 // License as published by the Free Software Foundation; either
  12 // version 2.1 of the License, or (at your option) any later version.
  13 //
  14 // This library is distributed in the hope that it will be useful,
  15 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  17 // Lesser General Public License for more details.
  18 //
  19 // You should have received a copy of the GNU Lesser General Public
  20 // License along with this library; if not, write to the Free Software
  21 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
  22 //
  23 // Contact: phylosoft @ gmail . com
  24 // WWW: https://sites.google.com/site/cmzmasek/home/software/forester
  25
  26 package org.forester.msa;
  27
  28 import java.util.ArrayList;
  29 import java.util.Iterator;
  30 import java.util.List;
  31 import java.util.Map;
  32 import java.util.SortedMap;
  33 import java.util.TreeMap;
  34
  35 import org.forester.sequence.BasicSequence;
  36 import org.forester.sequence.Sequence;
  37 import org.forester.util.BasicDescriptiveStatistics;
  38 import org.forester.util.DescriptiveStatistics;
  39
  40 public final class MsaMethods {
  41
  42     private ArrayList<String> _ignored_seqs_ids;
  43
  44     synchronized public ArrayList<String> getIgnoredSequenceIds() {
  45         return _ignored_seqs_ids;
  46     }
  47
  48     synchronized public static MsaMethods createInstance() {
  49         return new MsaMethods();
  50     }
  51
  52     private MsaMethods() {
  53         init();
  54     }
  55
  56     synchronized private void init() {
  57         _ignored_seqs_ids = new ArrayList<String>();
  58     }
  59
  60     @Override
  61     public Object clone() {
  62         throw new NoSuchMethodError();
  63     }
  64
  65     public static int calcGapSumPerColumn( final Msa msa, final int col ) {
  66         int gap_rows = 0;
  67         for( int j = 0; j < msa.getNumberOfSequences(); ++j ) {
  68             if ( msa.isGapAt( j, col ) ) {
  69                 gap_rows++;
  70             }
  71         }
  72         return gap_rows;
  73     }
  74
  75     final public static Msa removeSequence( final Msa msa, final String to_remove_id ) {
  76         final List<Sequence> seqs = new ArrayList<Sequence>();
  77         for( int row = 0; row < msa.getNumberOfSequences(); ++row ) {
  78             if ( !to_remove_id.equals( msa.getIdentifier( row ) ) ) {
  79                 seqs.add( msa.getSequence( row ) );
  80             }
  81         }
  82         if ( seqs.size() < 1 ) {
  83             return null;
  84         }
  85         return BasicMsa.createInstance( seqs );
  86     }
  87
  88     final public static Msa removeSequences( final Msa msa, final List<String> to_remove_ids ) {
  89         final List<Sequence> seqs = new ArrayList<Sequence>();
  90         for( int row = 0; row < msa.getNumberOfSequences(); ++row ) {
  91             if ( !to_remove_ids.contains( msa.getIdentifier( row ) ) ) {
  92                 seqs.add( msa.getSequence( row ) );
  93             }
  94         }
  95         if ( seqs.size() < 1 ) {
  96             return null;
  97         }
  98         return BasicMsa.createInstance( seqs );
  99     }
 100
 101     final public static Msa removeSequencesByRow( final Msa msa, final List<Integer> to_remove_rows ) {
 102         final List<Sequence> seqs = new ArrayList<Sequence>();
 103         for( int row = 0; row < msa.getNumberOfSequences(); ++row ) {
 104             if ( !to_remove_rows.contains( row ) ) {
 105                 seqs.add( msa.getSequence( row ) );
 106             }
 107         }
 108         if ( seqs.size() < 1 ) {
 109             return null;
 110         }
 111         return BasicMsa.createInstance( seqs );
 112     }
 113
 114     synchronized final public Msa deleteGapColumns( final double max_allowed_gap_ratio,
 115                                                     final int min_allowed_length,
 116                                                     final Msa msa ) {
 117         init();
 118         if ( ( max_allowed_gap_ratio < 0 ) || ( max_allowed_gap_ratio > 1 ) ) {
 119             throw new IllegalArgumentException( "max allowed gap ration is out of range: " + max_allowed_gap_ratio );
 120         }
 121         final boolean ignore_too_short_seqs = min_allowed_length > 0;
 122         final boolean[] delete_cols = new boolean[ msa.getLength() ];
 123         int new_length = 0;
 124         for( int col = 0; col < msa.getLength(); ++col ) {
 125             delete_cols[ col ] = ( ( double ) calcGapSumPerColumn( msa, col ) / msa.getNumberOfSequences() ) > max_allowed_gap_ratio;
 126             if ( !delete_cols[ col ] ) {
 127                 ++new_length;
 128             }
 129         }
 130         final List<Sequence> seqs = new ArrayList<Sequence>( msa.getNumberOfSequences() );
 131         for( int row = 0; row < msa.getNumberOfSequences(); ++row ) {
 132             final char[] mol_seq = new char[ new_length ];
 133             int new_col = 0;
 134             int non_gap_cols_sum = 0;
 135             for( int col = 0; col < msa.getLength(); ++col ) {
 136                 if ( !delete_cols[ col ] ) {
 137                     final char residue = msa.getResidueAt( row, col );
 138                     mol_seq[ new_col++ ] = ( residue );
 139                     if ( residue != Sequence.GAP ) {
 140                         ++non_gap_cols_sum;
 141                     }
 142                 }
 143             }
 144             if ( ignore_too_short_seqs ) {
 145                 if ( non_gap_cols_sum >= min_allowed_length ) {
 146                     seqs.add( new BasicSequence( msa.getIdentifier( row ), mol_seq, msa.getType() ) );
 147                 }
 148                 else {
 149                     _ignored_seqs_ids.add( msa.getIdentifier( row ).toString() );
 150                 }
 151             }
 152             else {
 153                 seqs.add( new BasicSequence( msa.getIdentifier( row ), mol_seq, msa.getType() ) );
 154             }
 155         }
 156         if ( seqs.size() < 1 ) {
 157             return null;
 158         }
 159         return BasicMsa.createInstance( seqs );
 160     }
 161
 162     final public static DescriptiveStatistics calculateIdentityRatio( final int from, final int to, final Msa msa ) {
 163         final DescriptiveStatistics stats = new BasicDescriptiveStatistics();
 164         for( int c = from; c <= to; ++c ) {
 165             stats.addValue( calculateIdentityRatio( msa, c ) );
 166         }
 167         return stats;
 168     }
 169
 170     final public static DescriptiveStatistics calculateEffectiveLengthStatistics( final Msa msa ) {
 171         final DescriptiveStatistics stats = new BasicDescriptiveStatistics();
 172         for( int row = 0; row < msa.getNumberOfSequences(); ++row ) {
 173             final Sequence s = msa.getSequence( row );
 174             stats.addValue( s.getLength() - s.getNumberOfGapResidues() );
 175         }
 176         return stats;
 177     }
 178
 179     final public static double calculateIdentityRatio( final Msa msa, final int column ) {
 180         final SortedMap<Character, Integer> dist = calculateResidueDestributionPerColumn( msa, column );
 181         int majority_count = 0;
 182         final Iterator<Map.Entry<Character, Integer>> it = dist.entrySet().iterator();
 183         while ( it.hasNext() ) {
 184             final Map.Entry<Character, Integer> pair = it.next();
 185             if ( pair.getValue() > majority_count ) {
 186                 majority_count = pair.getValue();
 187             }
 188         }
 189         return ( double ) majority_count / msa.getNumberOfSequences();
 190     }
 191
 192     public static SortedMap<Character, Integer> calculateResidueDestributionPerColumn( final Msa msa, final int column ) {
 193         final SortedMap<Character, Integer> map = new TreeMap<Character, Integer>();
 194         for( final Character r : msa.getColumnAt( column ) ) {
 195             if ( r != Sequence.GAP ) {
 196                 if ( !map.containsKey( r ) ) {
 197                     map.put( r, 1 );
 198                 }
 199                 else {
 200                     map.put( r, map.get( r ) + 1 );
 201                 }
 202             }
 203         }
 204         return map;
 205     }
 206
 207     public static DescriptiveStatistics calcBasicGapinessStatistics( final Msa msa ) {
 208         final DescriptiveStatistics stats = new BasicDescriptiveStatistics();
 209         for( int i = 0; i < msa.getLength(); ++i ) {
 210             stats.addValue( ( double ) calcGapSumPerColumn( msa, i ) / msa.getNumberOfSequences() );
 211         }
 212         return stats;
 213     }
 214
 215     public static Msa removeSequencesByMinimalLength( final Msa msa, final int min_effective_length ) {
 216         final List<Integer> to_remove_rows = new ArrayList<Integer>();
 217         for( int seq = 0; seq < msa.getNumberOfSequences(); ++seq ) {
 218             int eff_length = 0;
 219             for( int i = 0; i < msa.getLength(); ++i ) {
 220                 if ( msa.getResidueAt( seq, i ) != Sequence.GAP ) {
 221                     eff_length++;
 222                 }
 223             }
 224             if ( eff_length < min_effective_length ) {
 225                 to_remove_rows.add( seq );
 226             }
 227         }
 228         return removeSequencesByRow( msa, to_remove_rows );
 229     }
 230
 231     public static double calcGapRatio( final Msa msa ) {
 232         int gaps = 0;
 233         for( int seq = 0; seq < msa.getNumberOfSequences(); ++seq ) {
 234             for( int i = 0; i < msa.getLength(); ++i ) {
 235                 if ( msa.getResidueAt( seq, i ) == Sequence.GAP ) {
 236                     gaps++;
 237                 }
 238             }
 239         }
 240         return ( double ) gaps / ( msa.getLength() * msa.getNumberOfSequences() );
 241     }
 242 }