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