1 package jalview.ext.htsjdk;
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;
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;
22 * a source of sequence data accessed via the HTSJDK
27 public class HtsContigDb
32 private File dbLocation;
34 private htsjdk.samtools.reference.ReferenceSequenceFile refFile = null;
36 public HtsContigDb(String name, File descriptor) throws Exception
38 if (descriptor.isFile())
41 dbLocation = descriptor;
46 private void initSource() throws Exception
53 refFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(
55 if (refFile == null || refFile.getSequenceDictionary() == null)
57 // refFile = initSequenceDictionaryFor(dbLocation);
63 SAMSequenceDictionary rrefDict = null;
64 private ReferenceSequenceFile initSequenceDictionaryFor(File dbLocation2) throws Exception
66 rrefDict = getDictionary(dbLocation2, true);
69 ReferenceSequenceFile rrefFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(dbLocation2, true);
75 * code below hacked out from picard ----
77 * picard/src/java/picard/sam/CreateSequenceDictionary.java
79 * broadinstitute/picard/commit/270580d3e28123496576f0b91b3433179bb5d876
86 * Copyright (c) 2009 The Broad Institute
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:
95 * The above copyright notice and this permission notice shall be included in
96 * all copies or substantial portions of the Software.
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.
113 SAMSequenceDictionary getDictionary(File f, boolean truncate)
118 initCreateSequenceDictionary();
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)
127 if (sequenceNames.contains(refSeq.getName()))
130 "Sequence name appears more than once in reference: "
133 sequenceNames.add(refSeq.getName());
134 ret.add(makeSequenceRecord(refSeq));
136 return new SAMSequenceDictionary(ret);
139 public boolean isValid()
141 return dbLocation != null && refFile != null;
145 * Create one SAMSequenceRecord from a single fasta sequence
147 private SAMSequenceRecord makeSequenceRecord(
148 final ReferenceSequence refSeq)
151 final SAMSequenceRecord ret = new SAMSequenceRecord(refSeq.getName(),
154 // Compute MD5 of upcased bases
155 final byte[] bases = refSeq.getBases();
156 for (int i = 0; i < bases.length; ++i)
158 bases[i] = StringUtil.toUpperCase(bases[i]);
161 ret.setAttribute(SAMSequenceRecord.MD5_TAG, md5Hash(bases));
162 // if (GENOME_ASSEMBLY != null) {
163 // ret.setAttribute(SAMSequenceRecord.ASSEMBLY_TAG, GENOME_ASSEMBLY);
165 // ret.setAttribute(SAMSequenceRecord.URI_TAG, URI);
166 // if (SPECIES != null) {
167 // ret.setAttribute(SAMSequenceRecord.SPECIES_TAG, SPECIES);
172 private MessageDigest md5;
174 public void initCreateSequenceDictionary() throws Exception
178 md5 = MessageDigest.getInstance("MD5");
179 } catch (NoSuchAlgorithmException e)
181 throw new Exception("MD5 algorithm not found", e);
185 private String md5Hash(final byte[] bytes)
189 String s = new BigInteger(1, md5.digest()).toString(16);
190 if (s.length() != 32)
192 final String zeros = "00000000000000000000000000000000";
193 s = zeros.substring(0, 32 - s.length()) + s;
198 // ///// end of hts bits.
200 SequenceI getSequenceProxy(String id)
207 ReferenceSequence sseq = refFile.getSequence(id);
208 return new Sequence(sseq.getName(), new String(sseq.getBases()));