if (vc_input.isFiltered()) {
vcfWriter.add(vc_input);
return 1;
}
BeagleFeature beagleR2Feature = tracker.getFirstValue(beagleR2);
BeagleFeature beagleProbsFeature = tracker.getFirstValue(beagleProbs);
BeagleFeature beaglePhasedFeature = tracker.getFirstValue(beaglePhased);
// ignore places where we don't have a variant
if ( beagleR2Feature == null || beagleProbsFeature == null || beaglePhasedFeature == null)
{
vcfWriter.add(vc_input);
return 1;
}
// get reference base for current position
byte refByte = ref.getBase();
// make new Genotypes based on Beagle results
GenotypesContext genotypes = GenotypesContext.create(vc_input.getGenotypes().size());
// for each genotype, create a new object with Beagle information on it
int numGenotypesChangedByBeagle = 0;
Integer alleleCountH = 0, chrCountH = 0;
Double alleleFrequencyH = 0.0;
int beagleVarCounts = 0;
GenotypesContext hapmapGenotypes = null;
if (vc_comp != null) {
hapmapGenotypes = vc_comp.getGenotypes();
}
for ( final Genotype g : vc_input.getGenotypes() ) {
boolean genotypeIsPhased = true;
String sample = g.getSampleName();
// If we have a Hapmap (comp) ROD, compute Hapmap AC, AN and AF
// use sample as key into genotypes structure
if (vc_comp != null) {
if (vc_input.getGenotypes().containsSample(sample) && hapmapGenotypes.containsSample(sample)) {
Genotype hapmapGenotype = hapmapGenotypes.get(sample);
if (hapmapGenotype.isCalled()){
chrCountH += 2;
if (hapmapGenotype.isHet()) {
alleleCountH += 1;
} else if (hapmapGenotype.isHomVar()) {
alleleCountH += 2;
}
}
}
}
ArrayList<String> beagleProbabilities = beagleProbsFeature.getProbLikelihoods().get(sample);
ArrayList<String> beagleGenotypePairs = beaglePhasedFeature.getGenotypes().get(sample);
// original alleles at this genotype
Allele originalAlleleA = g.getAllele(0);
Allele originalAlleleB = (g.getAlleles().size() == 2) ? g.getAllele(1) : g.getAllele(0); // hack to deal with no-call genotypes