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