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