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.HashMap;
  30 import java.util.Iterator;
  31 import java.util.List;
  32 import java.util.Map;
  33 import java.util.SortedMap;
  34 import java.util.TreeMap;
  35
  36 import org.forester.sequence.BasicSequence;
  37 import org.forester.sequence.MolecularSequence;
  38 import org.forester.util.BasicDescriptiveStatistics;
  39 import org.forester.util.DescriptiveStatistics;
  40
  41 public final class MsaMethods {
  42
  43     private ArrayList<String> _ignored_seqs_ids;
  44
  45     private MsaMethods() {
  46         init();
  47     }
  48
  49     @Override
  50     public Object clone() {
  51         throw new NoSuchMethodError();
  52     }
  53
  54     synchronized final public Msa deleteGapColumns( final double max_allowed_gap_ratio,
  55                                                     final int min_allowed_length,
  56                                                     final Msa msa ) {
  57         init();
  58         if ( ( max_allowed_gap_ratio < 0 ) || ( max_allowed_gap_ratio > 1 ) ) {
  59             throw new IllegalArgumentException( "max allowed gap ration is out of range: " + max_allowed_gap_ratio );
  60         }
  61         final boolean ignore_too_short_seqs = min_allowed_length > 0;
  62         final boolean[] delete_cols = new boolean[ msa.getLength() ];
  63         int new_length = 0;
  64         for( int col = 0; col < msa.getLength(); ++col ) {
  65             delete_cols[ col ] = ( ( double ) calcGapSumPerColumn( msa, col ) / msa.getNumberOfSequences() ) > max_allowed_gap_ratio;
  66             if ( !delete_cols[ col ] ) {
  67                 ++new_length;
  68             }
  69         }
  70         final List<MolecularSequence> seqs = new ArrayList<MolecularSequence>( msa.getNumberOfSequences() );
  71         for( int row = 0; row < msa.getNumberOfSequences(); ++row ) {
  72             final char[] mol_seq = new char[ new_length ];
  73             int new_col = 0;
  74             int non_gap_cols_sum = 0;
  75             for( int col = 0; col < msa.getLength(); ++col ) {
  76                 if ( !delete_cols[ col ] ) {
  77                     final char residue = msa.getResidueAt( row, col );
  78                     mol_seq[ new_col++ ] = ( residue );
  79                     if ( residue != MolecularSequence.GAP ) {
  80                         ++non_gap_cols_sum;
  81                     }
  82                 }
  83             }
  84             if ( ignore_too_short_seqs ) {
  85                 if ( non_gap_cols_sum >= min_allowed_length ) {
  86                     seqs.add( new BasicSequence( msa.getIdentifier( row ), mol_seq, msa.getType() ) );
  87                 }
  88                 else {
  89                     _ignored_seqs_ids.add( msa.getIdentifier( row ).toString() );
  90                 }
  91             }
  92             else {
  93                 seqs.add( new BasicSequence( msa.getIdentifier( row ), mol_seq, msa.getType() ) );
  94             }
  95         }
  96         if ( seqs.size() < 1 ) {
  97             return null;
  98         }
  99         return BasicMsa.createInstance( seqs );
 100     }
 101
 102     synchronized public ArrayList<String> getIgnoredSequenceIds() {
 103         return _ignored_seqs_ids;
 104     }
 105
 106     synchronized private void init() {
 107         _ignored_seqs_ids = new ArrayList<String>();
 108     }
 109
 110     public static final DescriptiveStatistics calcNumberOfGapsStats( final Msa msa ) {
 111         final int[] gaps = calcNumberOfGapsInMsa( msa );
 112         final DescriptiveStatistics stats = new BasicDescriptiveStatistics();
 113         for( final int gap : gaps ) {
 114             stats.addValue(  gap );
 115         }
 116         return stats;
 117     }
 118
 119     public static final int[] calcNumberOfGapsInMsa( final Msa msa ) {
 120         final int seqs = msa.getNumberOfSequences();
 121         final int[]  gaps= new int[ seqs ];
 122         for( int i = 0; i < seqs; ++i ) {
 123             gaps[ i ] =  calcNumberOfGaps( msa.getSequence( i ) );
 124         }
 125         return gaps;
 126     }
 127
 128
 129
 130     public final static int calcNumberOfGaps( final MolecularSequence seq  ) {
 131         int gaps = 0;
 132         boolean was_gap = false;
 133         for( int i = 0; i < seq.getLength(); ++i ) {
 134             if ( seq.isGapAt( i ) ) {
 135                if ( !was_gap ) {
 136                    ++gaps;
 137                    was_gap = true;
 138                }
 139             }
 140             else {
 141                 was_gap = false;
 142             }
 143         }
 144         return gaps;
 145     }
 146
 147     public static DescriptiveStatistics calcBasicGapinessStatistics( final Msa msa ) {
 148         final DescriptiveStatistics stats = new BasicDescriptiveStatistics();
 149         for( int i = 0; i < msa.getLength(); ++i ) {
 150             stats.addValue( ( double ) calcGapSumPerColumn( msa, i ) / msa.getNumberOfSequences() );
 151         }
 152         return stats;
 153     }
 154
 155     public static double calcGapRatio( final Msa msa ) {
 156         int gaps = 0;
 157         for( int seq = 0; seq < msa.getNumberOfSequences(); ++seq ) {
 158             for( int i = 0; i < msa.getLength(); ++i ) {
 159                 if ( msa.getResidueAt( seq, i ) == MolecularSequence.GAP ) {
 160                     gaps++;
 161                 }
 162             }
 163         }
 164         return ( double ) gaps / ( msa.getLength() * msa.getNumberOfSequences() );
 165     }
 166
 167     public static int calcGapSumPerColumn( final Msa msa, final int col ) {
 168         int gap_rows = 0;
 169         for( int j = 0; j < msa.getNumberOfSequences(); ++j ) {
 170             if ( msa.isGapAt( j, col ) ) {
 171                 gap_rows++;
 172             }
 173         }
 174         return gap_rows;
 175     }
 176
 177     final public static double calcNormalizedShannonsEntropy( final int k, final Msa msa ) {
 178         double s = 0;
 179         for( int col = 0; col < msa.getLength(); ++col ) {
 180             s += calcNormalizedShannonsEntropy( k, msa, col );
 181         }
 182         return s / msa.getLength();
 183     }
 184
 185     final public static double calcNormalizedShannonsEntropy( final int k, final Msa msa, final int col ) {
 186         // http://www.ebi.ac.uk/thornton-srv/databases/valdarprograms/scorecons_server_help.html
 187         // n: number of residues in column
 188         // k: number of residue types
 189         // na: number of residues of type a
 190         // pa = na/n
 191         // S=-sum pa log2 pa
 192         double s = 0;
 193         final double n = msa.getNumberOfSequences();
 194         HashMap<Character, Integer> dist = null;
 195         if ( k == 6 ) {
 196             dist = calcResidueDistribution6( msa, col );
 197         }
 198         else if ( k == 7 ) {
 199             dist = calcResidueDistribution7( msa, col );
 200         }
 201         else if ( k == 20 ) {
 202             dist = calcResidueDistribution20( msa, col );
 203         }
 204         else if ( k == 21 ) {
 205             dist = calcResidueDistribution21( msa, col );
 206         }
 207         else {
 208             throw new IllegalArgumentException( "illegal value for k: " + k );
 209         }
 210         if ( dist.size() == 1 ) {
 211             return 0;
 212         }
 213         //        if ( dist.size() == n ) {
 214         //            return 0;
 215         //        }
 216         for( final int na : dist.values() ) {
 217             final double pa = na / n;
 218             s += pa * Math.log( pa );
 219         }
 220         if ( n < k ) {
 221             return -( s / ( Math.log( n ) ) );
 222         }
 223         else {
 224             return -( s / ( Math.log( k ) ) );
 225         }
 226     }
 227
 228     final public static DescriptiveStatistics calculateEffectiveLengthStatistics( final Msa msa ) {
 229         final DescriptiveStatistics stats = new BasicDescriptiveStatistics();
 230         for( int row = 0; row < msa.getNumberOfSequences(); ++row ) {
 231             final MolecularSequence s = msa.getSequence( row );
 232             stats.addValue( s.getLength() - s.getNumberOfGapResidues() );
 233         }
 234         return stats;
 235     }
 236
 237     final public static DescriptiveStatistics calculateIdentityRatio( final int from, final int to, final Msa msa ) {
 238         final DescriptiveStatistics stats = new BasicDescriptiveStatistics();
 239         for( int c = from; c <= to; ++c ) {
 240             stats.addValue( calculateIdentityRatio( msa, c ) );
 241         }
 242         return stats;
 243     }
 244
 245     final public static double calculateIdentityRatio( final Msa msa, final int column ) {
 246         final SortedMap<Character, Integer> dist = calculateResidueDestributionPerColumn( msa, column );
 247         int majority_count = 0;
 248         final Iterator<Map.Entry<Character, Integer>> it = dist.entrySet().iterator();
 249         while ( it.hasNext() ) {
 250             final Map.Entry<Character, Integer> pair = it.next();
 251             if ( pair.getValue() > majority_count ) {
 252                 majority_count = pair.getValue();
 253             }
 254         }
 255         return ( double ) majority_count / msa.getNumberOfSequences();
 256     }
 257
 258     public static SortedMap<Character, Integer> calculateResidueDestributionPerColumn( final Msa msa, final int column ) {
 259         final SortedMap<Character, Integer> map = new TreeMap<Character, Integer>();
 260         for( final Character r : msa.getColumnAt( column ) ) {
 261             if ( r != MolecularSequence.GAP ) {
 262                 if ( !map.containsKey( r ) ) {
 263                     map.put( r, 1 );
 264                 }
 265                 else {
 266                     map.put( r, map.get( r ) + 1 );
 267                 }
 268             }
 269         }
 270         return map;
 271     }
 272
 273     synchronized public static MsaMethods createInstance() {
 274         return new MsaMethods();
 275     }
 276
 277     final public static Msa removeSequence( final Msa msa, final String to_remove_id ) {
 278         final List<MolecularSequence> seqs = new ArrayList<MolecularSequence>();
 279         for( int row = 0; row < msa.getNumberOfSequences(); ++row ) {
 280             if ( !to_remove_id.equals( msa.getIdentifier( row ) ) ) {
 281                 seqs.add( msa.getSequence( row ) );
 282             }
 283         }
 284         if ( seqs.size() < 1 ) {
 285             return null;
 286         }
 287         return BasicMsa.createInstance( seqs );
 288     }
 289
 290     final public static Msa removeSequences( final Msa msa, final List<String> to_remove_ids ) {
 291         final List<MolecularSequence> seqs = new ArrayList<MolecularSequence>();
 292         for( int row = 0; row < msa.getNumberOfSequences(); ++row ) {
 293             if ( !to_remove_ids.contains( msa.getIdentifier( row ) ) ) {
 294                 seqs.add( msa.getSequence( row ) );
 295             }
 296         }
 297         if ( seqs.size() < 1 ) {
 298             return null;
 299         }
 300         return BasicMsa.createInstance( seqs );
 301     }
 302
 303     public static Msa removeSequencesByMinimalLength( final Msa msa, final int min_effective_length ) {
 304         final List<Integer> to_remove_rows = new ArrayList<Integer>();
 305         for( int seq = 0; seq < msa.getNumberOfSequences(); ++seq ) {
 306             int eff_length = 0;
 307             for( int i = 0; i < msa.getLength(); ++i ) {
 308                 if ( msa.getResidueAt( seq, i ) != MolecularSequence.GAP ) {
 309                     eff_length++;
 310                 }
 311             }
 312             if ( eff_length < min_effective_length ) {
 313                 to_remove_rows.add( seq );
 314             }
 315         }
 316         return removeSequencesByRow( msa, to_remove_rows );
 317     }
 318
 319     final public static Msa removeSequencesByRow( final Msa msa, final List<Integer> to_remove_rows ) {
 320         final List<MolecularSequence> seqs = new ArrayList<MolecularSequence>();
 321         for( int row = 0; row < msa.getNumberOfSequences(); ++row ) {
 322             if ( !to_remove_rows.contains( row ) ) {
 323                 seqs.add( msa.getSequence( row ) );
 324             }
 325         }
 326         if ( seqs.size() < 1 ) {
 327             return null;
 328         }
 329         return BasicMsa.createInstance( seqs );
 330     }
 331
 332     final private static HashMap<Character, Integer> calcResidueDistribution20( final Msa msa, final int col ) {
 333         final HashMap<Character, Integer> counts = new HashMap<Character, Integer>();
 334         for( int row = 0; row < msa.getNumberOfSequences(); ++row ) {
 335             final char c = msa.getResidueAt( row, col );
 336             if ( c != MolecularSequence.GAP ) {
 337                 if ( !counts.containsKey( c ) ) {
 338                     counts.put( c, 1 );
 339                 }
 340                 else {
 341                     counts.put( c, 1 + counts.get( c ) );
 342                 }
 343             }
 344         }
 345         return counts;
 346     }
 347
 348     final private static HashMap<Character, Integer> calcResidueDistribution21( final Msa msa, final int col ) {
 349         final HashMap<Character, Integer> counts = new HashMap<Character, Integer>();
 350         for( int row = 0; row < msa.getNumberOfSequences(); ++row ) {
 351             final char c = msa.getResidueAt( row, col );
 352             if ( !counts.containsKey( c ) ) {
 353                 counts.put( c, 1 );
 354             }
 355             else {
 356                 counts.put( c, 1 + counts.get( c ) );
 357             }
 358         }
 359         return counts;
 360     }
 361
 362     final private static HashMap<Character, Integer> calcResidueDistribution6( final Msa msa, final int col ) {
 363         // Residues are classified into one of tex2html_wrap199 types:
 364         // aliphatic [AVLIMC], aromatic [FWYH], polar [STNQ], positive [KR], negative [DE],
 365         // special conformations [GP] and gaps. This convention follows that
 366         // of Mirny & Shakhnovich (1999, J Mol Biol 291:177-196).
 367         final HashMap<Character, Integer> counts = new HashMap<Character, Integer>();
 368         for( int row = 0; row < msa.getNumberOfSequences(); ++row ) {
 369             final char c = msa.getResidueAt( row, col );
 370             char x;
 371             if ( ( c == 'A' ) || ( c == 'V' ) || ( c == 'L' ) || ( c == 'I' ) || ( c == 'M' ) || ( c == 'C' ) ) {
 372                 // aliphatic
 373                 x = 'a';
 374             }
 375             else if ( ( c == 'F' ) || ( c == 'W' ) || ( c == 'Y' ) || ( c == 'H' ) ) {
 376                 // aromatic
 377                 x = 'r';
 378             }
 379             else if ( ( c == 'S' ) || ( c == 'T' ) || ( c == 'N' ) || ( c == 'Q' ) ) {
 380                 // polar
 381                 x = 'p';
 382             }
 383             else if ( ( c == 'K' ) || ( c == 'R' ) ) {
 384                 // positive
 385                 x = 'o';
 386             }
 387             else if ( ( c == 'D' ) || ( c == 'E' ) ) {
 388                 // negative
 389                 x = 'e';
 390             }
 391             else if ( ( c == 'G' ) || ( c == 'P' ) ) {
 392                 // aliphatic - special conformation
 393                 x = 's';
 394             }
 395             else {
 396                 continue;
 397             }
 398             if ( !counts.containsKey( x ) ) {
 399                 counts.put( x, 1 );
 400             }
 401             else {
 402                 counts.put( x, 1 + counts.get( x ) );
 403             }
 404         }
 405         return counts;
 406     }
 407
 408     final private static HashMap<Character, Integer> calcResidueDistribution7( final Msa msa, final int col ) {
 409         // Residues are classified into one of tex2html_wrap199 types:
 410         // aliphatic [AVLIMC], aromatic [FWYH], polar [STNQ], positive [KR], negative [DE],
 411         // special conformations [GP] and gaps. This convention follows that
 412         // of Mirny & Shakhnovich (1999, J Mol Biol 291:177-196).
 413         final HashMap<Character, Integer> counts = new HashMap<Character, Integer>();
 414         for( int row = 0; row < msa.getNumberOfSequences(); ++row ) {
 415             final char c = msa.getResidueAt( row, col );
 416             char x = '-';
 417             if ( ( c == 'A' ) || ( c == 'V' ) || ( c == 'L' ) || ( c == 'I' ) || ( c == 'M' ) || ( c == 'C' ) ) {
 418                 // aliphatic
 419                 x = 'a';
 420             }
 421             else if ( ( c == 'F' ) || ( c == 'W' ) || ( c == 'Y' ) || ( c == 'H' ) ) {
 422                 // aromatic
 423                 x = 'r';
 424             }
 425             else if ( ( c == 'S' ) || ( c == 'T' ) || ( c == 'N' ) || ( c == 'Q' ) ) {
 426                 // polar
 427                 x = 'p';
 428             }
 429             else if ( ( c == 'K' ) || ( c == 'R' ) ) {
 430                 // positive
 431                 x = 'o';
 432             }
 433             else if ( ( c == 'D' ) || ( c == 'E' ) ) {
 434                 // negative
 435                 x = 'e';
 436             }
 437             else if ( ( c == 'G' ) || ( c == 'P' ) ) {
 438                 // aliphatic - special conformation
 439                 x = 's';
 440             }
 441             if ( !counts.containsKey( x ) ) {
 442                 counts.put( x, 1 );
 443             }
 444             else {
 445                 counts.put( x, 1 + counts.get( x ) );
 446             }
 447         }
 448         return counts;
 449     }
 450 }