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