svart - (structural) variant tool
Svart is a small library for representing genomic variants and regions. It attempts to solve several common issues:
- Coordinate system off-by-one errors
- Representation of VCF small, structural and breakend variants with a consistent API
- Different variant trimming strategies
The library provides a consistent API for creating, manipulating and comparing variation of different types by providing default implementations and extensible base classes and interfaces along with immutable default implementations which developers can utilise to gain maximum utility from a minimum amount of code without having to address issues in bioinformatics which are a common source of duplicated code and frequently errors.
The code is completely free of external dependencies.
Motivation
This is intended to be used as a standard library for various Monarch and Monarch-related projects involving genomic variation such as Jannovar, Exomiser, Squirls, PBGA.
These projects are inter-dependent to some extent and all require some level of ability to represent and manipulate genomic variation.
Given the inter-dependency there is substantial copying/re-implementation in the code-bases of these projects of various core concepts surrounding genomic variation. This library aims to provide a common set of interfaces and implementations which can be used to fulfill their collective use-cases and so reduce code duplication, bugs, and GC pressure due to object conversion between common types.
Design
There is a (mostly) compositional hierarchy of data types:
Contig
Contig -has_length-> int
Contig -has_identifier-> int
Contig -has_name-> String
GenomicAssembly
GenomicAssembly -has_some-> Contig
CoordinateSystem
ONE_BASED or ZERO_BASED
Coordinates
Coordinates -has_a-> CoordinateSystem
Coordinates -has_start-> int
Coordinates -has_end-> int
Region
Region -has_some-> Coordinates
Strand
POSITIVE or NEGATIVE
GenomicRegion
GenomicRegion -is_a-> Region
GenomicRegion -has_a-> Contig
GenomicRegion -has_a-> Strand
GenomicVariant
GenomicVariant -is_a-> GenomicRegion
GenomicVariant -has_a-> VariantType
GenomicBreakendVariant
GenomicBreakendVariant -is_a-> GenomicVariant
GenomicBreakendVariant -has_left-> GenomicBreakend
GenomicBreakendVariant -has_right-> GenomicBreakend
GenomicBreakend -is_a-> GenomicRegion
With this model, it is possible to express variation, using zero or one-based coordinates on a strand of a chromosome, for a
given genome assembly in any coordinate system and easily manipulate and compare them without having to worry about what
strand or which coordinate system another GenomicRegion uses. Svart will do all this for you eliminating off-by-one
errors and providing an extensible model to plug into your code, so you can focus on the task at hand. It will enforce
Example
Modelling a gene as a region on a chromosome. Here we're using the gene FGFR2.
Chromosome 10 (CM000672.1): 123,237,848-123,353,481 reverse strand. (1-based, positive strand coordinates for the FGFR2 gene on GRCh37)
Chromosome 10 (CM000672.2): 121,478,332-121,598,458 reverse strand. (1-based, positive strand coordinates for the FGFR2 gene on GRCh38)The code below shows how to create a GenomicAssembly, create a GenomicRegion to represent the region of chromosome
10 to which the longest transcript expressing the FGFR2 can be aligned. We then create a Variant object representing a
SNV and check to see if this is located within the gene. FGFR2 is located on the negative strand, whereas variants from
VCF are reported on the positive strand, so at first sight it might not appear that the variant is contained within the
gene.
class SnpAndGeneTest {
public void checkGeneContainsVariant() {
// Load the Human GRCh37.13 assembly from a NCBI assembly report
GenomicAssembly b37 = GenomicAssembly.readAssembly(Path.of("src/test/resources/GCF_000001405.25_GRCh37.p13_assembly_report.txt"));
// alternatively you can use svart to download one using the GenBank or RefSeq accession:
// GenomicAssembly grch37p13 = GenomicAssemblies.downloadAssembly("GCF_000001405.25");
// or use a default assembly:
// GenomicAssembly grch37p13 = GenomicAssemblies.GRCh37p13();
// FGFR2 gene is located on chromosome 10 (CM000672.1): 123,237,848-123_357_972 reverse strand. (1-based, positive strand coordinates)
Contig chr10b37 = b37.contigByName("10");
GenomicRegion fgfr2Gene = GenomicRegion.of(chr10b37, Strand.POSITIVE, CoordinateSystem.oneBased(), 123_237_848, 123_357_972);
// 10 123256215 . T G - a pathogenic missense variant (GRCh37 VCF coordinates - 1-based, positive strand)
GenomicVariant snv = GenomicVariant.of(chr10b37, "", Strand.POSITIVE, CoordinateSystem.oneBased(), 123_256_215, "T", "G");
// Because svart knows about coordinate systems and strands it is possible to...
// keep the gene on the positive strand:
// GenomicRegion{contig=10, strand=+, coordinateSystem=ONE_BASED, startPosition=123237848, endPosition=123357972}
assertThat(fgfr2Gene.contains(snv), equalTo(true));
// or use it on the negative strand:
// GenomicRegion{contig=10, strand=-, coordinateSystem=ONE_BASED, startPosition=12176776, endPosition=12296900}
assertThat(fgfr2Gene.toNegativeStrand().contains(snv), equalTo(true));
}
}In the above example the default implementations provided by the library were used. These are immutable classes which
can be used compositionally. Alternatively, for users requiring extra functionality on top of the default implementations
there are several Base classes which can be extended - BaseGenomicRegion, BaseVariant and BaseBreakendVariant.
Representing variants - VCF
Svart models variants in a similar manner to VCF, however it attempts to provide a unified and consistent model using
the Variant interface. In VCF precise variants of known sequence (sequence variants) described using CHR, ID, POS,
REF, ALT and also the symbolic variants where the ALT field denotes the type in angle brackets e.g. <INS> for an
insertion along with a change length and end (the SVLEN and END fields in the INFO column).
Internally Svart treats both the sequence and symbolic types in the same way modelling them as GenomicRegion so
that all variants have a start and end (like BED), but also have a VariantType a refLength and a changeLength derived
from the REF/ALT sequences and END if available. The Variant interface contains static factory methods to help users
easily create variants, requiring them to consider important, but often implicit values for the Strand and CoordinateSystem
with which a variant is described.
Sequence and symbolic variants are always intra-chromosomal changes (i.e. on the same chromosome/contig) whereas
breakend variants are somewhat anomalous as they describe the way one or two different chromosomes are broken and
re-arranged. These are handled by the BreakendVariant which is composed of a left and right Breakend derived from a
single VCF record (line). The full re-arrangement requires assembling the Breakends together by matching their id
with another Breakend.mateId. The VCF representation of breakends is quite esoteric and relies on a pattern in the
ALT column of a record to indicate the strand, chromosome, location and orientation of a mate end relative to the REF
position of the record e.g. C[ctg2:5[ indicates ctg2 position 5 (fully-closed coordinates) on the positive strand
is to the right of the reference allele C also on the positive strand. Svart provides a VcfBreakendResolver to handle the parsing and conversion of
these into BreakendVariant instances.
While svart is heavily influenced by VCF it does not provide a parser for VCF as this is better handled by the HTSJDK,
however much of the terminology and data is very close, for example the ConfidenceInterval class is used in conjunction
with the Position class to encapsulate the POS/CIPOS and END/CIEND fields from VCF. Svart does, however provide a
VcfConverter class which provides methods to easily convert VCF small, symbolic and breakend records into svart Variant
instances from, for example an HTSJDK VariantContext. The VcfConverter needs to be provided with the correct
GenomicAssembly for the VCF file and a VariantTrimmer to allow the user to specify how they wish any multi-allelic
sites to be trimmed. For example, this is how the VCF fields for a SNV called against GRCh37 would be converted:
class ConvertVcfTest {
@Test
public void convert() {
VcfConverter vcfConverter = new VcfConverter(b37, VariantTrimmer.leftShiftingTrimmer(VariantTrimmer.retainingCommonBase()));
// non-symbolic alleles
// CHR POS ID REF ALT
// chr1 12345 rs123456 C T 6 PASS .
Variant snv = vcfConverter.convert(vcfConverter.parseContig("chr1"), "rs123456", 12345, "C", "T");
assertThat(snv, equalTo(Variant.of(chr1, "rs123456", Strand.POSITIVE, CoordinateSystem.ONE_BASED, 12345, "C", "T")));
// symbolic alleles
// chr1 12345 . C <INS> 6 PASS SVTYPE=INS;END=12345;SVLEN=200
Variant ins = vcfConverter.convertSymbolic(vcfConverter.parseContig("chr1"), "", 12345, 12345, "C", "<INS>", 200);
assertThat(ins, equalTo(Variant.of(chr1, "rs123456", Strand.POSITIVE, CoordinateSystem.ONE_BASED, 12345, 12345, "C", "<INS>", 200)));
// breakend variants
// 1 12345 bnd_U C C[2:321682[ 6 PASS SVTYPE=BND;MATEID=bnd_V;EVENT=tra2
Variant bnd = vcfConverter.convertBreakend(vcfConverter.parseContig("chr1"), "bnd_U", 12345, "C", "C[2:321682[", ConfidenceInterval.precise(), "bnd_V", "tra2");
Breakend left = Breakend.of(chr1, "bnd_U", Strand.POSITIVE, CoordinateSystem.ONE_BASED, 12346, 12345);
Breakend right = Breakend.of(chr2, "bnd_V", Strand.POSITIVE, CoordinateSystem.ONE_BASED, 321682, 321681);
assertThat(bnd, equalTo(Variant.of("tra2", left, right, "C", "")));
}
}Multi-allelic sites
Another source of errors and confusion for the unprepared are VCF multi-allelic sites where a single VCF record lists two or more ALT alleles at a position. In order that the variants can be correctly represented they require 'trimming' to their shortest form which, depending on the untrimmed variant can have a different start to what is displayed in the VCF record.
Again there are several, occasionally opposing, schools of thought on how best to trim a variant left, then right ('right shift')
or right, then left ('left shift'). The HGVS recommend right shifting, whereas VCF requires left shifting. This is further
complicated by requiring a 'padding' base, as VCF does for the reference allele, or removing this base. The svart API
provides the VariantTrimmer class to enable users to perform this in a strand-safe way.
class TrimMultipleAlleleTests {
@Test
public void trimMultiAllelicSite() {
// Given the VCF record:
// chr1 225725424 . CTT C,CT 258.06 FreqFilter AC=1,1;AF=0.500,0.500;AN=2;DP=7;ExcessHet=3.0103;FS=0.000;MAX_FREQ=3.1360424;MLEAC=1,1;MLEAF=0.500,0.500;MQ=60.00;QD=29.21;SOR=0.941 GT:AD:DP:GQ:PL 1/2:0,4,3:7:58:275,76,58,106,0,85
Contig chr1 = TestContig.of(1, 249_250_621);
VariantTrimmer leftShiftingTrimmer = VariantTrimmer.leftShiftingTrimmer(VariantTrimmer.retainingCommonBase());
GenomicVariant firstAllele = trim(leftShiftingTrimmer, chr1, "", Strand.POSITIVE, CoordinateSystem.oneBased(), 225_725_424, "CTT", "C");
assertThat(firstAllele, equalTo(GenomicVariant.of(chr1, "", Strand.POSITIVE, CoordinateSystem.oneBased(), 225_725_424, "CTT", "C")));
GenomicVariant secondAllele = trim(leftShiftingTrimmer, chr1, "", Strand.POSITIVE, CoordinateSystem.oneBased(), 225_725_424, "CTT", "CT");
assertThat(secondAllele, equalTo(GenomicVariant.of(chr1, "", Strand.POSITIVE, CoordinateSystem.oneBased(), 225_725_424, "CT", "C")));
}
private GenomicVariant trim(VariantTrimmer variantTrimmer, Contig contig, String id, Strand strand, CoordinateSystem coordinateSystem, int start, String ref, String alt) {
VariantPosition trimmed = variantTrimmer.trim(strand, start, ref, alt);
return GenomicVariant.of(contig, id, strand, coordinateSystem, trimmed.start(), trimmed.ref(), trimmed.alt());
}
}Svart does not however perform full variant normalisation as this should have been performed as part of the VCF creation.
BED
What about BED? BED provides coordinates of a genomic region in
standard genomic coordinates i.e. left open intervals (a.k.a. zero-based) on the positive strand, but can optionally
indicate the strand on which the feature is located in column 6 with a + (forward/positive) or - (reverse/negative) character. The code
below shows how these records can be parsed into GenomicRegion instances.
class BedTests {
@Test
public void parseGenomicRegionFromBedFile() {
// Load the Human GRCh37.13 assembly
GenomicAssembly b37 = GenomicAssembly.GRCh37p13();
// BED uses left-open coordinates, with positions in standard genomic coordinates (i.e. positive strand), with
// the 6th column indicating the strand. Using the example from - https://grch37.ensembl.org/info/website/upload/bed.html
GenomicRegion pos1 = parseBedRecord(b37, "chr7\t127471196\t127472363\tPos1\t0\t+");
GenomicRegion neg1 = parseBedRecord(b37, "chr7\t127475864\t127477031\tNeg1\t0\t-");
assertThat(pos1.contigName(), equalTo("7"));
assertThat(pos1.startOnStrand(Strand.POSITIVE), equalTo(127471196));
assertThat(pos1.endOnStrand(Strand.POSITIVE), equalTo(127472363));
assertThat(pos1.strand(), equalTo(Strand.POSITIVE));
assertThat(neg1.contigName(), equalTo("7"));
// a new instance on the positive strand can be created using neg1.toPositiveStrand(), however this will create
// a new object which may not be wanted in the long term, so to avoid this you can use the start/endOnStrand method
// which will efficiently calculate the result of neg1.withStrand(POSITIVE).start()/.end()
assertThat(neg1.startOnStrand(Strand.POSITIVE), equalTo(127475864));
assertThat(neg1.endOnStrand(Strand.POSITIVE), equalTo(127477031));
assertThat(neg1.strand(), equalTo(Strand.NEGATIVE));
}
private GenomicRegion parseBedRecord(GenomicAssembly genomicAssembly, String bedRecord) {
String[] fields = bedRecord.split("\t");
Contig contig = genomicAssembly.contigByName(fields[0]);
Strand strand = Strand.parseStrand(fields[5]);
int start = strand == Strand.POSITIVE ? Integer.parseInt(fields[1]) : Coordinates.invertCoordinate(CoordinateSystem.zeroBased(), contig, Integer.parseInt(fields[2]));
int end = strand == Strand.POSITIVE ? Integer.parseInt(fields[2]) : Coordinates.invertCoordinate(CoordinateSystem.zeroBased(), contig, Integer.parseInt(fields[1]));
return GenomicRegion.of(contig, strand, CoordinateSystem.zeroBased(), start, end);
}
}