src/jalview/ext/htsjdk/HtsContigDb.java

   1 package jalview.ext.htsjdk;
   2
   3 import htsjdk.samtools.SAMSequenceDictionary;
   4 import htsjdk.samtools.SAMSequenceRecord;
   5 import htsjdk.samtools.reference.ReferenceSequence;
   6 import htsjdk.samtools.reference.ReferenceSequenceFile;
   7 import htsjdk.samtools.reference.ReferenceSequenceFileFactory;
   8 import htsjdk.samtools.util.StringUtil;
   9 import jalview.datamodel.Sequence;
  10 import jalview.datamodel.SequenceI;
  11
  12 import java.io.File;
  13 import java.math.BigInteger;
  14 import java.security.MessageDigest;
  15 import java.security.NoSuchAlgorithmException;
  16 import java.util.ArrayList;
  17 import java.util.HashSet;
  18 import java.util.List;
  19 import java.util.Set;
  20
  21 /**
  22  * a source of sequence data accessed via the HTSJDK
  23  *
  24  * @author jprocter
  25  *
  26  */
  27 public class HtsContigDb
  28 {
  29
  30   private String name;
  31
  32   private File dbLocation;
  33
  34   private htsjdk.samtools.reference.ReferenceSequenceFile refFile = null;
  35
  36   public HtsContigDb(String name, File descriptor) throws Exception
  37   {
  38     if (descriptor.isFile())
  39     {
  40       this.name = name;
  41       dbLocation = descriptor;
  42     }
  43     initSource();
  44   }
  45
  46   private void initSource() throws Exception
  47   {
  48     if (refFile != null)
  49     {
  50       return;
  51     }
  52
  53     refFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(
  54             dbLocation, true);
  55     if (refFile == null || refFile.getSequenceDictionary() == null)
  56     {
  57       // refFile = initSequenceDictionaryFor(dbLocation);
  58     }
  59
  60   }
  61
  62
  63   SAMSequenceDictionary rrefDict = null;
  64   private ReferenceSequenceFile initSequenceDictionaryFor(File dbLocation2) throws Exception
  65   {
  66     rrefDict = getDictionary(dbLocation2, true);
  67     if (rrefDict != null)
  68     {
  69       ReferenceSequenceFile rrefFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(dbLocation2, true);
  70       return rrefFile;
  71     }
  72     return null;
  73   }
  74   /**
  75    * code below hacked out from picard ----
  76    *
  77    * picard/src/java/picard/sam/CreateSequenceDictionary.java
  78    * https://github.com/
  79    * broadinstitute/picard/commit/270580d3e28123496576f0b91b3433179bb5d876
  80    */
  81
  82
  83   /*
  84    * The MIT License
  85    *
  86    * Copyright (c) 2009 The Broad Institute
  87    *
  88    * Permission is hereby granted, free of charge, to any person obtaining a
  89    * copy of this software and associated documentation files (the "Software"),
  90    * to deal in the Software without restriction, including without limitation
  91    * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  92    * and/or sell copies of the Software, and to permit persons to whom the
  93    * Software is furnished to do so, subject to the following conditions:
  94    *
  95    * The above copyright notice and this permission notice shall be included in
  96    * all copies or substantial portions of the Software.
  97    *
  98    * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  99    * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 100    * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 101    * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 102    * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 103    * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 104    * DEALINGS IN THE SOFTWARE.
 105    */
 106   /**
 107    *
 108    * @param f
 109    * @param truncate
 110    * @return
 111    * @throws Exception
 112    */
 113   SAMSequenceDictionary getDictionary(File f, boolean truncate)
 114           throws Exception
 115   {
 116     if (md5 == null)
 117     {
 118       initCreateSequenceDictionary();
 119     }
 120     final ReferenceSequenceFile refSeqFile = ReferenceSequenceFileFactory
 121             .getReferenceSequenceFile(f, truncate);
 122     ReferenceSequence refSeq;
 123     List<SAMSequenceRecord> ret = new ArrayList<SAMSequenceRecord>();
 124     Set<String> sequenceNames = new HashSet<String>();
 125     for (int numSequences = 0; (refSeq = refSeqFile.nextSequence()) != null; ++numSequences)
 126     {
 127       if (sequenceNames.contains(refSeq.getName()))
 128       {
 129         throw new Exception(
 130                 "Sequence name appears more than once in reference: "
 131                         + refSeq.getName());
 132       }
 133       sequenceNames.add(refSeq.getName());
 134       ret.add(makeSequenceRecord(refSeq));
 135     }
 136     return new SAMSequenceDictionary(ret);
 137   }
 138
 139   public boolean isValid()
 140   {
 141     return dbLocation != null && refFile != null;
 142   }
 143
 144   /**
 145    * Create one SAMSequenceRecord from a single fasta sequence
 146    */
 147   private SAMSequenceRecord makeSequenceRecord(
 148           final ReferenceSequence refSeq)
 149   {
 150
 151     final SAMSequenceRecord ret = new SAMSequenceRecord(refSeq.getName(),
 152             refSeq.length());
 153
 154     // Compute MD5 of upcased bases
 155     final byte[] bases = refSeq.getBases();
 156     for (int i = 0; i < bases.length; ++i)
 157     {
 158       bases[i] = StringUtil.toUpperCase(bases[i]);
 159     }
 160
 161     ret.setAttribute(SAMSequenceRecord.MD5_TAG, md5Hash(bases));
 162     // if (GENOME_ASSEMBLY != null) {
 163     // ret.setAttribute(SAMSequenceRecord.ASSEMBLY_TAG, GENOME_ASSEMBLY);
 164     // }
 165     // ret.setAttribute(SAMSequenceRecord.URI_TAG, URI);
 166     // if (SPECIES != null) {
 167     // ret.setAttribute(SAMSequenceRecord.SPECIES_TAG, SPECIES);
 168     // }
 169     return ret;
 170   }
 171
 172   private MessageDigest md5;
 173
 174   public void initCreateSequenceDictionary() throws Exception
 175   {
 176     try
 177     {
 178       md5 = MessageDigest.getInstance("MD5");
 179     } catch (NoSuchAlgorithmException e)
 180     {
 181       throw new Exception("MD5 algorithm not found", e);
 182     }
 183   }
 184
 185   private String md5Hash(final byte[] bytes)
 186   {
 187     md5.reset();
 188     md5.update(bytes);
 189     String s = new BigInteger(1, md5.digest()).toString(16);
 190     if (s.length() != 32)
 191     {
 192       final String zeros = "00000000000000000000000000000000";
 193       s = zeros.substring(0, 32 - s.length()) + s;
 194     }
 195     return s;
 196   }
 197
 198   // ///// end of hts bits.
 199
 200   SequenceI getSequenceProxy(String id)
 201   {
 202     if (!isValid())
 203     {
 204       return null;
 205     }
 206
 207     ReferenceSequence sseq = refFile.getSequence(id);
 208     return new Sequence(sseq.getName(), new String(sseq.getBases()));
 209   }
 210 }