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