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