src/jalview/io/StockholmFile.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.8)
   3  * Copyright (C) 2012 J Procter, AM Waterhouse, LM Lui, J Engelhardt, G Barton, M Clamp, S Searle
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
  10  *
  11  * Jalview is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty
  13  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  14  * PURPOSE.  See the GNU General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  17  */
  18 /*
  19  * This extension was written by Benjamin Schuster-Boeckler at sanger.ac.uk
  20  */
  21 package jalview.io;
  22
  23 import java.io.*;
  24 import java.util.*;
  25
  26 import com.stevesoft.pat.*;
  27 import jalview.datamodel.*;
  28 import jalview.analysis.Rna;
  29
  30 // import org.apache.log4j.*;
  31
  32 /**
  33  * This class is supposed to parse a Stockholm format file into Jalview There
  34  * are TODOs in this class: we do not know what the database source and version
  35  * is for the file when parsing the #GS= AC tag which associates accessions with
  36  * sequences. Database references are also not parsed correctly: a separate
  37  * reference string parser must be added to parse the database reference form
  38  * into Jalview's local representation.
  39  *
  40  * @author bsb at sanger.ac.uk
  41  * @version 0.3 + jalview mods
  42  *
  43  */
  44 public class StockholmFile extends AlignFile
  45 {
  46   // static Logger logger = Logger.getLogger("jalview.io.StockholmFile");
  47
  48   public StockholmFile()
  49   {
  50   }
  51
  52   public StockholmFile(String inFile, String type) throws IOException
  53   {
  54     super(inFile, type);
  55   }
  56
  57   public StockholmFile(FileParse source) throws IOException
  58   {
  59     super(source);
  60   }
  61
  62   public void initData()
  63   {
  64     super.initData();
  65   }
  66
  67   /**
  68    * Parse a file in Stockholm format into Jalview's data model. The file has to
  69    * be passed at construction time
  70    *
  71    * @throws IOException
  72    *           If there is an error with the input file
  73    */
  74   public void parse() throws IOException
  75   {
  76     StringBuffer treeString = new StringBuffer();
  77     String treeName = null;
  78     // --------------- Variable Definitions -------------------
  79     String line;
  80     String version;
  81     // String id;
  82     Hashtable seqAnn = new Hashtable(); // Sequence related annotations
  83     Hashtable seqs = new Hashtable();
  84     Regex p, r, rend, s, x;
  85
  86     // Temporary line for processing RNA annotation
  87     // String RNAannot = "";
  88
  89     // ------------------ Parsing File ----------------------
  90     // First, we have to check that this file has STOCKHOLM format, i.e. the
  91     // first line must match
  92     r = new Regex("# STOCKHOLM ([\\d\\.]+)");
  93     if (!r.search(nextLine()))
  94     {
  95       throw new IOException(
  96               "This file is not in valid STOCKHOLM format: First line does not contain '# STOCKHOLM'");
  97     }
  98     else
  99     {
 100       version = r.stringMatched(1);
 101       // logger.debug("Stockholm version: " + version);
 102     }
 103
 104     // We define some Regexes here that will be used regularily later
 105     rend = new Regex("^\\s*\\/\\/"); // Find the end of an alignment
 106     p = new Regex("(\\S+)\\/(\\d+)\\-(\\d+)"); // split sequence id in
 107     // id/from/to
 108     s = new Regex("(\\S+)\\s+(\\S*)\\s+(.*)"); // Parses annotation subtype
 109     r = new Regex("#=(G[FSRC]?)\\s+(.*)"); // Finds any annotation line
 110     x = new Regex("(\\S+)\\s+(\\S+)"); // split id from sequence
 111
 112     // Convert all bracket types to parentheses (necessary for passing to VARNA)
 113     Regex openparen = new Regex("(<|\\[)", "(");
 114     Regex closeparen = new Regex("(>|\\])", ")");
 115
 116     // Detect if file is RNA by looking for bracket types
 117     Regex detectbrackets = new Regex("(<|>|\\[|\\]|\\(|\\))");
 118
 119     rend.optimize();
 120     p.optimize();
 121     s.optimize();
 122     r.optimize();
 123     x.optimize();
 124     openparen.optimize();
 125     closeparen.optimize();
 126
 127     while ((line = nextLine()) != null)
 128     {
 129       if (line.length() == 0)
 130       {
 131         continue;
 132       }
 133       if (rend.search(line))
 134       {
 135         // End of the alignment, pass stuff back
 136
 137         this.noSeqs = seqs.size();
 138         // logger.debug("Number of sequences: " + this.noSeqs);
 139         Enumeration accs = seqs.keys();
 140         while (accs.hasMoreElements())
 141         {
 142           String acc = (String) accs.nextElement();
 143           // logger.debug("Processing sequence " + acc);
 144           String seq = (String) seqs.remove(acc);
 145           if (maxLength < seq.length())
 146           {
 147             maxLength = seq.length();
 148           }
 149           int start = 1;
 150           int end = -1;
 151           String sid = acc;
 152           /*
 153            * Retrieve hash of annotations for this accession Associate
 154            * Annotation with accession
 155            */
 156           Hashtable accAnnotations = null;
 157
 158           if (seqAnn != null && seqAnn.containsKey(acc))
 159           {
 160             accAnnotations = (Hashtable) seqAnn.remove(acc);
 161             // TODO: add structures to sequence
 162           }
 163
 164           // Split accession in id and from/to
 165           if (p.search(acc))
 166           {
 167             sid = p.stringMatched(1);
 168             start = Integer.parseInt(p.stringMatched(2));
 169             end = Integer.parseInt(p.stringMatched(3));
 170           }
 171           // logger.debug(sid + ", " + start + ", " + end);
 172
 173           Sequence seqO = new Sequence(sid, seq, start, end);
 174           // Add Description (if any)
 175           if (accAnnotations != null && accAnnotations.containsKey("DE"))
 176           {
 177             String desc = (String) accAnnotations.get("DE");
 178             seqO.setDescription((desc == null) ? "" : desc);
 179           }
 180           // Add DB References (if any)
 181           if (accAnnotations != null && accAnnotations.containsKey("DR"))
 182           {
 183             String dbr = (String) accAnnotations.get("DR");
 184             if (dbr != null && dbr.indexOf(";") > -1)
 185             {
 186               String src = dbr.substring(0, dbr.indexOf(";"));
 187               String acn = dbr.substring(dbr.indexOf(";") + 1);
 188               jalview.util.DBRefUtils.parseToDbRef(seqO, src, "0", acn);
 189               // seqO.addDBRef(dbref);
 190             }
 191           }
 192           if (accAnnotations != null && accAnnotations.containsKey("SS"))
 193           {
 194             Vector v = (Vector) accAnnotations.get("SS");
 195
 196             for (int i = 0; i < v.size(); i++)
 197             {
 198               AlignmentAnnotation an = (AlignmentAnnotation) v.elementAt(i);
 199               seqO.addAlignmentAnnotation(an);
 200               // annotations.add(an);
 201             }
 202           }
 203
 204           Hashtable features = null;
 205           // We need to adjust the positions of all features to account for gaps
 206           try
 207           {
 208             features = (Hashtable) accAnnotations.remove("features");
 209           } catch (java.lang.NullPointerException e)
 210           {
 211             // loggerwarn("Getting Features for " + acc + ": " +
 212             // e.getMessage());
 213             // continue;
 214           }
 215           // if we have features
 216           if (features != null)
 217           {
 218             int posmap[] = seqO.findPositionMap();
 219             Enumeration i = features.keys();
 220             while (i.hasMoreElements())
 221             {
 222               // TODO: parse out secondary structure annotation as annotation
 223               // row
 224               // TODO: parse out scores as annotation row
 225               // TODO: map coding region to core jalview feature types
 226               String type = i.nextElement().toString();
 227               Hashtable content = (Hashtable) features.remove(type);
 228               Enumeration j = content.keys();
 229               while (j.hasMoreElements())
 230               {
 231                 String desc = j.nextElement().toString();
 232                 String ns = content.get(desc).toString();
 233                 char[] byChar = ns.toCharArray();
 234                 for (int k = 0; k < byChar.length; k++)
 235                 {
 236                   char c = byChar[k];
 237                   if (!(c == ' ' || c == '_' || c == '-' || c == '.')) // PFAM
 238                   // uses
 239                   // '.'
 240                   // for
 241                   // feature
 242                   // background
 243                   {
 244                     int new_pos = posmap[k]; // look up nearest seqeunce
 245                     // position to this column
 246                     SequenceFeature feat = new SequenceFeature(type, desc,
 247                             new_pos, new_pos, 0f, null);
 248
 249                     seqO.addSequenceFeature(feat);
 250                   }
 251                 }
 252               }
 253
 254             }
 255
 256           }
 257           // garbage collect
 258
 259           // logger.debug("Adding seq " + acc + " from " + start + " to " + end
 260           // + ": " + seq);
 261           this.seqs.addElement(seqO);
 262         }
 263         return; // finished parsing this segment of source
 264       }
 265       else if (!r.search(line))
 266       {
 267         // System.err.println("Found sequence line: " + line);
 268
 269         // Split sequence in sequence and accession parts
 270         if (!x.search(line))
 271         {
 272           // logger.error("Could not parse sequence line: " + line);
 273           throw new IOException("Could not parse sequence line: " + line);
 274         }
 275         String ns = (String) seqs.get(x.stringMatched(1));
 276         if (ns == null)
 277         {
 278           ns = "";
 279         }
 280         ns += x.stringMatched(2);
 281
 282         seqs.put(x.stringMatched(1), ns);
 283       }
 284       else
 285       {
 286         String annType = r.stringMatched(1);
 287         String annContent = r.stringMatched(2);
 288
 289         // System.err.println("type:" + annType + " content: " + annContent);
 290
 291         if (annType.equals("GF"))
 292         {
 293           /*
 294            * Generic per-File annotation, free text Magic features: #=GF NH
 295            * <tree in New Hampshire eXtended format> #=GF TN <Unique identifier
 296            * for the next tree> Pfam descriptions: 7. DESCRIPTION OF FIELDS
 297            *
 298            * Compulsory fields: ------------------
 299            *
 300            * AC Accession number: Accession number in form PFxxxxx.version or
 301            * PBxxxxxx. ID Identification: One word name for family. DE
 302            * Definition: Short description of family. AU Author: Authors of the
 303            * entry. SE Source of seed: The source suggesting the seed members
 304            * belong to one family. GA Gathering method: Search threshold to
 305            * build the full alignment. TC Trusted Cutoff: Lowest sequence score
 306            * and domain score of match in the full alignment. NC Noise Cutoff:
 307            * Highest sequence score and domain score of match not in full
 308            * alignment. TP Type: Type of family -- presently Family, Domain,
 309            * Motif or Repeat. SQ Sequence: Number of sequences in alignment. AM
 310            * Alignment Method The order ls and fs hits are aligned to the model
 311            * to build the full align. // End of alignment.
 312            *
 313            * Optional fields: ----------------
 314            *
 315            * DC Database Comment: Comment about database reference. DR Database
 316            * Reference: Reference to external database. RC Reference Comment:
 317            * Comment about literature reference. RN Reference Number: Reference
 318            * Number. RM Reference Medline: Eight digit medline UI number. RT
 319            * Reference Title: Reference Title. RA Reference Author: Reference
 320            * Author RL Reference Location: Journal location. PI Previous
 321            * identifier: Record of all previous ID lines. KW Keywords: Keywords.
 322            * CC Comment: Comments. NE Pfam accession: Indicates a nested domain.
 323            * NL Location: Location of nested domains - sequence ID, start and
 324            * end of insert.
 325            *
 326            * Obsolete fields: ----------- AL Alignment method of seed: The
 327            * method used to align the seed members.
 328            */
 329           // Let's save the annotations, maybe we'll be able to do something
 330           // with them later...
 331           Regex an = new Regex("(\\w+)\\s*(.*)");
 332           if (an.search(annContent))
 333           {
 334             if (an.stringMatched(1).equals("NH"))
 335             {
 336               treeString.append(an.stringMatched(2));
 337             }
 338             else if (an.stringMatched(1).equals("TN"))
 339             {
 340               if (treeString.length() > 0)
 341               {
 342                 if (treeName == null)
 343                 {
 344                   treeName = "Tree " + (getTreeCount() + 1);
 345                 }
 346                 addNewickTree(treeName, treeString.toString());
 347               }
 348               treeName = an.stringMatched(2);
 349               treeString = new StringBuffer();
 350             }
 351             setAlignmentProperty(an.stringMatched(1), an.stringMatched(2));
 352           }
 353         }
 354         else if (annType.equals("GS"))
 355         {
 356           // Generic per-Sequence annotation, free text
 357           /*
 358            * Pfam uses these features: Feature Description ---------------------
 359            * ----------- AC <accession> ACcession number DE <freetext>
 360            * DEscription DR <db>; <accession>; Database Reference OS <organism>
 361            * OrganiSm (species) OC <clade> Organism Classification (clade, etc.)
 362            * LO <look> Look (Color, etc.)
 363            */
 364           if (s.search(annContent))
 365           {
 366             String acc = s.stringMatched(1);
 367             String type = s.stringMatched(2);
 368             String content = s.stringMatched(3);
 369             // TODO: store DR in a vector.
 370             // TODO: store AC according to generic file db annotation.
 371             Hashtable ann;
 372             if (seqAnn.containsKey(acc))
 373             {
 374               ann = (Hashtable) seqAnn.get(acc);
 375             }
 376             else
 377             {
 378               ann = new Hashtable();
 379             }
 380             ann.put(type, content);
 381             seqAnn.put(acc, ann);
 382           }
 383           else
 384           {
 385             throw new IOException("Error parsing " + line);
 386           }
 387         }
 388         else if (annType.equals("GC"))
 389         {
 390           // Generic per-Column annotation, exactly 1 char per column
 391           // always need a label.
 392           if (x.search(annContent))
 393           {
 394             // parse out and create alignment annotation directly.
 395             parseAnnotationRow(annotations, x.stringMatched(1),
 396                     x.stringMatched(2));
 397           }
 398         }
 399         else if (annType.equals("GR"))
 400         {
 401           // Generic per-Sequence AND per-Column markup, exactly 1 char per
 402           // column
 403           /*
 404            * Feature Description Markup letters ------- -----------
 405            * -------------- SS Secondary Structure [HGIEBTSCX] SA Surface
 406            * Accessibility [0-9X] (0=0%-10%; ...; 9=90%-100%) TM TransMembrane
 407            * [Mio] PP Posterior Probability [0-9*] (0=0.00-0.05; 1=0.05-0.15;
 408            * *=0.95-1.00) LI LIgand binding [*] AS Active Site [*] IN INtron (in
 409            * or after) [0-2]
 410            */
 411           if (s.search(annContent))
 412           {
 413             String acc = s.stringMatched(1);
 414             String type = s.stringMatched(2);
 415             String seq = new String(s.stringMatched(3));
 416             String description = null;
 417             // Check for additional information about the current annotation
 418             // We use a simple string tokenizer here for speed
 419             StringTokenizer sep = new StringTokenizer(seq, " \t");
 420             description = sep.nextToken();
 421             if (sep.hasMoreTokens())
 422             {
 423               seq = sep.nextToken();
 424             }
 425             else
 426             {
 427               seq = description;
 428               description = new String();
 429             }
 430             // sequence id with from-to fields
 431
 432             Hashtable ann;
 433             // Get an object with all the annotations for this sequence
 434             if (seqAnn.containsKey(acc))
 435             {
 436               // logger.debug("Found annotations for " + acc);
 437               ann = (Hashtable) seqAnn.get(acc);
 438             }
 439             else
 440             {
 441               // logger.debug("Creating new annotations holder for " + acc);
 442               ann = new Hashtable();
 443               seqAnn.put(acc, ann);
 444             }
 445             // TODO test structure, call parseAnnotationRow with vector from
 446             // hashtable for specific sequence
 447             Hashtable features;
 448             // Get an object with all the content for an annotation
 449             if (ann.containsKey("features"))
 450             {
 451               // logger.debug("Found features for " + acc);
 452               features = (Hashtable) ann.get("features");
 453             }
 454             else
 455             {
 456               // logger.debug("Creating new features holder for " + acc);
 457               features = new Hashtable();
 458               ann.put("features", features);
 459             }
 460
 461             Hashtable content;
 462             if (features.containsKey(this.id2type(type)))
 463             {
 464               // logger.debug("Found content for " + this.id2type(type));
 465               content = (Hashtable) features.get(this.id2type(type));
 466             }
 467             else
 468             {
 469               // logger.debug("Creating new content holder for " +
 470               // this.id2type(type));
 471               content = new Hashtable();
 472               features.put(this.id2type(type), content);
 473             }
 474             String ns = (String) content.get(description);
 475             if (ns == null)
 476             {
 477               ns = "";
 478             }
 479             ns += seq;
 480             content.put(description, ns);
 481
 482             if (type.equals("SS"))
 483             {
 484               Hashtable strucAnn;
 485               if (seqAnn.containsKey(acc))
 486               {
 487                 strucAnn = (Hashtable) seqAnn.get(acc);
 488               }
 489               else
 490               {
 491                 strucAnn = new Hashtable();
 492               }
 493
 494               Vector newStruc = new Vector();
 495               parseAnnotationRow(newStruc, type, ns);
 496
 497               strucAnn.put(type, newStruc);
 498               seqAnn.put(acc, strucAnn);
 499             }
 500           }
 501           else
 502           {
 503             System.err
 504                     .println("Warning - couldn't parse sequence annotation row line:\n"
 505                             + line);
 506             // throw new IOException("Error parsing " + line);
 507           }
 508         }
 509         else
 510         {
 511           throw new IOException("Unknown annotation detected: " + annType
 512                   + " " + annContent);
 513         }
 514       }
 515     }
 516     if (treeString.length() > 0)
 517     {
 518       if (treeName == null)
 519       {
 520         treeName = "Tree " + (1 + getTreeCount());
 521       }
 522       addNewickTree(treeName, treeString.toString());
 523     }
 524   }
 525
 526   protected static AlignmentAnnotation parseAnnotationRow(
 527           Vector annotation, String label, String annots)
 528   {
 529     String convert1, convert2 = null;
 530
 531     // Convert all bracket types to parentheses
 532     Regex openparen = new Regex("(<|\\[)", "(");
 533     Regex closeparen = new Regex("(>|\\])", ")");
 534
 535     // Detect if file is RNA by looking for bracket types
 536     Regex detectbrackets = new Regex("(<|>|\\[|\\]|\\(|\\))");
 537
 538     convert1 = openparen.replaceAll(annots);
 539     convert2 = closeparen.replaceAll(convert1);
 540     annots = convert2;
 541
 542     String type = (label.indexOf("_cons") == label.length() - 5) ? label
 543             .substring(0, label.length() - 5) : label;
 544     boolean ss = false;
 545     type = id2type(type);
 546     if (type.equals("secondary structure"))
 547     {
 548       ss = true;
 549     }
 550     // decide on secondary structure or not.
 551     Annotation[] els = new Annotation[annots.length()];
 552     for (int i = 0; i < annots.length(); i++)
 553     {
 554       String pos = annots.substring(i, i + 1);
 555       Annotation ann;
 556       ann = new Annotation(pos, "", ' ', 0f); // 0f is 'valid' null - will not
 557       // be written out
 558       if (ss)
 559       {
 560         if (detectbrackets.search(pos))
 561         {
 562           ann.secondaryStructure = jalview.schemes.ResidueProperties
 563                   .getRNASecStrucState(pos).charAt(0);
 564         }
 565         else
 566         {
 567           ann.secondaryStructure = jalview.schemes.ResidueProperties
 568                   .getDssp3state(pos).charAt(0);
 569         }
 570
 571         if (ann.secondaryStructure == pos.charAt(0) || pos.charAt(0) == 'C')
 572         {
 573           ann.displayCharacter = ""; // null; // " ";
 574         }
 575         else
 576         {
 577           ann.displayCharacter = " " + ann.displayCharacter;
 578         }
 579       }
 580
 581       els[i] = ann;
 582     }
 583     AlignmentAnnotation annot = null;
 584     Enumeration e = annotation.elements();
 585     while (e.hasMoreElements())
 586     {
 587       annot = (AlignmentAnnotation) e.nextElement();
 588       if (annot.label.equals(type))
 589         break;
 590       annot = null;
 591     }
 592     if (annot == null)
 593     {
 594       annot = new AlignmentAnnotation(type, type, els);
 595       annotation.addElement(annot);
 596     }
 597     else
 598     {
 599       Annotation[] anns = new Annotation[annot.annotations.length
 600               + els.length];
 601       System.arraycopy(annot.annotations, 0, anns, 0,
 602               annot.annotations.length);
 603       System.arraycopy(els, 0, anns, annot.annotations.length, els.length);
 604       annot.annotations = anns;
 605       // System.out.println("else: ");
 606     }
 607     return annot;
 608   }
 609
 610   public static String print(SequenceI[] s)
 611   {
 612     return "not yet implemented";
 613   }
 614
 615   public String print()
 616   {
 617     return print(getSeqsAsArray());
 618   }
 619
 620   private static Hashtable typeIds = null;
 621   static
 622   {
 623     if (typeIds == null)
 624     {
 625       typeIds = new Hashtable();
 626       typeIds.put("SS", "secondary structure");
 627       typeIds.put("SA", "surface accessibility");
 628       typeIds.put("TM", "transmembrane");
 629       typeIds.put("PP", "posterior probability");
 630       typeIds.put("LI", "ligand binding");
 631       typeIds.put("AS", "active site");
 632       typeIds.put("IN", "intron");
 633       typeIds.put("IR", "interacting residue");
 634       typeIds.put("AC", "accession");
 635       typeIds.put("OS", "organism");
 636       typeIds.put("CL", "class");
 637       typeIds.put("DE", "description");
 638       typeIds.put("DR", "reference");
 639       typeIds.put("LO", "look");
 640       typeIds.put("RF", "reference positions");
 641
 642     }
 643   }
 644
 645   protected static String id2type(String id)
 646   {
 647     if (typeIds.containsKey(id))
 648     {
 649       return (String) typeIds.get(id);
 650     }
 651     System.err.println("Warning : Unknown Stockholm annotation type code "
 652             + id);
 653     return id;
 654   }
 655   /**
 656    * //ssline is complete secondary structure line private AlignmentAnnotation
 657    * addHelices(Vector annotation, String label, String ssline) {
 658    *
 659    * // decide on secondary structure or not. Annotation[] els = new
 660    * Annotation[ssline.length()]; for (int i = 0; i < ssline.length(); i++) {
 661    * String pos = ssline.substring(i, i + 1); Annotation ann; ann = new
 662    * Annotation(pos, "", ' ', 0f); // 0f is 'valid' null - will not
 663    *
 664    * ann.secondaryStructure =
 665    * jalview.schemes.ResidueProperties.getRNAssState(pos).charAt(0);
 666    *
 667    * ann.displayCharacter = "x" + ann.displayCharacter;
 668    *
 669    * System.out.println(ann.displayCharacter);
 670    *
 671    * els[i] = ann; } AlignmentAnnotation helicesAnnot = null; Enumeration e =
 672    * annotation.elements(); while (e.hasMoreElements()) { helicesAnnot =
 673    * (AlignmentAnnotation) e.nextElement(); if (helicesAnnot.label.equals(type))
 674    * break; helicesAnnot = null; } if (helicesAnnot == null) { helicesAnnot =
 675    * new AlignmentAnnotation(type, type, els);
 676    * annotation.addElement(helicesAnnot); } else { Annotation[] anns = new
 677    * Annotation[helicesAnnot.annotations.length + els.length];
 678    * System.arraycopy(helicesAnnot.annotations, 0, anns, 0,
 679    * helicesAnnot.annotations.length); System.arraycopy(els, 0, anns,
 680    * helicesAnnot.annotations.length, els.length); helicesAnnot.annotations =
 681    * anns; }
 682    *
 683    * helicesAnnot.features = Rna.GetBasePairs(ssline);
 684    * Rna.HelixMap(helicesAnnot.features);
 685    *
 686    *
 687    * return helicesAnnot; }
 688    */
 689 }