/*
* Copyright 2016 Verily Life Sciences LLC.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package com.verily.genomewarp.utils;
import com.google.common.collect.ComparisonChain;
import com.google.genomics.v1.Position;
import com.google.genomics.v1.Range;
import com.google.genomics.v1.Variant;
import com.google.genomics.v1.VariantCall;
import com.verily.genomewarp.GenomeWarpSerial;
import com.verily.genomewarp.HomologousRangeOuterClass.HomologousRange;
import com.verily.genomewarp.HomologousRangeOuterClass.HomologousRange.RegionType;
import com.verily.genomewarp.HomologousRangeOuterClass.HomologousRange.TargetStrand;
import htsjdk.samtools.util.SequenceUtil;
import htsjdk.samtools.util.StringUtil;
import htsjdk.variant.variantcontext.VariantContext;
import htsjdk.variant.vcf.VCFFileReader;
import java.io.BufferedReader;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
import java.util.HashMap;
import java.util.HashSet;
import java.util.Iterator;
import java.util.List;
import java.util.ListIterator;
import java.util.Map;
import java.util.Set;
import java.util.logging.Level;
import java.util.logging.Logger;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
/**
* Utility functions for performing GenomeWarp
*/
public class GenomeWarpUtils {
// The size in base pairs of each bucket in which to bin HomologousRanges and Variants.
// Given ~4M variants in a whole human genome (~3B base pairs) implies that each bucket will have
// BUCKET_SORT_SIZE / (3e9 / 4e6) = 133 variants in it.
private static final long BUCKET_SORT_SIZE = 100_000L;
public static final Pattern dnaSequence = Pattern.compile("[ACTGactg]+");
public static final int VARIANT_CONTEXT_SIZE = 5;
public static void fail(Logger logger, String message) {
logger.log(Level.SEVERE, message);
throw new RuntimeException(message);
}
private enum GenotypeCategory {
REF_HOMOZYGOUS, // Diploid
HETEROZYGOUS, // Diploid
NON_REF_HOMOZYGOUS, // Diploid
REF, // Haploid
NON_REF // Haploid
}
private static final GenotypeCategory classifyGenotype(List<Integer> genotype) {
boolean isHaploid = genotype.size() < 2;
boolean hasRef = false;
boolean allRef = true;
for (Integer allele : genotype) {
if (allele == 0) {
hasRef = true;
} else {
allRef = false;
}
}
if (allRef) {
// Reference homozygous iff {0,0}, haploid reference if {0}
return isHaploid ? GenotypeCategory.REF : GenotypeCategory.REF_HOMOZYGOUS;
} else if (hasRef) {
return GenotypeCategory.HETEROZYGOUS;
}
// We group all alternative alleles into non-reference for QC purposes.
return isHaploid ? GenotypeCategory.NON_REF : GenotypeCategory.NON_REF_HOMOZYGOUS;
}
public static boolean hasVariation(Variant variant) {
for (VariantCall call : variant.getCallsList()) {
GenotypeCategory category = classifyGenotype(call.getGenotypeList());
if (category != GenotypeCategory.REF_HOMOZYGOUS && category != GenotypeCategory.REF) {
return true;
}
}
return false;
}
/**
* Returns true if the input string is a valid DNA string,
* which we take to mean only containing the characters ACTGactg.
*
* @param inString the input string sequence to test
* @return a boolean indicating whether the input is a DNA sequence
*/
public static boolean isValidDna(String inString) {
final byte[] in = StringUtil.stringToBytes(inString);
for (int i = 0; i < in.length; i++) {
if (!SequenceUtil.isValidBase(in[i])) {
return false;
}
}
return true;
}
/**
* Defines an ordering on variants by comparing on the chromosome then
* on the start position of the variant.
*/
public static class VariantComparator implements Comparator<Variant> {
@Override
public int compare(Variant varA, Variant varB) {
return ComparisonChain.start()
.compare(varA.getReferenceName(), varB.getReferenceName())
.compare(varA.getStart(), varB.getStart())
.compare(varA.getEnd(), varB.getEnd())
.result();
}
}
/**
* Defines an ordering on Range protos by comparing them on the reference
* name, the start, then the end
*/
public static class RangeComparator implements Comparator<Range> {
@Override
public int compare(Range a, Range b) {
return ComparisonChain.start()
.compare(a.getReferenceName(), b.getReferenceName())
.compare(a.getStart(), b.getStart())
.compare(a.getEnd(), b.getEnd())
.result();
}
}
/**
* Defines a ordering on HomologousRange protos by comparing them on the
* query chromosome, then target chromosome
*/
public static class HomologousRangeComparator implements Comparator<HomologousRange> {
@Override
public int compare(HomologousRange a, HomologousRange b) {
return ComparisonChain.start()
.compare(a.getQueryRange(), b.getQueryRange(), new RangeComparator())
.compare(a.getTargetRange(), b.getTargetRange(), new RangeComparator())
.result();
}
}
/**
* Takes a list of variants and returns a new list of variants,
* equal but without the info field present.
*
* @param variants the input list of variants
* @return a list of variants without the info field
*/
public static List<Variant> removeInfoField(List<Variant> variants) {
List<Variant> toReturn = new ArrayList<>();
for (Variant currVariant : variants) {
// Pulls out just the fields we need from verbose variant record
toReturn.add(Variant.newBuilder()
.setReferenceName(currVariant.getReferenceName())
.setStart(currVariant.getStart())
.setEnd(currVariant.getEnd())
.addAllNames(currVariant.getNamesList())
.setReferenceBases(currVariant.getReferenceBases())
.addAllAlternateBases(currVariant.getAlternateBasesList())
.setQuality(currVariant.getQuality())
.addAllFilter(currVariant.getFilterList())
.addAllCalls(currVariant.getCallsList())
.build());
}
return toReturn;
}
/**
* Returns a position proto encoding the input reference name
* and which chunk of size BUCKET_SORT_SIZE the start belongs in.
*
* @param refName the reference name of the returned position
* @param start the start coordinate used to decide the chunk index
* @return position
*/
public static Position getPosition(String refName, long start) {
return Position.newBuilder().setReferenceName(refName)
.setPosition(start / BUCKET_SORT_SIZE).build();
}
private static Map<Position, List<HomologousRange>> associateRegionsWithPosition(
List<HomologousRange> regions) {
Map<Position, List<HomologousRange>> positionToRegions = new HashMap<>();
for (HomologousRange currRegion : regions) {
String reference = currRegion.getQueryRange().getReferenceName();
long startBucket = currRegion.getQueryRange().getStart() / BUCKET_SORT_SIZE;
long endBucket = currRegion.getQueryRange().getEnd() / BUCKET_SORT_SIZE;
for (long i = startBucket; i <= endBucket; i++) {
Position bucket = Position.newBuilder().setReferenceName(reference).setPosition(i).build();
List<HomologousRange> regionList = null;
if (positionToRegions.containsKey(bucket)) {
regionList = positionToRegions.get(bucket);
} else {
regionList = new ArrayList<>();
}
regionList.add(currRegion);
positionToRegions.put(bucket, regionList);
}
}
return positionToRegions;
}
/**
* Given a BufferedReader encoding a file of variants, asssociate the string
* representation of the variants with a region. This association is done by
* mapping each region to all variants which share the same reference name
* (chromosome) and whose start coordinate is within the region.
*
* @param regions the list of regions to use as the map's keys
* @param vcfFile a BufferedReader encoding the input variants
* @return a map from regions to variant strings
*/
public static Map<HomologousRange, List<String>> associateVariantsWithRange(
List<HomologousRange> regions, BufferedReader vcfFile) throws IOException {
Map<HomologousRange, List<String>> toReturn = new HashMap<>();
// First associate each Region with one or more Positions
Map<Position, List<HomologousRange>> positionToRegions =
associateRegionsWithPosition(regions);
// Add initial empty lists
for (HomologousRange currRegion : regions) {
toReturn.put(currRegion, new ArrayList<String>());
}
// We iterate through the variants. For each variant we get the corresponding
// bucket, then check through the Regions in the bucket that suit our variant
//
// Furthermore, note that regions without any variants are still included, since
// we previously added every region with an initially empty list of variants to toReturn.
String line;
while ((line = vcfFile.readLine()) != null) {
if (line.length() == 0 || line.charAt(0) == '#') {
// Header or blank line
continue;
}
String[] vcfArray = line.split("\t");
String refName = vcfArray[0];
Long start = 0L;
try {
start = Long.parseLong(vcfArray[1]) - 1;
} catch (NumberFormatException ex) {
continue;
}
List<HomologousRange> correspondingRegions =
positionToRegions.get(getPosition(refName, start));
if (correspondingRegions == null) {
// Variant occurs outside the confidently-called regions. These cannot be accurately
// transformed, and so are ignored.
continue;
}
for (HomologousRange currRegion : correspondingRegions) {
if (start >= currRegion.getQueryRange().getStart()
&& start < currRegion.getQueryRange().getEnd()) {
List<String> vcfList = toReturn.get(currRegion);
vcfList.add(line);
toReturn.put(currRegion, vcfList);
}
}
}
return toReturn;
}
/**
* Duplication of associateVariantsWithRange. Required as we need to be able to deal
* with variants formatted as strings as well as in true variant format.
*
* @param regions the list of regions to use as the map's keys
* @param variants the input variants formatted as variant protos.
* @return a map from regions to variant protos
*/
public static Map<HomologousRange, List<Variant>> associateVariantsWithRange(
List<HomologousRange> regions, List<Variant> variants) {
Map<HomologousRange, List<Variant>> toReturn = new HashMap<>();
Map<Position, List<HomologousRange>> positionToRegions = associateRegionsWithPosition(regions);
// Add initial empty lists
for (HomologousRange currRegion : regions) {
toReturn.put(currRegion, new ArrayList<Variant>());
}
for (Variant currVariant : variants) {
String refName = currVariant.getReferenceName();
Long start = currVariant.getStart();
List<HomologousRange> correspondingRegions =
positionToRegions.get(getPosition(refName, start));
if (correspondingRegions == null) {
continue;
}
for (HomologousRange currRegion : correspondingRegions) {
if (start >= currRegion.getQueryRange().getStart()
&& start < currRegion.getQueryRange().getEnd()) {
List<Variant> variantList = toReturn.get(currRegion);
if (variantList == null) {
variantList = new ArrayList<>();
}
variantList.add(currVariant);
toReturn.put(currRegion, variantList);
}
}
}
return toReturn;
}
public static String homologousToString(HomologousRange range) {
Range queryRange = range.getQueryRange();
Range targetRange = range.getTargetRange();
return queryRange.getReferenceName() + "\t" + queryRange.getStart() + "\t" + queryRange.getEnd()
+ "\t" + targetRange.getReferenceName() + "\t" + targetRange.getStart() + "\t"
+ targetRange.getEnd() + "\t"
+ (range.getTargetStrand() == TargetStrand.POSITIVE_STRAND ? "+" : "-") + "\t"
+ range.getRegionType().toString();
}
/**
* Parses a HomologousRange proto from an input string. This string must be tab separated with
* each column defined as follows.
*
* <p> Col. 1: The query reference name (chromosome) <br>
* Col. 2: The start position of the query range (long) <br>
* Col. 3: The end position of the query range (long) <br>
* Col. 4: The target reference name (chromosome) <br>
* Col. 5: The start position of the target range (long) <br>
* Col. 6: The end position of the target range (long) <br>
* Col. 7: The strand (either + or -) <br>
* Col. 8: The region type, encoded as one of IDENTICAL, ALIGNMENT_REQUIRED, or
* MISMATCHED_BASES
*
* @param range the string encoding of the range
* @return a HomologousRange parsed from the string
*/
public static HomologousRange homologousFromString(String range) {
String[] parts = range.split("\t");
if (parts.length != 8) {
throw new IllegalArgumentException("input string must have 9 sections separated by tabs");
}
Long qStart = 0L, qEnd = 0L, tStart = 0L, tEnd = 0L;
try {
qStart = Long.parseLong(parts[1]);
qEnd = Long.parseLong(parts[2]);
tStart = Long.parseLong(parts[4]);
tEnd = Long.parseLong(parts[5]);
} catch (NumberFormatException ex) {
throw new IllegalArgumentException("Region.fromString: " + ex.getMessage());
}
boolean strand = parts[6].equals("+");
return HomologousRange.newBuilder()
.setQueryRange(Range.newBuilder().setReferenceName(parts[0]).setStart(qStart).setEnd(qEnd))
.setTargetRange(Range.newBuilder().setReferenceName(parts[3]).setStart(tStart).setEnd(tEnd))
.setTargetStrand(strand ? TargetStrand.POSITIVE_STRAND : TargetStrand.NEGATIVE_STRAND)
.setRegionType(RegionType.valueOf(parts[7]))
.build();
}
public static<T> T nextOrNull(Iterator<T> it) {
if (it.hasNext()) {
return it.next();
}
return null;
}
}
