Classes

Class	Line #	Total Statements	Complexity	Uncovered Elements	TOTAL Coverage	Actions
VCFLoader	52	355	134	176	0.6752767667.5%	Show methods
VCFLoader.VCFMap	64	3	2	2	0.660%	Show methods

Class VCFLoader

Class VCFLoader	Line # 52	Total Statements 355	Complexity 134	Uncovered Elements 176	TOTAL Coverage 0.6752767667.5%
VCFLoader(String)   VCFLoader(String)	199199	4.04	2.02	1.01	0.75 0.7575%
run() : void   run() : void	232232	1.01	1.01	1.01	0.0 0.00%
loadVCF(SequenceI[],AlignViewControllerGuiI) : void   loadVCF(SequenceI[],AlignViewControllerGuiI) : void	223223	3.03	2.02	5.05	0.0 0.00%
loadVCFContig(String) : SequenceI   loadVCFContig(String) : SequenceI	247247	12.012	3.03	4.04	0.75 0.7575%
doLoad(SequenceI[],AlignViewControllerGuiI) : void   doLoad(SequenceI[],AlignViewControllerGuiI) : void	281281	25.025	8.08	14.014	0.6 0.660%
initialise(String) : void   initialise(String) : void	346346	8.08	2.02	0.00	1.0 1.0100%
saveMetadata(String) : void   saveMetadata(String) : void	378378	29.029	7.07	7.07	0.7741935 0.774193577.4%
isFieldWanted(String,List<Pattern>) : boolean   isFieldWanted(String,List<Pattern>) : boolean	431431	4.04	2.02	2.02	0.6666667 0.666666766.7%
parseCsqHeader() : void   parseCsqHeader() : void	452452	23.023	8.08	4.04	0.8918919 0.891891989.2%
getFieldMatchers(String,String) : List<Pattern>   getFieldMatchers(String,String) : List<Pattern>	520520	8.08	2.02	1.01	0.875 0.87587.5%
transferAddedFeatures(SequenceI) : void   transferAddedFeatures(SequenceI) : void	544544	15.015	6.06	9.09	0.6086956 0.608695660.9%
loadSequenceVCF(SequenceI,String) : int   loadSequenceVCF(SequenceI,String) : int	594594	7.07	3.03	2.02	0.8181818 0.818181881.8%
getVcfMap(SequenceI,String) : VCFMap   getVcfMap(SequenceI,String) : VCFMap	620620	33.033	10.010	28.028	0.42857143 0.4285714342.9%
getContigMap(SequenceI) : VCFMap   getContigMap(SequenceI) : VCFMap	717717	8.08	3.03	3.03	0.75 0.7575%
vcfAssemblyMatchesSequence(String,String) : boolean   vcfAssemblyMatchesSequence(String,String) : boolean	743743	3.03	3.03	2.02	0.6 0.660%
vcfSpeciesMatchesSequence(String,String) : boolean   vcfSpeciesMatchesSequence(String,String) : boolean	765765	5.05	5.05	6.06	0.33333334 0.3333333433.3%
addVcfVariants(SequenceI,VCFMap) : int   addVcfVariants(SequenceI,VCFMap) : int	801801	15.015	3.03	0.00	1.0 1.0100%
getAlleleFrequency(VariantContext,int) : float   getAlleleFrequency(VariantContext,int) : float	845845	6.06	3.03	1.01	0.875 0.87587.5%
getAttributeValue(VariantContext,String,int) : String   getAttributeValue(VariantContext,String,int) : String	872872	6.06	3.03	2.02	0.8 0.880%
addAlleleFeatures(SequenceI,VariantContext,int,int,boolean) : int   addAlleleFeatures(SequenceI,VariantContext,int,int,boolean) : int	900900	5.05	2.02	0.00	1.0 1.0100%
addAlleleFeature(SequenceI,VariantContext,int,int,int,boolean) : int   addAlleleFeature(SequenceI,VariantContext,int,int,int,boolean) : int	933933	26.026	7.07	0.00	1.0 1.0100%
getOntologyTerm(String) : String   getOntologyTerm(String) : String	10121012	12.012	6.06	7.07	0.6818182 0.681818268.2%
getConsequenceForAlleleAndFeature(VariantContext,String,int,int,int,String,int) : String   getConsequenceForAlleleAndFeature(VariantContext,String,int,int,int,String,int) : String	10781078	14.014	9.09	3.03	0.875 0.87587.5%
matchAllele(VariantContext,int,String[],int,int) : boolean   matchAllele(VariantContext,int,String[],int,int) : boolean	11231123	10.010	5.05	9.09	0.4375 0.437543.8%
addAlleleProperties(VariantContext,SequenceFeature,int,String) : void   addAlleleProperties(VariantContext,SequenceFeature,int,String) : void	11651165	22.022	8.08	5.05	0.8611111 0.861111186.1%
addConsequences(VariantContext,SequenceFeature,String) : void   addConsequences(VariantContext,SequenceFeature,String) : void	12561256	17.017	9.09	5.05	0.8148148 0.814814881.5%
complement(byte[]) : String   complement(byte[]) : String	13131313	1.01	1.01	1.01	0.0 0.00%
mapReferenceRange(int[],String,String,String,String) : int[]   mapReferenceRange(int[],String,String,String,String) : int[]	13391339	12.012	4.04	18.018	0.0 0.00%
findSubsumedRangeMapping(int[],String,String,String,String) : int[]   findSubsumedRangeMapping(int[],String,String,String,String) : int[]	14011401	13.013	4.04	19.019	0.0 0.00%
transferFeature(SequenceFeature,SequenceI,MapList) : void   transferFeature(SequenceFeature,SequenceI,MapList) : void	14431443	7.07	2.02	9.09	0.0 0.00%
makeRangesKey(String,String,String,String) : String   makeRangesKey(String,String,String,String) : String	14681468	1.01	1.01	1.01	0.0 0.00%

 

Class VCFLoader.VCFMap

Class VCFLoader.VCFMap	Line # 64	Total Statements 3	Complexity 2	Uncovered Elements 2	TOTAL Coverage 0.660%
VCFMap(String,MapList)   VCFMap(String,MapList)	7070	2.02	1.01	0.00	1.0 1.0100%
toString() : String   toString() : String	7676	1.01	1.01	1.01	0.0 0.00%

 

Contributing tests

Source view

1		package jalview.io.vcf;
2
3		import jalview.analysis.AlignmentUtils;
4		import jalview.analysis.Dna;
5		import jalview.api.AlignViewControllerGuiI;
6		import jalview.bin.Cache;
7		import jalview.datamodel.DBRefEntry;
8		import jalview.datamodel.GeneLociI;
9		import jalview.datamodel.Mapping;
10		import jalview.datamodel.SequenceFeature;
11		import jalview.datamodel.SequenceI;
12		import jalview.datamodel.features.FeatureAttributeType;
13		import jalview.datamodel.features.FeatureSource;
14		import jalview.datamodel.features.FeatureSources;
15		import jalview.ext.ensembl.EnsemblMap;
16		import jalview.ext.htsjdk.HtsContigDb;
17		import jalview.ext.htsjdk.VCFReader;
18		import jalview.io.gff.Gff3Helper;
19		import jalview.io.gff.SequenceOntologyI;
20		import jalview.util.MapList;
21		import jalview.util.MappingUtils;
22		import jalview.util.MessageManager;
23
24		import java.io.File;
25		import java.io.IOException;
26		import java.util.ArrayList;
27		import java.util.HashMap;
28		import java.util.List;
29		import java.util.Map;
30		import java.util.Map.Entry;
31		import java.util.regex.Pattern;
32		import java.util.regex.PatternSyntaxException;
33
34		import htsjdk.samtools.SAMException;
35		import htsjdk.samtools.SAMSequenceDictionary;
36		import htsjdk.samtools.SAMSequenceRecord;
37		import htsjdk.samtools.util.CloseableIterator;
38		import htsjdk.variant.variantcontext.Allele;
39		import htsjdk.variant.variantcontext.VariantContext;
40		import htsjdk.variant.vcf.VCFHeader;
41		import htsjdk.variant.vcf.VCFHeaderLine;
42		import htsjdk.variant.vcf.VCFHeaderLineCount;
43		import htsjdk.variant.vcf.VCFHeaderLineType;
44		import htsjdk.variant.vcf.VCFInfoHeaderLine;
45
46		/**
47		* A class to read VCF data (using the htsjdk) and add variants as sequence
48		* features on dna and any related protein product sequences
49		*
50		* @author gmcarstairs
51		*/
		67.5% Uncovered Elements: 176 (542) Complexity: 134 Complexity Density: 0.38
52		public class VCFLoader
53		{
54		/**
55		* A class to model the mapping from sequence to VCF coordinates. Cases include
56		* <ul>
57		* <li>a direct 1:1 mapping where the sequence is one of the VCF contigs</li>
58		* <li>a mapping of sequence to chromosomal coordinates, where sequence and VCF
59		* use the same reference assembly</li>
60		* <li>a modified mapping of sequence to chromosomal coordinates, where sequence
61		* and VCF use different reference assembles</li>
62		* </ul>
63		*/
		60% Uncovered Elements: 2 (5) Complexity: 2 Complexity Density: 0.67
64		class VCFMap
65		{
66		final String chromosome;
67
68		final MapList map;
69
		100% Uncovered Elements: 0 (2) Complexity: 1 Complexity Density: 0.5
70	24	VCFMap(String chr, MapList m)...
71		{
72	24	chromosome = chr;
73	24	map = m;
74		}
75
		0% Uncovered Elements: 1 (1) Complexity: 1 Complexity Density: 1
76	0	@Override...
77		public String toString()
78		{
79	0	return chromosome + ":" + map.toString();
80		}
81		}
82
83		/*
84		* Lookup keys, and default values, for Preference entries that describe
85		* patterns for VCF and VEP fields to capture
86		*/
87		private static final String VEP_FIELDS_PREF = "VEP_FIELDS";
88
89		private static final String VCF_FIELDS_PREF = "VCF_FIELDS";
90
91		private static final String DEFAULT_VCF_FIELDS = ".*";
92
93		private static final String DEFAULT_VEP_FIELDS = ".*";// "Allele,Consequence,IMPACT,SWISSPROT,SIFT,PolyPhen,CLIN_SIG";
94
95		/*
96		* keys to fields of VEP CSQ consequence data
97		* see https://www.ensembl.org/info/docs/tools/vep/vep_formats.html
98		*/
99		private static final String CSQ_CONSEQUENCE_KEY = "Consequence";
100		private static final String CSQ_ALLELE_KEY = "Allele";
101		private static final String CSQ_ALLELE_NUM_KEY = "ALLELE_NUM"; // 0 (ref), 1...
102		private static final String CSQ_FEATURE_KEY = "Feature"; // Ensembl stable id
103
104		/*
105		* default VCF INFO key for VEP consequence data
106		* NB this can be overridden running VEP with --vcf_info_field
107		* - we don't handle this case (require identifier to be CSQ)
108		*/
109		private static final String CSQ_FIELD = "CSQ";
110
111		/*
112		* separator for fields in consequence data is '|'
113		*/
114		private static final String PIPE_REGEX = "\\|";
115
116		/*
117		* key for Allele Frequency output by VEP
118		* see http://www.ensembl.org/info/docs/tools/vep/vep_formats.html
119		*/
120		private static final String ALLELE_FREQUENCY_KEY = "AF";
121
122		/*
123		* delimiter that separates multiple consequence data blocks
124		*/
125		private static final String COMMA = ",";
126
127		/*
128		* the feature group assigned to a VCF variant in Jalview
129		*/
130		private static final String FEATURE_GROUP_VCF = "VCF";
131
132		/*
133		* internal delimiter used to build keys for assemblyMappings
134		*
135		*/
136		private static final String EXCL = "!";
137
138		/*
139		* the VCF file we are processing
140		*/
141		protected String vcfFilePath;
142
143		/*
144		* mappings between VCF and sequence reference assembly regions, as
145		* key = "species!chromosome!fromAssembly!toAssembly
146		* value = Map{fromRange, toRange}
147		*/
148		private Map<String, Map<int[], int[]>> assemblyMappings;
149
150		private VCFReader reader;
151
152		/*
153		* holds details of the VCF header lines (metadata)
154		*/
155		private VCFHeader header;
156
157		/*
158		* a Dictionary of contigs (if present) referenced in the VCF file
159		*/
160		private SAMSequenceDictionary dictionary;
161
162		/*
163		* the position (0...) of field in each block of
164		* CSQ (consequence) data (if declared in the VCF INFO header for CSQ)
165		* see http://www.ensembl.org/info/docs/tools/vep/vep_formats.html
166		*/
167		private int csqConsequenceFieldIndex = -1;
168		private int csqAlleleFieldIndex = -1;
169		private int csqAlleleNumberFieldIndex = -1;
170		private int csqFeatureFieldIndex = -1;
171
172		// todo the same fields for SnpEff ANN data if wanted
173		// see http://snpeff.sourceforge.net/SnpEff_manual.html#input
174
175		/*
176		* a unique identifier under which to save metadata about feature
177		* attributes (selected INFO field data)
178		*/
179		private String sourceId;
180
181		/*
182		* The INFO IDs of data that is both present in the VCF file, and
183		* also matched by any filters for data of interest
184		*/
185		List<String> vcfFieldsOfInterest;
186
187		/*
188		* The field offsets and identifiers for VEP (CSQ) data that is both present
189		* in the VCF file, and also matched by any filters for data of interest
190		* for example 0 -> Allele, 1 -> Consequence, ..., 36 -> SIFT, ...
191		*/
192		Map<Integer, String> vepFieldsOfInterest;
193
194		/**
195		* Constructor given a VCF file
196		*
197		* @param alignment
198		*/
		75% Uncovered Elements: 1 (4) Complexity: 2 Complexity Density: 0.5
199	4	public VCFLoader(String vcfFile)...
200		{
201	4	try
202		{
203	4	initialise(vcfFile);
204		} catch (IOException e)
205		{
206	0	System.err.println("Error opening VCF file: " + e.getMessage());
207		}
208
209		// map of species!chromosome!fromAssembly!toAssembly to {fromRange, toRange}
210	4	assemblyMappings = new HashMap<>();
211		}
212
213		/**
214		* Starts a new thread to query and load VCF variant data on to the given
215		* sequences
216		* <p>
217		* This method is not thread safe - concurrent threads should use separate
218		* instances of this class.
219		*
220		* @param seqs
221		* @param gui
222		*/
		0% Uncovered Elements: 5 (5) Complexity: 2 Complexity Density: 0.67
223	0	public void loadVCF(SequenceI[] seqs, final AlignViewControllerGuiI gui)...
224		{
225	0	if (gui != null)
226		{
227	0	gui.setStatus(MessageManager.getString("label.searching_vcf"));
228		}
229
230	0	new Thread()
231		{
		0% Uncovered Elements: 1 (1) Complexity: 1 Complexity Density: 1
232	0	@Override...
233		public void run()
234		{
235	0	VCFLoader.this.doLoad(seqs, gui);
236		}
237		}.start();
238		}
239
240		/**
241		* Reads the specified contig sequence and adds its VCF variants to it
242		*
243		* @param contig
244		* the id of a single sequence (contig) to load
245		* @return
246		*/
		75% Uncovered Elements: 4 (16) Complexity: 3 Complexity Density: 0.25
247	3	public SequenceI loadVCFContig(String contig)...
248		{
249	3	String ref = header.getOtherHeaderLine(VCFHeader.REFERENCE_KEY)
250		.getValue();
251	3	if (ref.startsWith("file://"))
252		{
253	0	ref = ref.substring(7);
254		}
255
256	3	SequenceI seq = null;
257	3	File dbFile = new File(ref);
258
259	3	if (dbFile.exists())
260		{
261	3	HtsContigDb db = new HtsContigDb("", dbFile);
262	3	seq = db.getSequenceProxy(contig);
263	3	loadSequenceVCF(seq, ref);
264	3	db.close();
265		}
266		else
267		{
268	0	System.err.println("VCF reference not found: " + ref);
269		}
270
271	3	return seq;
272		}
273
274		/**
275		* Loads VCF on to one or more sequences
276		*
277		* @param seqs
278		* @param gui
279		* optional callback handler for messages
280		*/
		60% Uncovered Elements: 14 (35) Complexity: 8 Complexity Density: 0.32
281	3	protected void doLoad(SequenceI[] seqs, AlignViewControllerGuiI gui)...
282		{
283	3	try
284		{
285	3	VCFHeaderLine ref = header
286		.getOtherHeaderLine(VCFHeader.REFERENCE_KEY);
287	3	String vcfAssembly = ref.getValue();
288
289	3	int varCount = 0;
290	3	int seqCount = 0;
291
292		/*
293		* query for VCF overlapping each sequence in turn
294		*/
295	3	for (SequenceI seq : seqs)
296		{
297	21	int added = loadSequenceVCF(seq, vcfAssembly);
298	21	if (added > 0)
299		{
300	11	seqCount++;
301	11	varCount += added;
302	11	transferAddedFeatures(seq);
303		}
304		}
305	3	if (gui != null)
306		{
307	0	String msg = MessageManager.formatMessage("label.added_vcf",
308		varCount, seqCount);
309	0	gui.setStatus(msg);
310	0	if (gui.getFeatureSettingsUI() != null)
311		{
312	0	gui.getFeatureSettingsUI().discoverAllFeatureData();
313		}
314		}
315		} catch (Throwable e)
316		{
317	0	System.err.println("Error processing VCF: " + e.getMessage());
318	0	e.printStackTrace();
319	0	if (gui != null)
320		{
321	0	gui.setStatus("Error occurred - see console for details");
322		}
323		} finally
324		{
325	3	if (reader != null)
326		{
327	3	try
328		{
329	3	reader.close();
330		} catch (IOException e)
331		{
332		// ignore
333		}
334		}
335	3	header = null;
336	3	dictionary = null;
337		}
338		}
339
340		/**
341		* Opens the VCF file and parses header data
342		*
343		* @param filePath
344		* @throws IOException
345		*/
		100% Uncovered Elements: 0 (8) Complexity: 2 Complexity Density: 0.25
346	4	private void initialise(String filePath) throws IOException...
347		{
348	4	vcfFilePath = filePath;
349
350	4	reader = new VCFReader(filePath);
351
352	4	header = reader.getFileHeader();
353
354	4	try
355		{
356	4	dictionary = header.getSequenceDictionary();
357		} catch (SAMException e)
358		{
359		// ignore - thrown if any contig line lacks length info
360		}
361
362	4	sourceId = filePath;
363
364	4	saveMetadata(sourceId);
365
366		/*
367		* get offset of CSQ ALLELE_NUM and Feature if declared
368		*/
369	4	parseCsqHeader();
370		}
371
372		/**
373		* Reads metadata (such as INFO field descriptions and datatypes) and saves
374		* them for future reference
375		*
376		* @param theSourceId
377		*/
		77.4% Uncovered Elements: 7 (31) Complexity: 7 Complexity Density: 0.24
378	4	void saveMetadata(String theSourceId)...
379		{
380	4	List<Pattern> vcfFieldPatterns = getFieldMatchers(VCF_FIELDS_PREF,
381		DEFAULT_VCF_FIELDS);
382	4	vcfFieldsOfInterest = new ArrayList<>();
383
384	4	FeatureSource metadata = new FeatureSource(theSourceId);
385
386	4	for (VCFInfoHeaderLine info : header.getInfoHeaderLines())
387		{
388	9	String attributeId = info.getID();
389	9	String desc = info.getDescription();
390	9	VCFHeaderLineType type = info.getType();
391	9	FeatureAttributeType attType = null;
392	9	switch (type)
393		{
394	0	case Character:
395	0	attType = FeatureAttributeType.Character;
396	0	break;
397	0	case Flag:
398	0	attType = FeatureAttributeType.Flag;
399	0	break;
400	5	case Float:
401	5	attType = FeatureAttributeType.Float;
402	5	break;
403	3	case Integer:
404	3	attType = FeatureAttributeType.Integer;
405	3	break;
406	1	case String:
407	1	attType = FeatureAttributeType.String;
408	1	break;
409		}
410	9	metadata.setAttributeName(attributeId, desc);
411	9	metadata.setAttributeType(attributeId, attType);
412
413	9	if (isFieldWanted(attributeId, vcfFieldPatterns))
414		{
415	9	vcfFieldsOfInterest.add(attributeId);
416		}
417		}
418
419	4	FeatureSources.getInstance().addSource(theSourceId, metadata);
420		}
421
422		/**
423		* Answers true if the field id is matched by any of the filter patterns, else
424		* false. Matching is against regular expression patterns, and is not
425		* case-sensitive.
426		*
427		* @param id
428		* @param filters
429		* @return
430		*/
		66.7% Uncovered Elements: 2 (6) Complexity: 2 Complexity Density: 0.5
431	18	private boolean isFieldWanted(String id, List<Pattern> filters)...
432		{
433	18	for (Pattern p : filters)
434		{
435	18	if (p.matcher(id.toUpperCase()).matches())
436		{
437	18	return true;
438		}
439		}
440	0	return false;
441		}
442
443		/**
444		* Records 'wanted' fields defined in the CSQ INFO header (if there is one).
445		* Also records the position of selected fields (Allele, ALLELE_NUM, Feature)
446		* required for processing.
447		* <p>
448		* CSQ fields are declared in the CSQ INFO Description e.g.
449		* <p>
450		* Description="Consequence ...from ... VEP. Format: Allele|Consequence|...
451		*/
		89.2% Uncovered Elements: 4 (37) Complexity: 8 Complexity Density: 0.35
452	4	protected void parseCsqHeader()...
453		{
454	4	List<Pattern> vepFieldFilters = getFieldMatchers(VEP_FIELDS_PREF,
455		DEFAULT_VEP_FIELDS);
456	4	vepFieldsOfInterest = new HashMap<>();
457
458	4	VCFInfoHeaderLine csqInfo = header.getInfoHeaderLine(CSQ_FIELD);
459	4	if (csqInfo == null)
460		{
461	3	return;
462		}
463
464		/*
465		* parse out the pipe-separated list of CSQ fields; we assume here that
466		* these form the last part of the description, and contain no spaces
467		*/
468	1	String desc = csqInfo.getDescription();
469	1	int spacePos = desc.lastIndexOf(" ");
470	1	desc = desc.substring(spacePos + 1);
471
472	1	if (desc != null)
473		{
474	1	String[] format = desc.split(PIPE_REGEX);
475	1	int index = 0;
476	1	for (String field : format)
477		{
478	9	if (CSQ_CONSEQUENCE_KEY.equals(field))
479		{
480	1	csqConsequenceFieldIndex = index;
481		}
482	9	if (CSQ_ALLELE_NUM_KEY.equals(field))
483		{
484	0	csqAlleleNumberFieldIndex = index;
485		}
486	9	if (CSQ_ALLELE_KEY.equals(field))
487		{
488	1	csqAlleleFieldIndex = index;
489		}
490	9	if (CSQ_FEATURE_KEY.equals(field))
491		{
492	1	csqFeatureFieldIndex = index;
493		}
494
495	9	if (isFieldWanted(field, vepFieldFilters))
496		{
497	9	vepFieldsOfInterest.put(index, field);
498		}
499
500	9	index++;
501		}
502		}
503		}
504
505		/**
506		* Reads the Preference value for the given key, with default specified if no
507		* preference set. The value is interpreted as a comma-separated list of
508		* regular expressions, and converted into a list of compiled patterns ready
509		* for matching. Patterns are forced to upper-case for non-case-sensitive
510		* matching.
511		* <p>
512		* This supports user-defined filters for fields of interest to capture while
513		* processing data. For example, VCF_FIELDS = AF,AC* would mean that VCF INFO
514		* fields with an ID of AF, or starting with AC, would be matched.
515		*
516		* @param key
517		* @param def
518		* @return
519		*/
		87.5% Uncovered Elements: 1 (8) Complexity: 2 Complexity Density: 0.25
520	8	private List<Pattern> getFieldMatchers(String key, String def)...
521		{
522	8	String pref = Cache.getDefault(key, def);
523	8	List<Pattern> patterns = new ArrayList<>();
524	8	String[] tokens = pref.split(",");
525	8	for (String token : tokens)
526		{
527	8	try
528		{
529	8	patterns.add(Pattern.compile(token.toUpperCase()));
530		} catch (PatternSyntaxException e)
531		{
532	0	System.err.println("Invalid pattern ignored: " + token);
533		}
534		}
535	8	return patterns;
536		}
537
538		/**
539		* Transfers VCF features to sequences to which this sequence has a mapping.
540		* If the mapping is 3:1, computes peptide variants from nucleotide variants.
541		*
542		* @param seq
543		*/
		60.9% Uncovered Elements: 9 (23) Complexity: 6 Complexity Density: 0.4
544	11	protected void transferAddedFeatures(SequenceI seq)...
545		{
546	11	DBRefEntry[] dbrefs = seq.getDBRefs();
547	11	if (dbrefs == null)
548		{
549	0	return;
550		}
551	11	for (DBRefEntry dbref : dbrefs)
552		{
553	16	Mapping mapping = dbref.getMap();
554	16	if (mapping == null || mapping.getTo() == null)
555		{
556	11	continue;
557		}
558
559	5	SequenceI mapTo = mapping.getTo();
560	5	MapList map = mapping.getMap();
561	5	if (map.getFromRatio() == 3)
562		{
563		/*
564		* dna-to-peptide product mapping
565		*/
566	5	AlignmentUtils.computeProteinFeatures(seq, mapTo, map);
567		}
568		else
569		{
570		/*
571		* nucleotide-to-nucleotide mapping e.g. transcript to CDS
572		*/
573	0	List<SequenceFeature> features = seq.getFeatures()
574		.getPositionalFeatures(SequenceOntologyI.SEQUENCE_VARIANT);
575	0	for (SequenceFeature sf : features)
576		{
577	0	if (FEATURE_GROUP_VCF.equals(sf.getFeatureGroup()))
578		{
579	0	transferFeature(sf, mapTo, map);
580		}
581		}
582		}
583		}
584		}
585
586		/**
587		* Tries to add overlapping variants read from a VCF file to the given sequence,
588		* and returns the number of variant features added
589		*
590		* @param seq
591		* @param vcfAssembly
592		* @return
593		*/
		81.8% Uncovered Elements: 2 (11) Complexity: 3 Complexity Density: 0.43
594	24	protected int loadSequenceVCF(SequenceI seq, String vcfAssembly)...
595		{
596	24	VCFMap vcfMap = getVcfMap(seq, vcfAssembly);
597	24	if (vcfMap == null)
598		{
599	0	return 0;
600		}
601
602		/*
603		* work with the dataset sequence here
604		*/
605	24	SequenceI dss = seq.getDatasetSequence();
606	24	if (dss == null)
607		{
608	3	dss = seq;
609		}
610	24	return addVcfVariants(dss, vcfMap);
611		}
612
613		/**
614		* Answers a map from sequence coordinates to VCF chromosome ranges
615		*
616		* @param seq
617		* @param vcfAssembly
618		* @return
619		*/
		42.9% Uncovered Elements: 28 (49) Complexity: 10 Complexity Density: 0.3
620	24	private VCFMap getVcfMap(SequenceI seq, String vcfAssembly)...
621		{
622		/*
623		* simplest case: sequence has id and length matching a VCF contig
624		*/
625	24	VCFMap vcfMap = null;
626	24	if (dictionary != null)
627		{
628	3	vcfMap = getContigMap(seq);
629		}
630	24	if (vcfMap != null)
631		{
632	3	return vcfMap;
633		}
634
635		/*
636		* otherwise, map to VCF from chromosomal coordinates
637		* of the sequence (if known)
638		*/
639	21	GeneLociI seqCoords = seq.getGeneLoci();
640	21	if (seqCoords == null)
641		{
642	0	Cache.log.warn(String.format(
643		"Can't query VCF for %s as chromosome coordinates not known",
644		seq.getName()));
645	0	return null;
646		}
647
648	21	String species = seqCoords.getSpeciesId();
649	21	String chromosome = seqCoords.getChromosomeId();
650	21	String seqRef = seqCoords.getAssemblyId();
651	21	MapList map = seqCoords.getMap();
652
653	21	if (!vcfSpeciesMatchesSequence(vcfAssembly, species))
654		{
655	0	return null;
656		}
657
658	21	if (vcfAssemblyMatchesSequence(vcfAssembly, seqRef))
659		{
660	21	return new VCFMap(chromosome, map);
661		}
662
663	0	if (!"GRCh38".equalsIgnoreCase(seqRef) // Ensembl
664		|| !vcfAssembly.contains("Homo_sapiens_assembly19")) // gnomAD
665		{
666	0	return null;
667		}
668
669		/*
670		* map chromosomal coordinates from sequence to VCF if the VCF
671		* data has a different reference assembly to the sequence
672		*/
673		// TODO generalise for cases other than GRCh38 -> GRCh37 !
674		// - or get the user to choose in a dialog
675
676	0	List<int[]> toVcfRanges = new ArrayList<>();
677	0	List<int[]> fromSequenceRanges = new ArrayList<>();
678	0	String toRef = "GRCh37";
679
680	0	for (int[] range : map.getToRanges())
681		{
682	0	int[] fromRange = map.locateInFrom(range[0], range[1]);
683	0	if (fromRange == null)
684		{
685		// corrupted map?!?
686	0	continue;
687		}
688
689	0	int[] newRange = mapReferenceRange(range, chromosome, "human", seqRef,
690		toRef);
691	0	if (newRange == null)
692		{
693	0	Cache.log.error(
694		String.format("Failed to map %s:%s:%s:%d:%d to %s", species,
695		chromosome, seqRef, range[0], range[1], toRef));
696	0	continue;
697		}
698		else
699		{
700	0	toVcfRanges.add(newRange);
701	0	fromSequenceRanges.add(fromRange);
702		}
703		}
704
705	0	return new VCFMap(chromosome,
706		new MapList(fromSequenceRanges, toVcfRanges, 1, 1));
707		}
708
709		/**
710		* If the sequence id matches a contig declared in the VCF file, and the
711		* sequence length matches the contig length, then returns a 1:1 map of the
712		* sequence to the contig, else returns null
713		*
714		* @param seq
715		* @return
716		*/
		75% Uncovered Elements: 3 (12) Complexity: 3 Complexity Density: 0.38
717	3	private VCFMap getContigMap(SequenceI seq)...
718		{
719	3	String id = seq.getName();
720	3	SAMSequenceRecord contig = dictionary.getSequence(id);
721	3	if (contig != null)
722		{
723	3	int len = seq.getLength();
724	3	if (len == contig.getSequenceLength())
725		{
726	3	MapList map = new MapList(new int[] { 1, len },
727		new int[]
728		{ 1, len }, 1, 1);
729	3	return new VCFMap(id, map);
730		}
731		}
732	0	return null;
733		}
734
735		/**
736		* Answers true if we determine that the VCF data uses the same reference
737		* assembly as the sequence, else false
738		*
739		* @param vcfAssembly
740		* @param seqRef
741		* @return
742		*/
		60% Uncovered Elements: 2 (5) Complexity: 3 Complexity Density: 1
743	21	private boolean vcfAssemblyMatchesSequence(String vcfAssembly,...
744		String seqRef)
745		{
746		// TODO improve on this stub, which handles gnomAD and
747		// hopes for the best for other cases
748
749	21	if ("GRCh38".equalsIgnoreCase(seqRef) // Ensembl
750		&& vcfAssembly.contains("Homo_sapiens_assembly19")) // gnomAD
751		{
752	0	return false;
753		}
754	21	return true;
755		}
756
757		/**
758		* Answers true if the species inferred from the VCF reference identifier
759		* matches that for the sequence
760		*
761		* @param vcfAssembly
762		* @param speciesId
763		* @return
764		*/
		33.3% Uncovered Elements: 6 (9) Complexity: 5 Complexity Density: 1
765	21	boolean vcfSpeciesMatchesSequence(String vcfAssembly, String speciesId)...
766		{
767		// PROBLEM 1
768		// there are many aliases for species - how to equate one with another?
769		// PROBLEM 2
770		// VCF ##reference header is an unstructured URI - how to extract species?
771		// perhaps check if ref includes any (Ensembl) alias of speciesId??
772		// TODO ask the user to confirm this??
773
774	21	if (vcfAssembly.contains("Homo_sapiens") // gnomAD exome data example
775		&& "HOMO_SAPIENS".equals(speciesId)) // Ensembl species id
776		{
777	21	return true;
778		}
779
780	0	if (vcfAssembly.contains("c_elegans") // VEP VCF response example
781		&& "CAENORHABDITIS_ELEGANS".equals(speciesId)) // Ensembl
782		{
783	0	return true;
784		}
785
786		// this is not a sustainable solution...
787
788	0	return false;
789		}
790
791		/**
792		* Queries the VCF reader for any variants that overlap the mapped chromosome
793		* ranges of the sequence, and adds as variant features. Returns the number of
794		* overlapping variants found.
795		*
796		* @param seq
797		* @param map
798		* mapping from sequence to VCF coordinates
799		* @return
800		*/
		100% Uncovered Elements: 0 (19) Complexity: 3 Complexity Density: 0.2
801	24	protected int addVcfVariants(SequenceI seq, VCFMap map)...
802		{
803	24	boolean forwardStrand = map.map.isToForwardStrand();
804
805		/*
806		* query the VCF for overlaps of each contiguous chromosomal region
807		*/
808	24	int count = 0;
809
810	24	for (int[] range : map.map.getToRanges())
811		{
812	39	int vcfStart = Math.min(range[0], range[1]);
813	39	int vcfEnd = Math.max(range[0], range[1]);
814	39	CloseableIterator<VariantContext> variants = reader
815		.query(map.chromosome, vcfStart, vcfEnd);
816	67	while (variants.hasNext())
817		{
818	28	VariantContext variant = variants.next();
819
820	28	int[] featureRange = map.map.locateInFrom(variant.getStart(),
821		variant.getEnd());
822
823	28	if (featureRange != null)
824		{
825	26	int featureStart = Math.min(featureRange[0], featureRange[1]);
826	26	int featureEnd = Math.max(featureRange[0], featureRange[1]);
827	26	count += addAlleleFeatures(seq, variant, featureStart, featureEnd,
828		forwardStrand);
829		}
830		}
831	39	variants.close();
832		}
833
834	24	return count;
835		}
836
837		/**
838		* A convenience method to get the AF value for the given alternate allele
839		* index
840		*
841		* @param variant
842		* @param alleleIndex
843		* @return
844		*/
		87.5% Uncovered Elements: 1 (8) Complexity: 3 Complexity Density: 0.5
845	43	protected float getAlleleFrequency(VariantContext variant, int alleleIndex)...
846		{
847	43	float score = 0f;
848	43	String attributeValue = getAttributeValue(variant,
849		ALLELE_FREQUENCY_KEY, alleleIndex);
850	43	if (attributeValue != null)
851		{
852	43	try
853		{
854	43	score = Float.parseFloat(attributeValue);
855		} catch (NumberFormatException e)
856		{
857		// leave as 0
858		}
859		}
860
861	43	return score;
862		}
863
864		/**
865		* A convenience method to get an attribute value for an alternate allele
866		*
867		* @param variant
868		* @param attributeName
869		* @param alleleIndex
870		* @return
871		*/
		80% Uncovered Elements: 2 (10) Complexity: 3 Complexity Density: 0.5
872	119	protected String getAttributeValue(VariantContext variant,...
873		String attributeName, int alleleIndex)
874		{
875	119	Object att = variant.getAttribute(attributeName);
876
877	119	if (att instanceof String)
878		{
879	40	return (String) att;
880		}
881	79	else if (att instanceof ArrayList)
882		{
883	79	return ((List<String>) att).get(alleleIndex);
884		}
885
886	0	return null;
887		}
888
889		/**
890		* Adds one variant feature for each allele in the VCF variant record, and
891		* returns the number of features added.
892		*
893		* @param seq
894		* @param variant
895		* @param featureStart
896		* @param featureEnd
897		* @param forwardStrand
898		* @return
899		*/
		100% Uncovered Elements: 0 (7) Complexity: 2 Complexity Density: 0.4
900	26	protected int addAlleleFeatures(SequenceI seq, VariantContext variant,...
901		int featureStart, int featureEnd, boolean forwardStrand)
902		{
903	26	int added = 0;
904
905		/*
906		* Javadoc says getAlternateAlleles() imposes no order on the list returned
907		* so we proceed defensively to get them in strict order
908		*/
909	26	int altAlleleCount = variant.getAlternateAlleles().size();
910	69	for (int i = 0; i < altAlleleCount; i++)
911		{
912	43	added += addAlleleFeature(seq, variant, i, featureStart, featureEnd,
913		forwardStrand);
914		}
915	26	return added;
916		}
917
918		/**
919		* Inspects one allele and attempts to add a variant feature for it to the
920		* sequence. The additional data associated with this allele is extracted to
921		* store in the feature's key-value map. Answers the number of features added (0
922		* or 1).
923		*
924		* @param seq
925		* @param variant
926		* @param altAlleleIndex
927		* (0, 1..)
928		* @param featureStart
929		* @param featureEnd
930		* @param forwardStrand
931		* @return
932		*/
		100% Uncovered Elements: 0 (34) Complexity: 7 Complexity Density: 0.27
933	43	protected int addAlleleFeature(SequenceI seq, VariantContext variant,...
934		int altAlleleIndex, int featureStart, int featureEnd,
935		boolean forwardStrand)
936		{
937	43	String reference = variant.getReference().getBaseString();
938	43	Allele alt = variant.getAlternateAllele(altAlleleIndex);
939	43	String allele = alt.getBaseString();
940
941		/*
942		* insertion after a genomic base, if on reverse strand, has to be
943		* converted to insertion of complement after the preceding position
944		*/
945	43	int referenceLength = reference.length();
946	43	if (!forwardStrand && allele.length() > referenceLength
947		&& allele.startsWith(reference))
948		{
949	4	featureStart -= referenceLength;
950	4	featureEnd = featureStart;
951	4	char insertAfter = seq.getCharAt(featureStart - seq.getStart());
952	4	reference = Dna.reverseComplement(String.valueOf(insertAfter));
953	4	allele = allele.substring(referenceLength) + reference;
954		}
955
956		/*
957		* build the ref,alt allele description e.g. "G,A", using the base
958		* complement if the sequence is on the reverse strand
959		*/
960	43	StringBuilder sb = new StringBuilder();
961	43	sb.append(forwardStrand ? reference : Dna.reverseComplement(reference));
962	43	sb.append(COMMA);
963	43	sb.append(forwardStrand ? allele : Dna.reverseComplement(allele));
964	43	String alleles = sb.toString(); // e.g. G,A
965
966		/*
967		* pick out the consequence data (if any) that is for the current allele
968		* and feature (transcript) that matches the current sequence
969		*/
970	43	String consequence = getConsequenceForAlleleAndFeature(variant, CSQ_FIELD,
971		altAlleleIndex, csqAlleleFieldIndex,
972		csqAlleleNumberFieldIndex, seq.getName().toLowerCase(),
973		csqFeatureFieldIndex);
974
975		/*
976		* pick out the ontology term for the consequence type
977		*/
978	43	String type = SequenceOntologyI.SEQUENCE_VARIANT;
979	43	if (consequence != null)
980		{
981	7	type = getOntologyTerm(consequence);
982		}
983
984	43	float score = getAlleleFrequency(variant, altAlleleIndex);
985
986	43	SequenceFeature sf = new SequenceFeature(type, alleles, featureStart,
987		featureEnd, score, FEATURE_GROUP_VCF);
988	43	sf.setSource(sourceId);
989
990	43	sf.setValue(Gff3Helper.ALLELES, alleles);
991
992	43	addAlleleProperties(variant, sf, altAlleleIndex, consequence);
993
994	43	seq.addSequenceFeature(sf);
995
996	43	return 1;
997		}
998
999		/**
1000		* Determines the Sequence Ontology term to use for the variant feature type in
1001		* Jalview. The default is 'sequence_variant', but a more specific term is used
1002		* if:
1003		* <ul>
1004		* <li>VEP (or SnpEff) Consequence annotation is included in the VCF</li>
1005		* <li>sequence id can be matched to VEP Feature (or SnpEff Feature_ID)</li>
1006		* </ul>
1007		*
1008		* @param consequence
1009		* @return
1010		* @see http://www.sequenceontology.org/browser/current_svn/term/SO:0001060
1011		*/
		68.2% Uncovered Elements: 7 (22) Complexity: 6 Complexity Density: 0.5
1012	7	String getOntologyTerm(String consequence)...
1013		{
1014	7	String type = SequenceOntologyI.SEQUENCE_VARIANT;
1015
1016		/*
1017		* could we associate Consequence data with this allele and feature (transcript)?
1018		* if so, prefer the consequence term from that data
1019		*/
1020	7	if (csqAlleleFieldIndex == -1) // && snpEffAlleleFieldIndex == -1
1021		{
1022		/*
1023		* no Consequence data so we can't refine the ontology term
1024		*/
1025	0	return type;
1026		}
1027
1028	7	if (consequence != null)
1029		{
1030	7	String[] csqFields = consequence.split(PIPE_REGEX);
1031	7	if (csqFields.length > csqConsequenceFieldIndex)
1032		{
1033	7	type = csqFields[csqConsequenceFieldIndex];
1034		}
1035		}
1036		else
1037		{
1038		// todo the same for SnpEff consequence data matching if wanted
1039		}
1040
1041		/*
1042		* if of the form (e.g.) missense_variant&splice_region_variant,
1043		* just take the first ('most severe') consequence
1044		*/
1045	7	if (type != null)
1046		{
1047	7	int pos = type.indexOf('&');
1048	7	if (pos > 0)
1049		{
1050	0	type = type.substring(0, pos);
1051		}
1052		}
1053	7	return type;
1054		}
1055
1056		/**
1057		* Returns matched consequence data if it can be found, else null.
1058		* <ul>
1059		* <li>inspects the VCF data for key 'vcfInfoId'</li>
1060		* <li>splits this on comma (to distinct consequences)</li>
1061		* <li>returns the first consequence (if any) where</li>
1062		* <ul>
1063		* <li>the allele matches the altAlleleIndex'th allele of variant</li>
1064		* <li>the feature matches the sequence name (e.g. transcript id)</li>
1065		* </ul>
1066		* </ul>
1067		* If matched, the consequence is returned (as pipe-delimited fields).
1068		*
1069		* @param variant
1070		* @param vcfInfoId
1071		* @param altAlleleIndex
1072		* @param alleleFieldIndex
1073		* @param alleleNumberFieldIndex
1074		* @param seqName
1075		* @param featureFieldIndex
1076		* @return
1077		*/
		87.5% Uncovered Elements: 3 (24) Complexity: 9 Complexity Density: 0.64
1078	43	private String getConsequenceForAlleleAndFeature(VariantContext variant,...
1079		String vcfInfoId, int altAlleleIndex, int alleleFieldIndex,
1080		int alleleNumberFieldIndex,
1081		String seqName, int featureFieldIndex)
1082		{
1083	43	if (alleleFieldIndex == -1 || featureFieldIndex == -1)
1084		{
1085	29	return null;
1086		}
1087	14	Object value = variant.getAttribute(vcfInfoId);
1088
1089	14	if (value == null || !(value instanceof List<?>))
1090		{
1091	0	return null;
1092		}
1093
1094		/*
1095		* inspect each consequence in turn (comma-separated blocks
1096		* extracted by htsjdk)
1097		*/
1098	14	List<String> consequences = (List<String>) value;
1099
1100	14	for (String consequence : consequences)
1101		{
1102	41	String[] csqFields = consequence.split(PIPE_REGEX);
1103	41	if (csqFields.length > featureFieldIndex)
1104		{
1105	41	String featureIdentifier = csqFields[featureFieldIndex];
1106	41	if (featureIdentifier.length() > 4
1107		&& seqName.indexOf(featureIdentifier.toLowerCase()) > -1)
1108		{
1109		/*
1110		* feature (transcript) matched - now check for allele match
1111		*/
1112	10	if (matchAllele(variant, altAlleleIndex, csqFields,
1113		alleleFieldIndex, alleleNumberFieldIndex))
1114		{
1115	7	return consequence;
1116		}
1117		}
1118		}
1119		}
1120	7	return null;
1121		}
1122
		43.8% Uncovered Elements: 9 (16) Complexity: 5 Complexity Density: 0.5
1123	10	private boolean matchAllele(VariantContext variant, int altAlleleIndex,...
1124		String[] csqFields, int alleleFieldIndex,
1125		int alleleNumberFieldIndex)
1126		{
1127		/*
1128		* if ALLELE_NUM is present, it must match altAlleleIndex
1129		* NB first alternate allele is 1 for ALLELE_NUM, 0 for altAlleleIndex
1130		*/
1131	10	if (alleleNumberFieldIndex > -1)
1132		{
1133	0	if (csqFields.length <= alleleNumberFieldIndex)
1134		{
1135	0	return false;
1136		}
1137	0	String alleleNum = csqFields[alleleNumberFieldIndex];
1138	0	return String.valueOf(altAlleleIndex + 1).equals(alleleNum);
1139		}
1140
1141		/*
1142		* else consequence allele must match variant allele
1143		*/
1144	10	if (alleleFieldIndex > -1 && csqFields.length > alleleFieldIndex)
1145		{
1146	10	String csqAllele = csqFields[alleleFieldIndex];
1147	10	String vcfAllele = variant.getAlternateAllele(altAlleleIndex)
1148		.getBaseString();
1149	10	return csqAllele.equals(vcfAllele);
1150		}
1151	0	return false;
1152		}
1153
1154		/**
1155		* Add any allele-specific VCF key-value data to the sequence feature
1156		*
1157		* @param variant
1158		* @param sf
1159		* @param altAlelleIndex
1160		* (0, 1..)
1161		* @param consequence
1162		* if not null, the consequence specific to this sequence (transcript
1163		* feature) and allele
1164		*/
		86.1% Uncovered Elements: 5 (36) Complexity: 8 Complexity Density: 0.36
1165	43	protected void addAlleleProperties(VariantContext variant,...
1166		SequenceFeature sf, final int altAlelleIndex, String consequence)
1167		{
1168	43	Map<String, Object> atts = variant.getAttributes();
1169
1170	43	for (Entry<String, Object> att : atts.entrySet())
1171		{
1172	142	String key = att.getKey();
1173
1174		/*
1175		* extract Consequence data (if present) that we are able to
1176		* associated with the allele for this variant feature
1177		*/
1178	142	if (CSQ_FIELD.equals(key))
1179		{
1180	14	addConsequences(variant, sf, consequence);
1181	14	continue;
1182		}
1183
1184		/*
1185		* filter out fields we don't want to capture
1186		*/
1187	128	if (!vcfFieldsOfInterest.contains(key))
1188		{
1189	38	continue;
1190		}
1191
1192		/*
1193		* filter out fields we don't want to capture
1194		*/
1195	90	if (!vcfFieldsOfInterest.contains(key))
1196		{
1197	0	continue;
1198		}
1199
1200		/*
1201		* we extract values for other data which are allele-specific;
1202		* these may be per alternate allele (INFO[key].Number = 'A')
1203		* or per allele including reference (INFO[key].Number = 'R')
1204		*/
1205	90	VCFInfoHeaderLine infoHeader = header.getInfoHeaderLine(key);
1206	90	if (infoHeader == null)
1207		{
1208		/*
1209		* can't be sure what data belongs to this allele, so
1210		* play safe and don't take any
1211		*/
1212	0	continue;
1213		}
1214
1215	90	VCFHeaderLineCount number = infoHeader.getCountType();
1216	90	int index = altAlelleIndex;
1217	90	if (number == VCFHeaderLineCount.R)
1218		{
1219		/*
1220		* one value per allele including reference, so bump index
1221		* e.g. the 3rd value is for the 2nd alternate allele
1222		*/
1223	14	index++;
1224		}
1225	76	else if (number != VCFHeaderLineCount.A)
1226		{
1227		/*
1228		* don't save other values as not allele-related
1229		*/
1230	14	continue;
1231		}
1232
1233		/*
1234		* take the index'th value
1235		*/
1236	76	String value = getAttributeValue(variant, key, index);
1237	76	if (value != null)
1238		{
1239	76	sf.setValue(key, value);
1240		}
1241		}
1242		}
1243
1244		/**
1245		* Inspects CSQ data blocks (consequences) and adds attributes on the sequence
1246		* feature.
1247		* <p>
1248		* If <code>myConsequence</code> is not null, then this is the specific
1249		* consequence data (pipe-delimited fields) that is for the current allele and
1250		* transcript (sequence) being processed)
1251		*
1252		* @param variant
1253		* @param sf
1254		* @param myConsequence
1255		*/
		81.5% Uncovered Elements: 5 (27) Complexity: 9 Complexity Density: 0.53
1256	14	protected void addConsequences(VariantContext variant, SequenceFeature sf,...
1257		String myConsequence)
1258		{
1259	14	Object value = variant.getAttribute(CSQ_FIELD);
1260
1261	14	if (value == null || !(value instanceof List<?>))
1262		{
1263	0	return;
1264		}
1265
1266	14	List<String> consequences = (List<String>) value;
1267
1268		/*
1269		* inspect CSQ consequences; restrict to the consequence
1270		* associated with the current transcript (Feature)
1271		*/
1272	14	Map<String, String> csqValues = new HashMap<>();
1273
1274	14	for (String consequence : consequences)
1275		{
1276	50	if (myConsequence == null || myConsequence.equals(consequence))
1277		{
1278	31	String[] csqFields = consequence.split(PIPE_REGEX);
1279
1280		/*
1281		* inspect individual fields of this consequence, copying non-null
1282		* values which are 'fields of interest'
1283		*/
1284	31	int i = 0;
1285	31	for (String field : csqFields)
1286		{
1287	279	if (field != null && field.length() > 0)
1288		{
1289	279	String id = vepFieldsOfInterest.get(i);
1290	279	if (id != null)
1291		{
1292	279	csqValues.put(id, field);
1293		}
1294		}
1295	279	i++;
1296		}
1297		}
1298		}
1299
1300	14	if (!csqValues.isEmpty())
1301		{
1302	14	sf.setValue(CSQ_FIELD, csqValues);
1303		}
1304		}
1305
1306		/**
1307		* A convenience method to complement a dna base and return the string value
1308		* of its complement
1309		*
1310		* @param reference
1311		* @return
1312		*/
		0% Uncovered Elements: 1 (1) Complexity: 1 Complexity Density: 1
1313	0	protected String complement(byte[] reference)...
1314		{
1315	0	return String.valueOf(Dna.getComplement((char) reference[0]));
1316		}
1317
1318		/**
1319		* Determines the location of the query range (chromosome positions) in a
1320		* different reference assembly.
1321		* <p>
1322		* If the range is just a subregion of one for which we already have a mapping
1323		* (for example, an exon sub-region of a gene), then the mapping is just
1324		* computed arithmetically.
1325		* <p>
1326		* Otherwise, calls the Ensembl REST service that maps from one assembly
1327		* reference's coordinates to another's
1328		*
1329		* @param queryRange
1330		* start-end chromosomal range in 'fromRef' coordinates
1331		* @param chromosome
1332		* @param species
1333		* @param fromRef
1334		* assembly reference for the query coordinates
1335		* @param toRef
1336		* assembly reference we wish to translate to
1337		* @return the start-end range in 'toRef' coordinates
1338		*/
		0% Uncovered Elements: 18 (18) Complexity: 4 Complexity Density: 0.33
1339	0	protected int[] mapReferenceRange(int[] queryRange, String chromosome,...
1340		String species, String fromRef, String toRef)
1341		{
1342		/*
1343		* first try shorcut of computing the mapping as a subregion of one
1344		* we already have (e.g. for an exon, if we have the gene mapping)
1345		*/
1346	0	int[] mappedRange = findSubsumedRangeMapping(queryRange, chromosome,
1347		species, fromRef, toRef);
1348	0	if (mappedRange != null)
1349		{
1350	0	return mappedRange;
1351		}
1352
1353		/*
1354		* call (e.g.) http://rest.ensembl.org/map/human/GRCh38/17:45051610..45109016:1/GRCh37
1355		*/
1356	0	EnsemblMap mapper = new EnsemblMap();
1357	0	int[] mapping = mapper.getAssemblyMapping(species, chromosome, fromRef,
1358		toRef, queryRange);
1359
1360	0	if (mapping == null)
1361		{
1362		// mapping service failure
1363	0	return null;
1364		}
1365
1366		/*
1367		* save mapping for possible future re-use
1368		*/
1369	0	String key = makeRangesKey(chromosome, species, fromRef, toRef);
1370	0	if (!assemblyMappings.containsKey(key))
1371		{
1372	0	assemblyMappings.put(key, new HashMap<int[], int[]>());
1373		}
1374
1375	0	assemblyMappings.get(key).put(queryRange, mapping);
1376
1377	0	return mapping;
1378		}
1379
1380		/**
1381		* If we already have a 1:1 contiguous mapping which subsumes the given query
1382		* range, this method just calculates and returns the subset of that mapping,
1383		* else it returns null. In practical terms, if a gene has a contiguous
1384		* mapping between (for example) GRCh37 and GRCh38, then we assume that its
1385		* subsidiary exons occupy unchanged relative positions, and just compute
1386		* these as offsets, rather than do another lookup of the mapping.
1387		* <p>
1388		* If in future these assumptions prove invalid (e.g. for bacterial dna?!),
1389		* simply remove this method or let it always return null.
1390		* <p>
1391		* Warning: many rapid calls to the /map service map result in a 429 overload
1392		* error response
1393		*
1394		* @param queryRange
1395		* @param chromosome
1396		* @param species
1397		* @param fromRef
1398		* @param toRef
1399		* @return
1400		*/
		0% Uncovered Elements: 19 (19) Complexity: 4 Complexity Density: 0.31
1401	0	protected int[] findSubsumedRangeMapping(int[] queryRange, String chromosome,...
1402		String species, String fromRef, String toRef)
1403		{
1404	0	String key = makeRangesKey(chromosome, species, fromRef, toRef);
1405	0	if (assemblyMappings.containsKey(key))
1406		{
1407	0	Map<int[], int[]> mappedRanges = assemblyMappings.get(key);
1408	0	for (Entry<int[], int[]> mappedRange : mappedRanges.entrySet())
1409		{
1410	0	int[] fromRange = mappedRange.getKey();
1411	0	int[] toRange = mappedRange.getValue();
1412	0	if (fromRange[1] - fromRange[0] == toRange[1] - toRange[0])
1413		{
1414		/*
1415		* mapping is 1:1 in length, so we trust it to have no discontinuities
1416		*/
1417	0	if (MappingUtils.rangeContains(fromRange, queryRange))
1418		{
1419		/*
1420		* fromRange subsumes our query range
1421		*/
1422	0	int offset = queryRange[0] - fromRange[0];
1423	0	int mappedRangeFrom = toRange[0] + offset;
1424	0	int mappedRangeTo = mappedRangeFrom + (queryRange[1] - queryRange[0]);
1425	0	return new int[] { mappedRangeFrom, mappedRangeTo };
1426		}
1427		}
1428		}
1429		}
1430	0	return null;
1431		}
1432
1433		/**
1434		* Transfers the sequence feature to the target sequence, locating its start
1435		* and end range based on the mapping. Features which do not overlap the
1436		* target sequence are ignored.
1437		*
1438		* @param sf
1439		* @param targetSequence
1440		* @param mapping
1441		* mapping from the feature's coordinates to the target sequence
1442		*/
		0% Uncovered Elements: 9 (9) Complexity: 2 Complexity Density: 0.29
1443	0	protected void transferFeature(SequenceFeature sf,...
1444		SequenceI targetSequence, MapList mapping)
1445		{
1446	0	int[] mappedRange = mapping.locateInTo(sf.getBegin(), sf.getEnd());
1447
1448	0	if (mappedRange != null)
1449		{
1450	0	String group = sf.getFeatureGroup();
1451	0	int newBegin = Math.min(mappedRange[0], mappedRange[1]);
1452	0	int newEnd = Math.max(mappedRange[0], mappedRange[1]);
1453	0	SequenceFeature copy = new SequenceFeature(sf, newBegin, newEnd,
1454		group, sf.getScore());
1455	0	targetSequence.addSequenceFeature(copy);
1456		}
1457		}
1458
1459		/**
1460		* Formats a ranges map lookup key
1461		*
1462		* @param chromosome
1463		* @param species
1464		* @param fromRef
1465		* @param toRef
1466		* @return
1467		*/
		0% Uncovered Elements: 1 (1) Complexity: 1 Complexity Density: 1
1468	0	protected static String makeRangesKey(String chromosome, String species,...
1469		String fromRef, String toRef)
1470		{
1471	0	return species + EXCL + chromosome + EXCL + fromRef + EXCL
1472		+ toRef;
1473		}
1474		}