src/jalview/ws/sifts/SiftsClient.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   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
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 package jalview.ws.sifts;
  22
  23 import jalview.analysis.AlignSeq;
  24 import jalview.analysis.scoremodels.ScoreMatrix;
  25 import jalview.analysis.scoremodels.ScoreModels;
  26 import jalview.api.DBRefEntryI;
  27 import jalview.api.SiftsClientI;
  28 import jalview.datamodel.DBRefEntry;
  29 import jalview.datamodel.DBRefSource;
  30 import jalview.datamodel.SequenceI;
  31 import jalview.io.StructureFile;
  32 import jalview.schemes.ResidueProperties;
  33 import jalview.structure.StructureMapping;
  34 import jalview.util.Comparison;
  35 import jalview.util.DBRefUtils;
  36 import jalview.util.Format;
  37 import jalview.xml.binding.sifts.Entry;
  38 import jalview.xml.binding.sifts.Entry.Entity;
  39 import jalview.xml.binding.sifts.Entry.Entity.Segment;
  40 import jalview.xml.binding.sifts.Entry.Entity.Segment.ListMapRegion.MapRegion;
  41 import jalview.xml.binding.sifts.Entry.Entity.Segment.ListResidue.Residue;
  42 import jalview.xml.binding.sifts.Entry.Entity.Segment.ListResidue.Residue.CrossRefDb;
  43 import jalview.xml.binding.sifts.Entry.Entity.Segment.ListResidue.Residue.ResidueDetail;
  44
  45 import java.io.File;
  46 import java.io.FileInputStream;
  47 import java.io.FileOutputStream;
  48 import java.io.IOException;
  49 import java.io.InputStream;
  50 import java.io.PrintStream;
  51 import java.net.URL;
  52 import java.net.URLConnection;
  53 import java.nio.file.Files;
  54 import java.nio.file.Path;
  55 import java.nio.file.attribute.BasicFileAttributes;
  56 import java.util.ArrayList;
  57 import java.util.Arrays;
  58 import java.util.Collection;
  59 import java.util.Collections;
  60 import java.util.Date;
  61 import java.util.HashMap;
  62 import java.util.HashSet;
  63 import java.util.List;
  64 import java.util.Map;
  65 import java.util.Set;
  66 import java.util.TreeMap;
  67 import java.util.zip.GZIPInputStream;
  68
  69 import javax.xml.bind.JAXBContext;
  70 import javax.xml.bind.Unmarshaller;
  71 import javax.xml.stream.XMLInputFactory;
  72 import javax.xml.stream.XMLStreamReader;
  73
  74 import MCview.Atom;
  75 import MCview.PDBChain;
  76
  77 public class SiftsClient implements SiftsClientI
  78 {
  79   /*
  80    * for use in mocking out file fetch for tests only
  81    * - reset to null after testing!
  82    */
  83   private static File mockSiftsFile;
  84
  85   private Entry siftsEntry;
  86
  87   private StructureFile pdb;
  88
  89   private String pdbId;
  90
  91   private String structId;
  92
  93   private CoordinateSys seqCoordSys = CoordinateSys.UNIPROT;
  94
  95   /**
  96    * PDB sequence position to sequence coordinate mapping as derived from SIFTS
  97    * record for the identified SeqCoordSys Used for lift-over from sequence
  98    * derived from PDB (with first extracted PDBRESNUM as 'start' to the sequence
  99    * being annotated with PDB data
 100    */
 101   private jalview.datamodel.Mapping seqFromPdbMapping;
 102
 103   private static final int BUFFER_SIZE = 4096;
 104
 105   public static final int UNASSIGNED = Integer.MIN_VALUE;
 106
 107   private static final int PDB_RES_POS = 0;
 108
 109   private static final int PDB_ATOM_POS = 1;
 110
 111   private static final int PDBE_POS = 2;
 112
 113   private static final String NOT_OBSERVED = "Not_Observed";
 114
 115   private static final String SIFTS_FTP_BASE_URL = "http://ftp.ebi.ac.uk/pub/databases/msd/sifts/xml/";
 116
 117   private final static String NEWLINE = System.lineSeparator();
 118
 119   private String curSourceDBRef;
 120
 121   private HashSet<String> curDBRefAccessionIdsString;
 122
 123   private enum CoordinateSys
 124   {
 125     UNIPROT("UniProt"), PDB("PDBresnum"), PDBe("PDBe");
 126     private String name;
 127
 128     private CoordinateSys(String name)
 129     {
 130       this.name = name;
 131     }
 132
 133     public String getName()
 134     {
 135       return name;
 136     }
 137   };
 138
 139   private enum ResidueDetailType
 140   {
 141     NAME_SEC_STRUCTURE("nameSecondaryStructure"),
 142     CODE_SEC_STRUCTURE("codeSecondaryStructure"), ANNOTATION("Annotation");
 143     private String code;
 144
 145     private ResidueDetailType(String code)
 146     {
 147       this.code = code;
 148     }
 149
 150     public String getCode()
 151     {
 152       return code;
 153     }
 154   };
 155
 156   /**
 157    * Fetch SIFTs file for the given PDBfile and construct an instance of
 158    * SiftsClient
 159    *
 160    * @param pdbId
 161    * @throws SiftsException
 162    */
 163   public SiftsClient(StructureFile pdb) throws SiftsException
 164   {
 165     this.pdb = pdb;
 166     this.pdbId = pdb.getId();
 167     File siftsFile = getSiftsFile(pdbId);
 168     siftsEntry = parseSIFTs(siftsFile);
 169   }
 170
 171   /**
 172    * Parse the given SIFTs File and return a JAXB POJO of parsed data
 173    *
 174    * @param siftFile
 175    *          - the GZipped SIFTs XML file to parse
 176    * @return
 177    * @throws Exception
 178    *           if a problem occurs while parsing the SIFTs XML
 179    */
 180   private Entry parseSIFTs(File siftFile) throws SiftsException
 181   {
 182     try (InputStream in = new FileInputStream(siftFile);
 183             GZIPInputStream gzis = new GZIPInputStream(in);)
 184     {
 185       // System.out.println("File : " + siftFile.getAbsolutePath());
 186       JAXBContext jc = JAXBContext.newInstance("jalview.xml.binding.sifts");
 187       XMLStreamReader streamReader = XMLInputFactory.newInstance()
 188               .createXMLStreamReader(gzis);
 189       Unmarshaller um = jc.createUnmarshaller();
 190       return (Entry) um.unmarshal(streamReader);
 191     } catch (Exception e)
 192     {
 193       e.printStackTrace();
 194       throw new SiftsException(e.getMessage());
 195     }
 196   }
 197
 198   /**
 199    * Get a SIFTs XML file for a given PDB Id from Cache or download from FTP
 200    * repository if not found in cache
 201    *
 202    * @param pdbId
 203    * @return SIFTs XML file
 204    * @throws SiftsException
 205    */
 206   public static File getSiftsFile(String pdbId) throws SiftsException
 207   {
 208     /*
 209      * return mocked file if it has been set
 210      */
 211     if (mockSiftsFile != null)
 212     {
 213       return mockSiftsFile;
 214     }
 215
 216     String siftsFileName = SiftsSettings.getSiftDownloadDirectory()
 217             + pdbId.toLowerCase() + ".xml.gz";
 218     File siftsFile = new File(siftsFileName);
 219     if (siftsFile.exists())
 220     {
 221       // The line below is required for unit testing... don't comment it out!!!
 222       System.out.println(">>> SIFTS File already downloaded for " + pdbId);
 223
 224       if (isFileOlderThanThreshold(siftsFile,
 225               SiftsSettings.getCacheThresholdInDays()))
 226       {
 227         File oldSiftsFile = new File(siftsFileName + "_old");
 228         siftsFile.renameTo(oldSiftsFile);
 229         try
 230         {
 231           siftsFile = downloadSiftsFile(pdbId.toLowerCase());
 232           oldSiftsFile.delete();
 233           return siftsFile;
 234         } catch (IOException e)
 235         {
 236           e.printStackTrace();
 237           oldSiftsFile.renameTo(siftsFile);
 238           return new File(siftsFileName);
 239         }
 240       }
 241       else
 242       {
 243         return siftsFile;
 244       }
 245     }
 246     try
 247     {
 248       siftsFile = downloadSiftsFile(pdbId.toLowerCase());
 249     } catch (IOException e)
 250     {
 251       throw new SiftsException(e.getMessage());
 252     }
 253     return siftsFile;
 254   }
 255
 256   /**
 257    * This method enables checking if a cached file has exceeded a certain
 258    * threshold(in days)
 259    *
 260    * @param file
 261    *          the cached file
 262    * @param noOfDays
 263    *          the threshold in days
 264    * @return
 265    */
 266   public static boolean isFileOlderThanThreshold(File file, int noOfDays)
 267   {
 268     Path filePath = file.toPath();
 269     BasicFileAttributes attr;
 270     int diffInDays = 0;
 271     try
 272     {
 273       attr = Files.readAttributes(filePath, BasicFileAttributes.class);
 274       diffInDays = (int) ((new Date().getTime()
 275               - attr.lastModifiedTime().toMillis())
 276               / (1000 * 60 * 60 * 24));
 277       // System.out.println("Diff in days : " + diffInDays);
 278     } catch (IOException e)
 279     {
 280       e.printStackTrace();
 281     }
 282     return noOfDays <= diffInDays;
 283   }
 284
 285   /**
 286    * Download a SIFTs XML file for a given PDB Id from an FTP repository
 287    *
 288    * @param pdbId
 289    * @return downloaded SIFTs XML file
 290    * @throws SiftsException
 291    * @throws IOException
 292    */
 293   public static File downloadSiftsFile(String pdbId)
 294           throws SiftsException, IOException
 295   {
 296     if (pdbId.contains(".cif"))
 297     {
 298       pdbId = pdbId.replace(".cif", "");
 299     }
 300     String siftFile = pdbId + ".xml.gz";
 301     String siftsFileFTPURL = SIFTS_FTP_BASE_URL + siftFile;
 302     String downloadedSiftsFile = SiftsSettings.getSiftDownloadDirectory()
 303             + siftFile;
 304     File siftsDownloadDir = new File(
 305             SiftsSettings.getSiftDownloadDirectory());
 306     if (!siftsDownloadDir.exists())
 307     {
 308       siftsDownloadDir.mkdirs();
 309     }
 310     // System.out.println(">> Download ftp url : " + siftsFileFTPURL);
 311     // long now = System.currentTimeMillis();
 312     URL url = new URL(siftsFileFTPURL);
 313     URLConnection conn = url.openConnection();
 314     InputStream inputStream = conn.getInputStream();
 315     FileOutputStream outputStream = new FileOutputStream(
 316             downloadedSiftsFile);
 317     byte[] buffer = new byte[BUFFER_SIZE];
 318     int bytesRead = -1;
 319     while ((bytesRead = inputStream.read(buffer)) != -1)
 320     {
 321       outputStream.write(buffer, 0, bytesRead);
 322     }
 323     outputStream.close();
 324     inputStream.close();
 325     // System.out.println(">>> File downloaded : " + downloadedSiftsFile
 326     // + " took " + (System.currentTimeMillis() - now) + "ms");
 327     return new File(downloadedSiftsFile);
 328   }
 329
 330   /**
 331    * Delete the SIFTs file for the given PDB Id in the local SIFTs download
 332    * directory
 333    *
 334    * @param pdbId
 335    * @return true if the file was deleted or doesn't exist
 336    */
 337   public static boolean deleteSiftsFileByPDBId(String pdbId)
 338   {
 339     File siftsFile = new File(SiftsSettings.getSiftDownloadDirectory()
 340             + pdbId.toLowerCase() + ".xml.gz");
 341     if (siftsFile.exists())
 342     {
 343       return siftsFile.delete();
 344     }
 345     return true;
 346   }
 347
 348   /**
 349    * Get a valid SIFTs DBRef for the given sequence current SIFTs entry
 350    *
 351    * @param seq
 352    *          - the target sequence for the operation
 353    * @return a valid DBRefEntry that is SIFTs compatible
 354    * @throws Exception
 355    *           if no valid source DBRefEntry was found for the given sequences
 356    */
 357   public DBRefEntryI getValidSourceDBRef(SequenceI seq)
 358           throws SiftsException
 359   {
 360     List<DBRefEntry> dbRefs = seq.getPrimaryDBRefs();
 361     if (dbRefs == null || dbRefs.size() < 1)
 362     {
 363       throw new SiftsException(
 364               "Source DBRef could not be determined. DBRefs might not have been retrieved.");
 365     }
 366
 367     for (DBRefEntry dbRef : dbRefs)
 368     {
 369       if (dbRef == null || dbRef.getAccessionId() == null
 370               || dbRef.getSource() == null)
 371       {
 372         continue;
 373       }
 374       String canonicalSource = DBRefUtils
 375               .getCanonicalName(dbRef.getSource());
 376       if (isValidDBRefEntry(dbRef)
 377               && (canonicalSource.equalsIgnoreCase(DBRefSource.UNIPROT)
 378                       || canonicalSource.equalsIgnoreCase(DBRefSource.PDB)))
 379       {
 380         return dbRef;
 381       }
 382     }
 383     throw new SiftsException("Could not get source DB Ref");
 384   }
 385
 386   /**
 387    * Check that the DBRef Entry is properly populated and is available in this
 388    * SiftClient instance
 389    *
 390    * @param entry
 391    *          - DBRefEntry to validate
 392    * @return true validation is successful otherwise false is returned.
 393    */
 394   boolean isValidDBRefEntry(DBRefEntryI entry)
 395   {
 396     return entry != null && entry.getAccessionId() != null
 397             && isFoundInSiftsEntry(entry.getAccessionId());
 398   }
 399
 400   @Override
 401   public HashSet<String> getAllMappingAccession()
 402   {
 403     HashSet<String> accessions = new HashSet<String>();
 404     List<Entity> entities = siftsEntry.getEntity();
 405     for (Entity entity : entities)
 406     {
 407       List<Segment> segments = entity.getSegment();
 408       for (Segment segment : segments)
 409       {
 410         List<MapRegion> mapRegions = segment.getListMapRegion()
 411                 .getMapRegion();
 412         for (MapRegion mapRegion : mapRegions)
 413         {
 414           accessions
 415                   .add(mapRegion.getDb().getDbAccessionId().toLowerCase());
 416         }
 417       }
 418     }
 419     return accessions;
 420   }
 421
 422   @Override
 423   public StructureMapping getSiftsStructureMapping(SequenceI seq,
 424           String pdbFile, String chain) throws SiftsException
 425   {
 426     SequenceI aseq = seq;
 427     while (seq.getDatasetSequence() != null)
 428     {
 429       seq = seq.getDatasetSequence();
 430     }
 431     structId = (chain == null) ? pdbId : pdbId + "|" + chain;
 432     System.out.println("Getting SIFTS mapping for " + structId + ": seq "
 433             + seq.getName());
 434
 435     final StringBuilder mappingDetails = new StringBuilder(128);
 436     PrintStream ps = new PrintStream(System.out)
 437     {
 438       @Override
 439       public void print(String x)
 440       {
 441         mappingDetails.append(x);
 442       }
 443
 444       @Override
 445       public void println()
 446       {
 447         mappingDetails.append(NEWLINE);
 448       }
 449     };
 450     HashMap<Integer, int[]> mapping = getGreedyMapping(chain, seq, ps);
 451
 452     String mappingOutput = mappingDetails.toString();
 453     StructureMapping siftsMapping = new StructureMapping(aseq, pdbFile,
 454             pdbId, chain, mapping, mappingOutput, seqFromPdbMapping);
 455
 456     return siftsMapping;
 457   }
 458
 459   @Override
 460   public HashMap<Integer, int[]> getGreedyMapping(String entityId,
 461           SequenceI seq, java.io.PrintStream os) throws SiftsException
 462   {
 463     List<Integer> omitNonObserved = new ArrayList<>();
 464     int nonObservedShiftIndex = 0,pdbeNonObserved=0;
 465     // System.out.println("Generating mappings for : " + entityId);
 466     Entity entity = null;
 467     entity = getEntityById(entityId);
 468     String originalSeq = AlignSeq.extractGaps(
 469             jalview.util.Comparison.GapChars, seq.getSequenceAsString());
 470     HashMap<Integer, int[]> mapping = new HashMap<Integer, int[]>();
 471     DBRefEntryI sourceDBRef;
 472     sourceDBRef = getValidSourceDBRef(seq);
 473     // TODO ensure sequence start/end is in the same coordinate system and
 474     // consistent with the choosen sourceDBRef
 475
 476     // set sequence coordinate system - default value is UniProt
 477     if (sourceDBRef.getSource().equalsIgnoreCase(DBRefSource.PDB))
 478     {
 479       seqCoordSys = CoordinateSys.PDB;
 480     }
 481
 482     HashSet<String> dbRefAccessionIdsString = new HashSet<String>();
 483     for (DBRefEntry dbref : seq.getDBRefs())
 484     {
 485       dbRefAccessionIdsString.add(dbref.getAccessionId().toLowerCase());
 486     }
 487     dbRefAccessionIdsString.add(sourceDBRef.getAccessionId().toLowerCase());
 488
 489     curDBRefAccessionIdsString = dbRefAccessionIdsString;
 490     curSourceDBRef = sourceDBRef.getAccessionId();
 491
 492     TreeMap<Integer, String> resNumMap = new TreeMap<Integer, String>();
 493     List<Segment> segments = entity.getSegment();
 494     SegmentHelperPojo shp = new SegmentHelperPojo(seq, mapping, resNumMap,
 495             omitNonObserved, nonObservedShiftIndex,pdbeNonObserved);
 496     processSegments(segments, shp);
 497     try
 498     {
 499       populateAtomPositions(entityId, mapping);
 500     } catch (Exception e)
 501     {
 502       e.printStackTrace();
 503     }
 504     if (seqCoordSys == CoordinateSys.UNIPROT)
 505     {
 506       padWithGaps(resNumMap, omitNonObserved);
 507     }
 508     int seqStart = UNASSIGNED;
 509     int seqEnd = UNASSIGNED;
 510     int pdbStart = UNASSIGNED;
 511     int pdbEnd = UNASSIGNED;
 512
 513     if (mapping.isEmpty())
 514     {
 515       throw new SiftsException("SIFTS mapping failed");
 516     }
 517     // also construct a mapping object between the seq-coord sys and the PDB seq's coord sys
 518
 519     Integer[] keys = mapping.keySet().toArray(new Integer[0]);
 520     Arrays.sort(keys);
 521     seqStart = keys[0];
 522     seqEnd = keys[keys.length - 1];
 523     List<int[]> from=new ArrayList<>(),to=new ArrayList<>();
 524     int[]_cfrom=null,_cto=null;
 525     String matchedSeq = originalSeq;
 526     if (seqStart != UNASSIGNED) // fixme! seqStart can map to -1 for a pdb sequence that starts <-1
 527     {
 528       for (int seqps:keys)
 529       {
 530         int pdbpos = mapping.get(seqps)[PDBE_POS];
 531         if (pdbpos == UNASSIGNED)
 532         {
 533           // not correct - pdbpos might be -1, but leave it for now
 534           continue;
 535         }
 536         if (_cfrom==null || seqps!=_cfrom[1]+1)
 537         {
 538           _cfrom = new int[] { seqps,seqps};
 539           from.add(_cfrom);
 540           _cto = null; // discontinuity
 541         } else {
 542           _cfrom[1]= seqps;
 543         }
 544         if (_cto==null || pdbpos!=1+_cto[1])
 545         {
 546           _cto = new int[] { pdbpos,pdbpos};
 547           to.add(_cto);
 548         } else {
 549           _cto[1] = pdbpos;
 550         }
 551       }
 552       _cfrom = new int[from.size() * 2];
 553       _cto = new int[to.size() * 2];
 554       int p = 0;
 555       for (int[] range : from)
 556       {
 557         _cfrom[p++] = range[0];
 558         _cfrom[p++] = range[1];
 559       }
 560       ;
 561       p = 0;
 562       for (int[] range : to)
 563       {
 564         _cto[p++] = range[0];
 565         _cto[p++] = range[1];
 566       }
 567       ;
 568
 569       seqFromPdbMapping = new jalview.datamodel.Mapping(null, _cto, _cfrom,
 570               1,
 571               1);
 572       pdbStart = mapping.get(seqStart)[PDB_RES_POS];
 573       pdbEnd = mapping.get(seqEnd)[PDB_RES_POS];
 574       int orignalSeqStart = seq.getStart();
 575       if (orignalSeqStart >= 1)
 576       {
 577         int subSeqStart = (seqStart >= orignalSeqStart)
 578                 ? seqStart - orignalSeqStart
 579                 : 0;
 580         int subSeqEnd = seqEnd - (orignalSeqStart - 1);
 581         subSeqEnd = originalSeq.length() < subSeqEnd ? originalSeq.length()
 582                 : subSeqEnd;
 583         matchedSeq = originalSeq.substring(subSeqStart, subSeqEnd);
 584       }
 585       else
 586       {
 587         matchedSeq = originalSeq.substring(1, originalSeq.length());
 588       }
 589     }
 590
 591     StringBuilder targetStrucSeqs = new StringBuilder();
 592     for (String res : resNumMap.values())
 593     {
 594       targetStrucSeqs.append(res);
 595     }
 596
 597     if (os != null)
 598     {
 599       MappingOutputPojo mop = new MappingOutputPojo();
 600       mop.setSeqStart(seqStart);
 601       mop.setSeqEnd(seqEnd);
 602       mop.setSeqName(seq.getName());
 603       mop.setSeqResidue(matchedSeq);
 604
 605       mop.setStrStart(pdbStart);
 606       mop.setStrEnd(pdbEnd);
 607       mop.setStrName(structId);
 608       mop.setStrResidue(targetStrucSeqs.toString());
 609
 610       mop.setType("pep");
 611       os.print(getMappingOutput(mop).toString());
 612       os.println();
 613     }
 614     return mapping;
 615   }
 616
 617   void processSegments(List<Segment> segments, SegmentHelperPojo shp)
 618   {
 619     SequenceI seq = shp.getSeq();
 620     HashMap<Integer, int[]> mapping = shp.getMapping();
 621     TreeMap<Integer, String> resNumMap = shp.getResNumMap();
 622     List<Integer> omitNonObserved = shp.getOmitNonObserved();
 623     int nonObservedShiftIndex = shp.getNonObservedShiftIndex();
 624     int pdbeNonObservedCount = shp.getPdbeNonObserved();
 625     int firstPDBResNum = UNASSIGNED;
 626     for (Segment segment : segments)
 627     {
 628       // System.out.println("Mapping segments : " + segment.getSegId() + "\\"s
 629       // + segStartEnd);
 630       List<Residue> residues = segment.getListResidue().getResidue();
 631       for (Residue residue : residues)
 632       {
 633         boolean isObserved = isResidueObserved(residue);
 634         int pdbeIndex = getLeadingIntegerValue(residue.getDbResNum(),
 635                 UNASSIGNED);
 636         int currSeqIndex = UNASSIGNED;
 637         List<CrossRefDb> cRefDbs = residue.getCrossRefDb();
 638         CrossRefDb pdbRefDb = null;
 639         for (CrossRefDb cRefDb : cRefDbs)
 640         {
 641           if (cRefDb.getDbSource().equalsIgnoreCase(DBRefSource.PDB))
 642           {
 643             pdbRefDb = cRefDb;
 644             if (firstPDBResNum == UNASSIGNED)
 645             {
 646               firstPDBResNum = getLeadingIntegerValue(cRefDb.getDbResNum(),
 647                       UNASSIGNED);
 648             }
 649             else
 650             {
 651               if (isObserved)
 652               {
 653                 // after we find the first observed residue we just increment
 654                 firstPDBResNum++;
 655               }
 656             }
 657           }
 658           if (cRefDb.getDbCoordSys().equalsIgnoreCase(seqCoordSys.getName())
 659                   && isAccessionMatched(cRefDb.getDbAccessionId()))
 660           {
 661             currSeqIndex = getLeadingIntegerValue(cRefDb.getDbResNum(),
 662                     UNASSIGNED);
 663             if (pdbRefDb != null)
 664             {
 665               break;// exit loop if pdb and uniprot are already found
 666             }
 667           }
 668         }
 669         if (!isObserved)
 670         {
 671           ++pdbeNonObservedCount;
 672         }
 673         if (seqCoordSys == seqCoordSys.PDB) // FIXME: is seqCoordSys ever PDBe
 674                                             // ???
 675         {
 676           // if the sequence has a primary reference to the PDB, then we are
 677           // dealing with a sequence extracted directly from the PDB. In that
 678           // case, numbering is PDBe - non-observed residues
 679           currSeqIndex = seq.getStart() - 1 + pdbeIndex;
 680         }
 681         if (!isObserved)
 682         {
 683           if (seqCoordSys != CoordinateSys.UNIPROT) // FIXME: PDB or PDBe only
 684                                                     // here
 685           {
 686             // mapping to PDB or PDBe so we need to bookkeep for the
 687             // non-observed
 688             // SEQRES positions
 689             omitNonObserved.add(currSeqIndex);
 690             ++nonObservedShiftIndex;
 691           }
 692         }
 693         // if (currSeqIndex >= seq.getStart() && currSeqIndex <= seqlength) //
 694         // true
 695                                                                          // numbering
 696                                                                          // is
 697                                                                          // not
 698                                                                          // up
 699                                                                          // to
 700                                                                          // seq.getEnd()
 701         {
 702
 703           int resNum = (pdbRefDb == null)
 704                   ? getLeadingIntegerValue(residue.getDbResNum(),
 705                           UNASSIGNED)
 706                   : getLeadingIntegerValue(pdbRefDb.getDbResNum(),
 707                           UNASSIGNED);
 708
 709           if (isResidueObserved(residue)
 710                   || seqCoordSys == CoordinateSys.UNIPROT)
 711           {
 712             char resCharCode = ResidueProperties
 713                     .getSingleCharacterCode(ResidueProperties
 714                             .getCanonicalAminoAcid(residue.getDbResName()));
 715             resNumMap.put(currSeqIndex, String.valueOf(resCharCode));
 716           }
 717           // TODO: mapping should have SeqCoordSysNum => int[] {PDBeNum,
 718           // PDBRESNUM, ATOMNUM }
 719           mapping.put(currSeqIndex - nonObservedShiftIndex,
 720                   new int[]
 721                   { Integer.valueOf(resNum), UNASSIGNED,
 722                       isObserved ? firstPDBResNum
 723                               : UNASSIGNED });
 724         }
 725       }
 726     }
 727   }
 728
 729   /**
 730    * Get the leading integer part of a string that begins with an integer.
 731    *
 732    * @param input
 733    *          - the string input to process
 734    * @param failValue
 735    *          - value returned if unsuccessful
 736    * @return
 737    */
 738   static int getLeadingIntegerValue(String input, int failValue)
 739   {
 740     if (input == null)
 741     {
 742       return failValue;
 743     }
 744     String[] parts = input.split("(?=\\D)(?<=\\d)");
 745     if (parts != null && parts.length > 0 && parts[0].matches("[0-9]+"))
 746     {
 747       return Integer.valueOf(parts[0]);
 748     }
 749     return failValue;
 750   }
 751
 752   /**
 753    *
 754    * @param chainId
 755    *          Target chain to populate mapping of its atom positions.
 756    * @param mapping
 757    *          Two dimension array of residue index versus atom position
 758    * @throws IllegalArgumentException
 759    *           Thrown if chainId or mapping is null
 760    * @throws SiftsException
 761    */
 762   void populateAtomPositions(String chainId, Map<Integer, int[]> mapping)
 763           throws IllegalArgumentException, SiftsException
 764   {
 765     try
 766     {
 767       PDBChain chain = pdb.findChain(chainId);
 768
 769       if (chain == null || mapping == null)
 770       {
 771         throw new IllegalArgumentException(
 772                 "Chain id or mapping must not be null.");
 773       }
 774       for (int[] map : mapping.values())
 775       {
 776         if (map[PDB_RES_POS] != UNASSIGNED)
 777         {
 778           map[PDB_ATOM_POS] = getAtomIndex(map[PDB_RES_POS], chain.atoms);
 779         }
 780       }
 781     } catch (NullPointerException e)
 782     {
 783       throw new SiftsException(e.getMessage());
 784     } catch (Exception e)
 785     {
 786       throw new SiftsException(e.getMessage());
 787     }
 788   }
 789
 790   /**
 791    *
 792    * @param residueIndex
 793    *          The residue index used for the search
 794    * @param atoms
 795    *          A collection of Atom to search
 796    * @return atom position for the given residue index
 797    */
 798   int getAtomIndex(int residueIndex, Collection<Atom> atoms)
 799   {
 800     if (atoms == null)
 801     {
 802       throw new IllegalArgumentException(
 803               "atoms collection must not be null!");
 804     }
 805     for (Atom atom : atoms)
 806     {
 807       if (atom.resNumber == residueIndex)
 808       {
 809         return atom.atomIndex;
 810       }
 811     }
 812     return UNASSIGNED;
 813   }
 814
 815   /**
 816    * Checks if the residue instance is marked 'Not_observed' or not
 817    *
 818    * @param residue
 819    * @return
 820    */
 821   private boolean isResidueObserved(Residue residue)
 822   {
 823     Set<String> annotations = getResidueAnnotaitons(residue,
 824             ResidueDetailType.ANNOTATION);
 825     if (annotations == null || annotations.isEmpty())
 826     {
 827       return true;
 828     }
 829     for (String annotation : annotations)
 830     {
 831       if (annotation.equalsIgnoreCase(NOT_OBSERVED))
 832       {
 833         return false;
 834       }
 835     }
 836     return true;
 837   }
 838
 839   /**
 840    * Get annotation String for a given residue and annotation type
 841    *
 842    * @param residue
 843    * @param type
 844    * @return
 845    */
 846   private Set<String> getResidueAnnotaitons(Residue residue,
 847           ResidueDetailType type)
 848   {
 849     HashSet<String> foundAnnotations = new HashSet<String>();
 850     List<ResidueDetail> resDetails = residue.getResidueDetail();
 851     for (ResidueDetail resDetail : resDetails)
 852     {
 853       if (resDetail.getProperty().equalsIgnoreCase(type.getCode()))
 854       {
 855         foundAnnotations.add(resDetail.getContent());
 856       }
 857     }
 858     return foundAnnotations;
 859   }
 860
 861   @Override
 862   public boolean isAccessionMatched(String accession)
 863   {
 864     boolean isStrictMatch = true;
 865     return isStrictMatch ? curSourceDBRef.equalsIgnoreCase(accession)
 866             : curDBRefAccessionIdsString.contains(accession.toLowerCase());
 867   }
 868
 869   private boolean isFoundInSiftsEntry(String accessionId)
 870   {
 871     Set<String> siftsDBRefs = getAllMappingAccession();
 872     return accessionId != null
 873             && siftsDBRefs.contains(accessionId.toLowerCase());
 874   }
 875
 876   /**
 877    * Pad omitted residue positions in PDB sequence with gaps
 878    *
 879    * @param resNumMap
 880    */
 881   void padWithGaps(Map<Integer, String> resNumMap,
 882           List<Integer> omitNonObserved)
 883   {
 884     if (resNumMap == null || resNumMap.isEmpty())
 885     {
 886       return;
 887     }
 888     Integer[] keys = resNumMap.keySet().toArray(new Integer[0]);
 889     // Arrays.sort(keys);
 890     int firstIndex = keys[0];
 891     int lastIndex = keys[keys.length - 1];
 892     // System.out.println("Min value " + firstIndex);
 893     // System.out.println("Max value " + lastIndex);
 894     for (int x = firstIndex; x <= lastIndex; x++)
 895     {
 896       if (!resNumMap.containsKey(x) && !omitNonObserved.contains(x))
 897       {
 898         resNumMap.put(x, "-");
 899       }
 900     }
 901   }
 902
 903   @Override
 904   public Entity getEntityById(String id) throws SiftsException
 905   {
 906     // Determines an entity to process by performing a heuristic matching of all
 907     // Entities with the given chainId and choosing the best matching Entity
 908     Entity entity = getEntityByMostOptimalMatchedId(id);
 909     if (entity != null)
 910     {
 911       return entity;
 912     }
 913     throw new SiftsException("Entity " + id + " not found");
 914   }
 915
 916   /**
 917    * This method was added because EntityId is NOT always equal to ChainId.
 918    * Hence, it provides the logic to greedily detect the "true" Entity for a
 919    * given chainId where discrepancies exist.
 920    *
 921    * @param chainId
 922    * @return
 923    */
 924   public Entity getEntityByMostOptimalMatchedId(String chainId)
 925   {
 926     // System.out.println("---> advanced greedy entityId matching block
 927     // entered..");
 928     List<Entity> entities = siftsEntry.getEntity();
 929     SiftsEntitySortPojo[] sPojo = new SiftsEntitySortPojo[entities.size()];
 930     int count = 0;
 931     for (Entity entity : entities)
 932     {
 933       sPojo[count] = new SiftsEntitySortPojo();
 934       sPojo[count].entityId = entity.getEntityId();
 935
 936       List<Segment> segments = entity.getSegment();
 937       for (Segment segment : segments)
 938       {
 939         List<Residue> residues = segment.getListResidue().getResidue();
 940         for (Residue residue : residues)
 941         {
 942           List<CrossRefDb> cRefDbs = residue.getCrossRefDb();
 943           for (CrossRefDb cRefDb : cRefDbs)
 944           {
 945             if (!cRefDb.getDbSource().equalsIgnoreCase("PDB"))
 946             {
 947               continue;
 948             }
 949             ++sPojo[count].resCount;
 950             if (cRefDb.getDbChainId().equalsIgnoreCase(chainId))
 951             {
 952               ++sPojo[count].chainIdFreq;
 953             }
 954           }
 955         }
 956       }
 957       sPojo[count].pid = (100 * sPojo[count].chainIdFreq)
 958               / sPojo[count].resCount;
 959       ++count;
 960     }
 961     Arrays.sort(sPojo, Collections.reverseOrder());
 962     // System.out.println("highest matched entity : " + sPojo[0].entityId);
 963     // System.out.println("highest matched pid : " + sPojo[0].pid);
 964
 965     if (sPojo[0].entityId != null)
 966     {
 967       if (sPojo[0].pid < 1)
 968       {
 969         return null;
 970       }
 971       for (Entity entity : entities)
 972       {
 973         if (!entity.getEntityId().equalsIgnoreCase(sPojo[0].entityId))
 974         {
 975           continue;
 976         }
 977         return entity;
 978       }
 979     }
 980     return null;
 981   }
 982
 983   private class SiftsEntitySortPojo
 984           implements Comparable<SiftsEntitySortPojo>
 985   {
 986     public String entityId;
 987
 988     public int chainIdFreq;
 989
 990     public int pid;
 991
 992     public int resCount;
 993
 994     @Override
 995     public int compareTo(SiftsEntitySortPojo o)
 996     {
 997       return this.pid - o.pid;
 998     }
 999   }
1000
1001   private class SegmentHelperPojo
1002   {
1003     private SequenceI seq;
1004
1005     private HashMap<Integer, int[]> mapping;
1006
1007     private TreeMap<Integer, String> resNumMap;
1008
1009     private List<Integer> omitNonObserved;
1010
1011     private int nonObservedShiftIndex;
1012
1013     /**
1014      * count of number of 'not observed' positions in the PDB record's SEQRES
1015      * (total number of residues with coordinates == length(SEQRES) -
1016      * pdbeNonObserved
1017      */
1018     private int pdbeNonObserved;
1019
1020     public SegmentHelperPojo(SequenceI seq, HashMap<Integer, int[]> mapping,
1021             TreeMap<Integer, String> resNumMap,
1022             List<Integer> omitNonObserved, int nonObservedShiftIndex,
1023             int pdbeNonObserved)
1024     {
1025       setSeq(seq);
1026       setMapping(mapping);
1027       setResNumMap(resNumMap);
1028       setOmitNonObserved(omitNonObserved);
1029       setNonObservedShiftIndex(nonObservedShiftIndex);
1030       setPdbeNonObserved(pdbeNonObserved);
1031
1032     }
1033
1034     public void setPdbeNonObserved(int pdbeNonObserved2)
1035     {
1036       this.pdbeNonObserved = pdbeNonObserved2;
1037     }
1038
1039     public int getPdbeNonObserved()
1040     {
1041       return pdbeNonObserved;
1042     }
1043     public SequenceI getSeq()
1044     {
1045       return seq;
1046     }
1047
1048     public void setSeq(SequenceI seq)
1049     {
1050       this.seq = seq;
1051     }
1052
1053     public HashMap<Integer, int[]> getMapping()
1054     {
1055       return mapping;
1056     }
1057
1058     public void setMapping(HashMap<Integer, int[]> mapping)
1059     {
1060       this.mapping = mapping;
1061     }
1062
1063     public TreeMap<Integer, String> getResNumMap()
1064     {
1065       return resNumMap;
1066     }
1067
1068     public void setResNumMap(TreeMap<Integer, String> resNumMap)
1069     {
1070       this.resNumMap = resNumMap;
1071     }
1072
1073     public List<Integer> getOmitNonObserved()
1074     {
1075       return omitNonObserved;
1076     }
1077
1078     public void setOmitNonObserved(List<Integer> omitNonObserved)
1079     {
1080       this.omitNonObserved = omitNonObserved;
1081     }
1082
1083     public int getNonObservedShiftIndex()
1084     {
1085       return nonObservedShiftIndex;
1086     }
1087
1088     public void setNonObservedShiftIndex(int nonObservedShiftIndex)
1089     {
1090       this.nonObservedShiftIndex = nonObservedShiftIndex;
1091     }
1092
1093   }
1094
1095   @Override
1096   public StringBuilder getMappingOutput(MappingOutputPojo mp)
1097           throws SiftsException
1098   {
1099     String seqRes = mp.getSeqResidue();
1100     String seqName = mp.getSeqName();
1101     int sStart = mp.getSeqStart();
1102     int sEnd = mp.getSeqEnd();
1103
1104     String strRes = mp.getStrResidue();
1105     String strName = mp.getStrName();
1106     int pdbStart = mp.getStrStart();
1107     int pdbEnd = mp.getStrEnd();
1108
1109     String type = mp.getType();
1110
1111     int maxid = (seqName.length() >= strName.length()) ? seqName.length()
1112             : strName.length();
1113     int len = 72 - maxid - 1;
1114
1115     int nochunks = ((seqRes.length()) / len)
1116             + ((seqRes.length()) % len > 0 ? 1 : 0);
1117     // output mappings
1118     StringBuilder output = new StringBuilder(512);
1119     output.append(NEWLINE);
1120     output.append("Sequence \u27f7 Structure mapping details")
1121             .append(NEWLINE);
1122     output.append("Method: SIFTS");
1123     output.append(NEWLINE).append(NEWLINE);
1124
1125     output.append(new Format("%" + maxid + "s").form(seqName));
1126     output.append(" :  ");
1127     output.append(String.valueOf(sStart));
1128     output.append(" - ");
1129     output.append(String.valueOf(sEnd));
1130     output.append(" Maps to ");
1131     output.append(NEWLINE);
1132     output.append(new Format("%" + maxid + "s").form(structId));
1133     output.append(" :  ");
1134     output.append(String.valueOf(pdbStart));
1135     output.append(" - ");
1136     output.append(String.valueOf(pdbEnd));
1137     output.append(NEWLINE).append(NEWLINE);
1138
1139     ScoreMatrix pam250 = ScoreModels.getInstance().getPam250();
1140     int matchedSeqCount = 0;
1141     for (int j = 0; j < nochunks; j++)
1142     {
1143       // Print the first aligned sequence
1144       output.append(new Format("%" + (maxid) + "s").form(seqName))
1145               .append(" ");
1146
1147       for (int i = 0; i < len; i++)
1148       {
1149         if ((i + (j * len)) < seqRes.length())
1150         {
1151           output.append(seqRes.charAt(i + (j * len)));
1152         }
1153       }
1154
1155       output.append(NEWLINE);
1156       output.append(new Format("%" + (maxid) + "s").form(" ")).append(" ");
1157
1158       /*
1159        * Print out the match symbols:
1160        * | for exact match (ignoring case)
1161        * . if PAM250 score is positive
1162        * else a space
1163        */
1164       for (int i = 0; i < len; i++)
1165       {
1166         try
1167         {
1168           if ((i + (j * len)) < seqRes.length())
1169           {
1170             char c1 = seqRes.charAt(i + (j * len));
1171             char c2 = strRes.charAt(i + (j * len));
1172             boolean sameChar = Comparison.isSameResidue(c1, c2, false);
1173             if (sameChar && !Comparison.isGap(c1))
1174             {
1175               matchedSeqCount++;
1176               output.append("|");
1177             }
1178             else if (type.equals("pep"))
1179             {
1180               if (pam250.getPairwiseScore(c1, c2) > 0)
1181               {
1182                 output.append(".");
1183               }
1184               else
1185               {
1186                 output.append(" ");
1187               }
1188             }
1189             else
1190             {
1191               output.append(" ");
1192             }
1193           }
1194         } catch (IndexOutOfBoundsException e)
1195         {
1196           continue;
1197         }
1198       }
1199       // Now print the second aligned sequence
1200       output = output.append(NEWLINE);
1201       output = output.append(new Format("%" + (maxid) + "s").form(strName))
1202               .append(" ");
1203       for (int i = 0; i < len; i++)
1204       {
1205         if ((i + (j * len)) < strRes.length())
1206         {
1207           output.append(strRes.charAt(i + (j * len)));
1208         }
1209       }
1210       output.append(NEWLINE).append(NEWLINE);
1211     }
1212     float pid = (float) matchedSeqCount / seqRes.length() * 100;
1213     if (pid < SiftsSettings.getFailSafePIDThreshold())
1214     {
1215       throw new SiftsException(">>> Low PID detected for SIFTs mapping...");
1216     }
1217     output.append("Length of alignment = " + seqRes.length())
1218             .append(NEWLINE);
1219     output.append(new Format("Percentage ID = %2.2f").form(pid));
1220     return output;
1221   }
1222
1223   @Override
1224   public int getEntityCount()
1225   {
1226     return siftsEntry.getEntity().size();
1227   }
1228
1229   @Override
1230   public String getDbAccessionId()
1231   {
1232     return siftsEntry.getDbAccessionId();
1233   }
1234
1235   @Override
1236   public String getDbCoordSys()
1237   {
1238     return siftsEntry.getDbCoordSys();
1239   }
1240
1241   @Override
1242   public String getDbSource()
1243   {
1244     return siftsEntry.getDbSource();
1245   }
1246
1247   @Override
1248   public String getDbVersion()
1249   {
1250     return siftsEntry.getDbVersion();
1251   }
1252
1253   public static void setMockSiftsFile(File file)
1254   {
1255     mockSiftsFile = file;
1256   }
1257
1258 }