X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fio%2Fvcf%2FVCFLoader.java;h=ea5b8e06f723f3ac512a903c780508e2158184f5;hb=57738a1f3c19b1c3a00bd3ac5108f8cd0af32f99;hp=168f1c6499bfcaf4a1d9bd437bb5afe4dc9eeda5;hpb=006890b02106eb31841e6e84d75f1027434823e0;p=jalview.git diff --git a/src/jalview/io/vcf/VCFLoader.java b/src/jalview/io/vcf/VCFLoader.java index 168f1c6..ea5b8e0 100644 --- a/src/jalview/io/vcf/VCFLoader.java +++ b/src/jalview/io/vcf/VCFLoader.java @@ -1,25 +1,26 @@ +/* + * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$) + * Copyright (C) $$Year-Rel$$ The Jalview Authors + * + * This file is part of Jalview. + * + * Jalview is free software: you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation, either version 3 + * of the License, or (at your option) any later version. + * + * Jalview is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR + * PURPOSE. See the GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with Jalview. If not, see . + * The Jalview Authors are detailed in the 'AUTHORS' file. + */ package jalview.io.vcf; -import jalview.analysis.Dna; -import jalview.api.AlignViewControllerGuiI; -import jalview.bin.Cache; -import jalview.datamodel.DBRefEntry; -import jalview.datamodel.GeneLociI; -import jalview.datamodel.Mapping; -import jalview.datamodel.SequenceFeature; -import jalview.datamodel.SequenceI; -import jalview.datamodel.features.FeatureAttributeType; -import jalview.datamodel.features.FeatureSource; -import jalview.datamodel.features.FeatureSources; -import jalview.ext.ensembl.EnsemblMap; -import jalview.ext.htsjdk.HtsContigDb; -import jalview.ext.htsjdk.VCFReader; -import jalview.io.gff.Gff3Helper; -import jalview.io.gff.SequenceOntologyI; -import jalview.util.MapList; -import jalview.util.MappingUtils; -import jalview.util.MessageManager; -import jalview.util.StringUtils; +import java.util.Locale; import java.io.File; import java.io.IOException; @@ -47,6 +48,27 @@ import htsjdk.variant.vcf.VCFHeaderLine; import htsjdk.variant.vcf.VCFHeaderLineCount; import htsjdk.variant.vcf.VCFHeaderLineType; import htsjdk.variant.vcf.VCFInfoHeaderLine; +import jalview.analysis.Dna; +import jalview.api.AlignViewControllerGuiI; +import jalview.bin.Cache; +import jalview.bin.Console; +import jalview.datamodel.DBRefEntry; +import jalview.datamodel.GeneLociI; +import jalview.datamodel.Mapping; +import jalview.datamodel.SequenceFeature; +import jalview.datamodel.SequenceI; +import jalview.datamodel.features.FeatureAttributeType; +import jalview.datamodel.features.FeatureSource; +import jalview.datamodel.features.FeatureSources; +import jalview.ext.ensembl.EnsemblMap; +import jalview.ext.htsjdk.HtsContigDb; +import jalview.ext.htsjdk.VCFReader; +import jalview.io.gff.Gff3Helper; +import jalview.io.gff.SequenceOntologyI; +import jalview.util.MapList; +import jalview.util.MappingUtils; +import jalview.util.MessageManager; +import jalview.util.StringUtils; /** * A class to read VCF data (using the htsjdk) and add variants as sequence @@ -74,13 +96,14 @@ public class VCFLoader private static final String DEFAULT_SPECIES = "homo_sapiens"; /** - * A class to model the mapping from sequence to VCF coordinates. Cases include + * A class to model the mapping from sequence to VCF coordinates. Cases + * include * */ class VCFMap @@ -124,15 +147,20 @@ public class VCFLoader private static final String VCF_SPECIES = "VCF_SPECIES"; // default is human - private static final String DEFAULT_REFERENCE = "grch37"; // fallback default is human GRCh37 + private static final String DEFAULT_REFERENCE = "grch37"; // fallback default + // is human GRCh37 /* * keys to fields of VEP CSQ consequence data * see https://www.ensembl.org/info/docs/tools/vep/vep_formats.html */ private static final String CSQ_CONSEQUENCE_KEY = "Consequence"; + private static final String CSQ_ALLELE_KEY = "Allele"; - private static final String CSQ_ALLELE_NUM_KEY = "ALLELE_NUM"; // 0 (ref), 1... + + private static final String CSQ_ALLELE_NUM_KEY = "ALLELE_NUM"; // 0 (ref), + // 1... + private static final String CSQ_FEATURE_KEY = "Feature"; // Ensembl stable id /* @@ -203,8 +231,11 @@ public class VCFLoader * see http://www.ensembl.org/info/docs/tools/vep/vep_formats.html */ private int csqConsequenceFieldIndex = -1; + private int csqAlleleFieldIndex = -1; + private int csqAlleleNumberFieldIndex = -1; + private int csqFeatureFieldIndex = -1; // todo the same fields for SnpEff ANN data if wanted @@ -290,10 +321,11 @@ public class VCFLoader */ public SequenceI loadVCFContig(String contig) { - VCFHeaderLine headerLine = header.getOtherHeaderLine(VCFHeader.REFERENCE_KEY); + VCFHeaderLine headerLine = header + .getOtherHeaderLine(VCFHeader.REFERENCE_KEY); if (headerLine == null) { - Cache.log.error("VCF reference header not found"); + Console.error("VCF reference header not found"); return null; } String ref = headerLine.getValue(); @@ -315,7 +347,7 @@ public class VCFLoader } else { - Cache.log.error("VCF reference not found: " + ref); + Console.error("VCF reference not found: " + ref); } return seq; @@ -391,9 +423,9 @@ public class VCFLoader /** * Attempts to determine and save the species and genome assembly version to - * which the VCF data applies. This may be done by parsing the {@code reference} - * header line, configured in a property file, or (potentially) confirmed - * interactively by the user. + * which the VCF data applies. This may be done by parsing the + * {@code reference} header line, configured in a property file, or + * (potentially) confirmed interactively by the user. *

* The saved values should be identifiers valid for Ensembl's REST service * {@code map} endpoint, so they can be used (if necessary) to retrieve the @@ -408,11 +440,11 @@ public class VCFLoader { if (reference == null) { - Cache.log.error("No VCF ##reference found, defaulting to " + Console.error("No VCF ##reference found, defaulting to " + DEFAULT_REFERENCE + ":" + DEFAULT_SPECIES); reference = DEFAULT_REFERENCE; // default to GRCh37 if not specified } - reference = reference.toLowerCase(); + reference = reference.toLowerCase(Locale.ROOT); /* * for a non-human species, or other assembly identifier, @@ -427,7 +459,7 @@ public class VCFLoader String[] tokens = token.split("="); if (tokens.length == 2) { - if (reference.contains(tokens[0].trim().toLowerCase())) + if (reference.contains(tokens[0].trim().toLowerCase(Locale.ROOT))) { vcfAssembly = tokens[1].trim(); break; @@ -444,7 +476,7 @@ public class VCFLoader String[] tokens = token.split("="); if (tokens.length == 2) { - if (reference.contains(tokens[0].trim().toLowerCase())) + if (reference.contains(tokens[0].trim().toLowerCase(Locale.ROOT))) { vcfSpecies = tokens[1].trim(); break; @@ -549,7 +581,7 @@ public class VCFLoader { for (Pattern p : filters) { - if (p.matcher(id.toUpperCase()).matches()) + if (p.matcher(id.toUpperCase(Locale.ROOT)).matches()) { return true; } @@ -643,7 +675,7 @@ public class VCFLoader { try { - patterns.add(Pattern.compile(token.toUpperCase())); + patterns.add(Pattern.compile(token.toUpperCase(Locale.ROOT))); } catch (PatternSyntaxException e) { System.err.println("Invalid pattern ignored: " + token); @@ -654,7 +686,6 @@ public class VCFLoader /** * Transfers VCF features to sequences to which this sequence has a mapping. - * If the mapping is 3:1, computes peptide variants from nucleotide variants. * * @param seq */ @@ -702,8 +733,8 @@ public class VCFLoader } /** - * Tries to add overlapping variants read from a VCF file to the given sequence, - * and returns the number of variant features added + * Tries to add overlapping variants read from a VCF file to the given + * sequence, and returns the number of variant features added * * @param seq * @return @@ -755,7 +786,7 @@ public class VCFLoader GeneLociI seqCoords = seq.getGeneLoci(); if (seqCoords == null) { - Cache.log.warn(String.format( + Console.warn(String.format( "Can't query VCF for %s as chromosome coordinates not known", seq.getName())); return null; @@ -770,7 +801,7 @@ public class VCFLoader // returned with the Ensembl sequence; todo: support aliases? if (!vcfSpecies.equalsIgnoreCase(species)) { - Cache.log.warn("No VCF loaded to " + seq.getName() + Console.warn("No VCF loaded to " + seq.getName() + " as species not matched"); return null; } @@ -800,10 +831,9 @@ public class VCFLoader vcfAssembly); if (newRange == null) { - Cache.log.error( - String.format("Failed to map %s:%s:%s:%d:%d to %s", species, - chromosome, seqRef, range[0], range[1], - vcfAssembly)); + Console.error(String.format("Failed to map %s:%s:%s:%d:%d to %s", + species, chromosome, seqRef, range[0], range[1], + vcfAssembly)); continue; } else @@ -877,12 +907,19 @@ public class VCFLoader int[] featureRange = map.map.locateInFrom(variant.getStart(), variant.getEnd()); + /* + * only take features whose range is fully mappable to sequence positions + */ if (featureRange != null) { int featureStart = Math.min(featureRange[0], featureRange[1]); int featureEnd = Math.max(featureRange[0], featureRange[1]); - count += addAlleleFeatures(seq, variant, featureStart, - featureEnd, forwardStrand); + if (featureEnd - featureStart == variant.getEnd() + - variant.getStart()) + { + count += addAlleleFeatures(seq, variant, featureStart, + featureEnd, forwardStrand); + } } } variants.close(); @@ -892,8 +929,8 @@ public class VCFLoader * RuntimeException throwable by htsjdk */ String msg = String.format("Error reading VCF for %s:%d-%d: %s ", - map.chromosome, vcfStart, vcfEnd); - Cache.log.error(msg); + map.chromosome, vcfStart, vcfEnd, e.getLocalizedMessage()); + Console.error(msg); } } @@ -957,8 +994,8 @@ public class VCFLoader /** * Inspects one allele and attempts to add a variant feature for it to the * sequence. The additional data associated with this allele is extracted to - * store in the feature's key-value map. Answers the number of features added (0 - * or 1). + * store in the feature's key-value map. Answers the number of features added + * (0 or 1). * * @param seq * @param variant @@ -1006,10 +1043,10 @@ public class VCFLoader * pick out the consequence data (if any) that is for the current allele * and feature (transcript) that matches the current sequence */ - String consequence = getConsequenceForAlleleAndFeature(variant, CSQ_FIELD, - altAlleleIndex, csqAlleleFieldIndex, - csqAlleleNumberFieldIndex, seq.getName().toLowerCase(), - csqFeatureFieldIndex); + String consequence = getConsequenceForAlleleAndFeature(variant, + CSQ_FIELD, altAlleleIndex, csqAlleleFieldIndex, + csqAlleleNumberFieldIndex, + seq.getName().toLowerCase(Locale.ROOT), csqFeatureFieldIndex); /* * pick out the ontology term for the consequence type @@ -1049,8 +1086,8 @@ public class VCFLoader /** * Answers the VCF FILTER value for the variant - or an approximation to it. * This field is either PASS, or a semi-colon separated list of filters not - * passed. htsjdk saves filters as a HashSet, so the order when reassembled into - * a list may be different. + * passed. htsjdk saves filters as a HashSet, so the order when reassembled + * into a list may be different. * * @param variant * @return @@ -1094,9 +1131,9 @@ public class VCFLoader } /** - * Determines the Sequence Ontology term to use for the variant feature type in - * Jalview. The default is 'sequence_variant', but a more specific term is used - * if: + * Determines the Sequence Ontology term to use for the variant feature type + * in Jalview. The default is 'sequence_variant', but a more specific term is + * used if: *