zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Moving UnifiedGenotyper and VariantAnnotator over to VariantContext system. 

Removing obsolete genotyping classes.
First stage of removing dependence on old Genotype class.
More changes to come.



git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@2960 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
ebanks 2010-03-09 03:41:07 +00:00
parent 02f48b6457
commit 9f3b99c11b
86 changed files with 748 additions and 2494 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
java/src/org/broadinstitute/sting/gatk/contexts/variantcontext/InferredGeneticContext.java
View File

@ -11,7 +11,7 @@ import java.util.*;
* @author depristo
*/
final class InferredGeneticContext {
public static final double NO_NEG_LOG_10PERROR = Double.MAX_VALUE; // todo -- is this really safe?
public static final double NO_NEG_LOG_10PERROR = -1.0;
private double negLog10PError = NO_NEG_LOG_10PERROR;
private String name = null;

11
java/src/org/broadinstitute/sting/gatk/contexts/variantcontext/MutableGenotype.java
View File

@ -12,11 +12,18 @@ public class MutableGenotype extends Genotype {
super(parent.getSampleName(), parent.getAlleles(), parent.getNegLog10PError(), parent.getFilters(), parent.getAttributes(), parent.genotypesArePhased());
}
// todo -- add rest of genotype constructors here
public MutableGenotype(String sampleName, List<Allele> alleles, double negLog10PError, Set<String> filters, Map<String, ?> attributes, boolean genotypesArePhased) {
super(sampleName, alleles, negLog10PError, filters, attributes, genotypesArePhased);
}
public MutableGenotype(String sampleName, List<Allele> alleles, double negLog10PError) {
super(sampleName, alleles, negLog10PError);
}
public MutableGenotype(String sampleName, List<Allele> alleles) {
super(sampleName, alleles);
}
public Genotype unmodifiableGenotype() {
return new Genotype(getSampleName(), getAlleles(), getNegLog10PError(), getFilters(), getAttributes(), genotypesArePhased());
}
@ -24,7 +31,7 @@ public class MutableGenotype extends Genotype {
/**
*
* @param alleles
* @param alleles list of alleles
*/
public void setAlleles(List<Allele> alleles) {
this.alleles.clear();

2
java/src/org/broadinstitute/sting/gatk/contexts/variantcontext/VariantContext.java
View File

@ -463,7 +463,7 @@ public class VariantContext {
public double getNegLog10PError() { return commonInfo.getNegLog10PError(); }
public double getPhredScaledQual() { return commonInfo.getPhredScaledQual(); }
public Map<String, Object> getAttributes() { return commonInfo.getAttributes(); }
public Map<String, Object> getAttributes() { return commonInfo.getAttributes(); }
public boolean hasAttribute(String key) { return commonInfo.hasAttribute(key); }
public Object getAttribute(String key) { return commonInfo.getAttribute(key); }

1
java/src/org/broadinstitute/sting/gatk/io/storage/GLFGenotypeWriterStorage.java
View File

@ -54,7 +54,6 @@ public class GLFGenotypeWriterStorage extends GenotypeWriterStorage<GLFGenotypeW
* Merges the stream backing up this temporary storage into the target.
* @param target Target stream for the temporary storage. May not be null.
*/
@Override
public void mergeInto(GLFGenotypeWriter target) {
GLFReader reader = new GLFReader(file);
while ( reader.hasNext() ) {

1
java/src/org/broadinstitute/sting/gatk/io/storage/GeliBinaryGenotypeWriterStorage.java
View File

@ -54,7 +54,6 @@ public class GeliBinaryGenotypeWriterStorage extends GenotypeWriterStorage<GeliG
* Merges the stream backing up this temporary storage into the target.
* @param target Target stream for the temporary storage. May not be null.
*/
@Override
public void mergeInto(GeliGenotypeWriter target) {
GeliFileReader reader = new GeliFileReader(file);
while ( reader.hasNext() )

1
java/src/org/broadinstitute/sting/gatk/io/storage/GeliTextGenotypeWriterStorage.java
View File

@ -53,7 +53,6 @@ public class GeliTextGenotypeWriterStorage extends GenotypeWriterStorage<GeliGen
* Merges the stream backing up this temporary storage into the target.
* @param target Target stream for the temporary storage. May not be null.
*/
@Override
public void mergeInto(GeliGenotypeWriter target) {
GeliFileReader reader = new GeliFileReader(file);
while ( reader.hasNext() )

18
java/src/org/broadinstitute/sting/gatk/io/storage/GenotypeWriterStorage.java
View File

@ -26,11 +26,11 @@
package org.broadinstitute.sting.gatk.io.storage;
import java.io.*;
import java.util.List;
import java.util.Set;
import java.util.HashSet;
import org.broadinstitute.sting.gatk.io.stubs.GenotypeWriterStub;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.utils.genotype.*;
import org.broadinstitute.sting.utils.genotype.vcf.*;
import org.broadinstitute.sting.utils.SampleUtils;
@ -76,20 +76,8 @@ public abstract class GenotypeWriterStorage<T extends GenotypeWriter> implements
GenotypeWriterFactory.writeHeader(writer, stub.getSAMFileHeader(), samples, new HashSet<VCFHeaderLine>());
}
public void addGenotypeCall(Genotype call) {
writer.addGenotypeCall(call);
}
public void addNoCall(int position) {
writer.addNoCall(position);
}
public void addMultiSampleCall(List<Genotype> genotypes, VariationCall variation) {
writer.addMultiSampleCall(genotypes, variation);
}
public boolean supportsMultiSample() {
return writer.supportsMultiSample();
public void addCall(VariantContext vc) {
writer.addCall(vc);
}
public void close() {

29
java/src/org/broadinstitute/sting/gatk/io/stubs/GenotypeWriterStub.java
View File

@ -27,13 +27,11 @@ package org.broadinstitute.sting.gatk.io.stubs;
import java.io.File;
import java.io.PrintStream;
import java.util.List;
import org.broadinstitute.sting.gatk.io.OutputTracker;
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.utils.genotype.GenotypeWriter;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.genotype.VariationCall;
import org.broadinstitute.sting.utils.genotype.GenotypeWriterFactory;
import net.sf.samtools.SAMFileHeader;
@ -131,29 +129,8 @@ public abstract class GenotypeWriterStub<T extends GenotypeWriter> implements St
/**
* @{inheritDoc}
*/
public void addGenotypeCall(Genotype call) {
outputTracker.getStorage(this).addGenotypeCall(call);
}
/**
* @{inheritDoc}
*/
public void addNoCall(int position) {
outputTracker.getStorage(this).addNoCall(position);
}
/**
* @{inheritDoc}
*/
public void addMultiSampleCall(List<Genotype> genotypes, VariationCall variation) {
outputTracker.getStorage(this).addMultiSampleCall(genotypes, variation);
}
/**
* @{inheritDoc}
*/
public boolean supportsMultiSample() {
return outputTracker.getStorage(this).supportsMultiSample();
public void addCall(VariantContext vc) {
outputTracker.getStorage(this).addCall(vc);
}
/**

55
java/src/org/broadinstitute/sting/gatk/refdata/RodGeliText.java
View File

@ -28,17 +28,14 @@ package org.broadinstitute.sting.gatk.refdata;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.GenomeLocParser;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.genotype.DiploidGenotype;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.genotype.VariantBackedByGenotype;
import org.broadinstitute.sting.utils.genotype.geli.GeliGenotypeCall;
import org.broadinstitute.sting.utils.genotype.Variation;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
public class RodGeliText extends BasicReferenceOrderedDatum implements VariationRod, VariantBackedByGenotype {
public class RodGeliText extends BasicReferenceOrderedDatum implements Variation {
public enum Genotype_Strings {
AA, AC, AG, AT, CC, CG, CT, GG, GT, TT
}
@ -119,7 +116,6 @@ public class RodGeliText extends BasicReferenceOrderedDatum implements Variation
*
* @return the reference base or bases, as a string
*/
@Override
public String getReference() {
return String.valueOf(this.refBase);
}
@ -130,7 +126,6 @@ public class RodGeliText extends BasicReferenceOrderedDatum implements Variation
*
* @return the log based error estimate
*/
@Override
public double getNegLog10PError() {
return Math.abs(lodBtr);
}
@ -143,7 +138,6 @@ public class RodGeliText extends BasicReferenceOrderedDatum implements Variation
*
* @return an alternate allele list
*/
@Override
public List<String> getAlternateAlleleList() {
List<String> list = new ArrayList<String>();
for (char base : bestGenotype.toCharArray())
@ -159,7 +153,6 @@ public class RodGeliText extends BasicReferenceOrderedDatum implements Variation
*
* @return an alternate allele list
*/
@Override
public List<String> getAlleleList() {
List<String> list = new ArrayList<String>();
if (this.bestGenotype.contains(getReference())) list.add(getReference());
@ -197,15 +190,13 @@ public class RodGeliText extends BasicReferenceOrderedDatum implements Variation
*
* @return VariantFrequency with the stored frequency
*/
@Override
public double getNonRefAlleleFrequency() {
return 1.0;
}
/** @return the VARIANT_TYPE of the current variant */
@Override
public VARIANT_TYPE getType() {
return VARIANT_TYPE.SNP;
public Variation.VARIANT_TYPE getType() {
return Variation.VARIANT_TYPE.SNP;
}
public boolean isSNP() {
@ -232,7 +223,6 @@ public class RodGeliText extends BasicReferenceOrderedDatum implements Variation
*
* @return a char, representing the alternate base
*/
@Override
public char getAlternativeBaseForSNP() {
if (!this.isSNP()) throw new IllegalStateException("we're not a SNP");
// we know that if we're a SNP, the alt is a single base
@ -247,7 +237,6 @@ public class RodGeliText extends BasicReferenceOrderedDatum implements Variation
*
* @return a char, representing the alternate base
*/
@Override
public char getReferenceForSNP() {
if (!isSNP()) throw new IllegalStateException("This site is not a SNP");
// we know that if we're a SNP, the reference is a single base
@ -357,40 +346,4 @@ public class RodGeliText extends BasicReferenceOrderedDatum implements Variation
this.lodBtr = (bestLikelihood - refLikelihood);
this.lodBtnb = (bestLikelihood - nextBestLikelihood);
}
/**
* get the genotype
*
* @return a map in lexigraphical order of the genotypes
*/
@Override
public Genotype getCalledGenotype() {
return new GeliGenotypeCall(Character.toString(refBase), getLocation(), bestGenotype, lodBtnb);
}
/**
* get the likelihoods
*
* @return an array in lexigraphical order of the likelihoods
*/
@Override
public List<Genotype> getGenotypes() {
List<Genotype> ret = new ArrayList<Genotype>();
ret.add(new GeliGenotypeCall(Character.toString(refBase), getLocation(), bestGenotype, lodBtnb));
return ret;
}
/**
* do we have the specified genotype? not all backedByGenotypes
* have all the genotype data.
*
* @param x the genotype
*
* @return true if available, false otherwise
*/
@Override
public boolean hasGenotype(DiploidGenotype x) {
return (x.toString().equals(this.getAltBasesFWD()));
}
}

54
java/src/org/broadinstitute/sting/gatk/refdata/RodVCF.java
View File

@ -2,9 +2,7 @@ package org.broadinstitute.sting.gatk.refdata;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.StingException;
import org.broadinstitute.sting.utils.genotype.DiploidGenotype;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.genotype.VariantBackedByGenotype;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.vcf.*;
import java.io.File;
@ -20,7 +18,7 @@ import java.util.*;
* <p/>
* An implementation of the ROD for VCF.
*/
public class RodVCF extends BasicReferenceOrderedDatum implements VariationRod, VariantBackedByGenotype, Iterator<RodVCF> {
public class RodVCF extends BasicReferenceOrderedDatum implements Variation, Iterator<RodVCF> {
public VCFReader getReader() {
return mReader;
}
@ -107,7 +105,7 @@ public class RodVCF extends BasicReferenceOrderedDatum implements VariationRod,
}
/** @return the VARIANT_TYPE of the current variant */
public VARIANT_TYPE getType() {
public Variation.VARIANT_TYPE getType() {
assertNotNull();
return mCurrentRecord.getType();
}
@ -130,8 +128,9 @@ public class RodVCF extends BasicReferenceOrderedDatum implements VariationRod,
/**
* are we a novel site? Is there a DBSNP identifier
* or a hapmap entry for the site?
*
* @return true if this site is novel
*/
public boolean isNovel() {
assertNotNull();
return mCurrentRecord.isNovel();
@ -173,7 +172,9 @@ public class RodVCF extends BasicReferenceOrderedDatum implements VariationRod,
return mCurrentRecord.getReference();
}
/** are we bi-allelic? */
/** are we bi-allelic?
* @return true if we're bi-allelic
*/
public boolean isBiallelic() {
assertNotNull();
return mCurrentRecord.isBiallelic();
@ -270,28 +271,6 @@ public class RodVCF extends BasicReferenceOrderedDatum implements VariationRod,
return mCurrentRecord.hasGenotypeData();
}
/**
* get the genotype
*
* // todo -- WTF is this? This is a deeply unsafe call
*
* @return a map in lexigraphical order of the genotypes
*/
public Genotype getCalledGenotype() {
assertNotNull();
return mCurrentRecord.getCalledGenotype();
}
/**
* get the genotypes
*
* @return a list of the genotypes
*/
public List<Genotype> getGenotypes() {
assertNotNull();
return mCurrentRecord.getGenotypes();
}
/**
* get the genotypes
*
@ -306,25 +285,12 @@ public class RodVCF extends BasicReferenceOrderedDatum implements VariationRod,
* Returns the genotype associated with sample, or null if the genotype is missing
*
* @param sampleName the name of the sample genotype to fetch
* @return
* @return genotype record
*/
public Genotype getGenotype(final String sampleName) {
public VCFGenotypeRecord getGenotype(final String sampleName) {
return mCurrentRecord.getGenotype(sampleName);
}
/**
* do we have the specified genotype? not all backedByGenotypes
* have all the genotype data.
*
* @param x the genotype
*
* @return true if available, false otherwise
*/
public boolean hasGenotype(DiploidGenotype x) {
assertNotNull();
return mCurrentRecord.hasGenotype(x);
}
public String[] getSampleNames() {
assertNotNull();
return mCurrentRecord.getSampleNames();

97
java/src/org/broadinstitute/sting/gatk/refdata/VariantContextAdaptors.java
View File

@ -1,10 +1,9 @@
package org.broadinstitute.sting.gatk.refdata;
import org.broadinstitute.sting.utils.genotype.vcf.*;
import org.broadinstitute.sting.utils.StingException;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.GenomeLocParser;
import org.broadinstitute.sting.utils.genotype.CalledGenotype;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.Allele;
import org.broadinstitute.sting.gatk.contexts.variantcontext.Genotype;
@ -155,7 +154,7 @@ public class VariantContextAdaptors {
}
Set<String> filters = vcf.isFiltered() ? new HashSet<String>(Arrays.asList(vcf.getFilteringCodes())) : null;
Map<String, String> attributes = vcf.getInfoValues();
Map<String, String> attributes = new HashMap<String, String>(vcf.getInfoValues());
attributes.put("ID", vcf.getID());
// add all of the alt alleles
@ -181,8 +180,6 @@ public class VariantContextAdaptors {
genotypeAlleles.add(a);
}
double pError = vcfG.getNegLog10PError() == VCFGenotypeRecord.MISSING_GENOTYPE_QUALITY ? Genotype.NO_NEG_LOG_10PERROR : vcfG.getNegLog10PError();
Map<String, String> fields = new HashMap<String, String>();
for ( Map.Entry<String, String> e : vcfG.getFields().entrySet() ) {
// todo -- fixme if we put GQ and GF into key itself
@ -194,7 +191,7 @@ public class VariantContextAdaptors {
if ( vcfG.isFiltered() ) // setup the FL genotype filter fields
genotypeFilters.addAll(Arrays.asList(vcfG.getFields().get(VCFGenotypeRecord.GENOTYPE_FILTER_KEY).split(";")));
Genotype g = new Genotype(vcfG.getSampleName(), genotypeAlleles, pError, genotypeFilters, fields, vcfG.getPhaseType() == VCFGenotypeRecord.PHASE.PHASED);
Genotype g = new Genotype(vcfG.getSampleName(), genotypeAlleles, vcfG.getNegLog10PError(), genotypeFilters, fields, vcfG.getPhaseType() == VCFGenotypeRecord.PHASE.PHASED);
genotypes.put(g.getSampleName(), g);
}
@ -216,16 +213,20 @@ public class VariantContextAdaptors {
return new VCFHeader(hInfo == null ? new HashSet<VCFHeaderLine>() : hInfo, names);
}
public static VCFRecord toVCF(VariantContext vc, char previousRefBase) {
public static VCFRecord toVCF(VariantContext vc, char vcfRefBase) {
return toVCF(vc, vcfRefBase, null, true);
}
public static VCFRecord toVCF(VariantContext vc, char vcfRefBase, List<String> vcfGenotypeAttributeKeys, boolean filtersWereAppliedToContext) {
// deal with the reference
String referenceBases = new String(vc.getReference().getBases());
String contig = vc.getLocation().getContig();
long position = vc.getLocation().getStart();
String ID = vc.hasAttribute("ID") ? vc.getAttributeAsString("ID") : ".";
String ID = vc.hasAttribute("ID") ? vc.getAttributeAsString("ID") : VCFRecord.EMPTY_ID_FIELD;
double qual = vc.hasNegLog10PError() ? vc.getPhredScaledQual() : -1;
String filters = vc.isFiltered() ? Utils.join(";", Utils.sorted(vc.getFilters())) : VCFRecord.PASSES_FILTERS;
String filters = vc.isFiltered() ? Utils.join(";", Utils.sorted(vc.getFilters())) : (filtersWereAppliedToContext ? VCFRecord.PASSES_FILTERS : VCFRecord.UNFILTERED);
// TODO - fixme: temporarily using Strings as keys because Alleles aren't hashing correctly
Map<String, VCFGenotypeEncoding> alleleMap = new HashMap<String, VCFGenotypeEncoding>();
@ -247,21 +248,21 @@ public class VariantContextAdaptors {
if ( a.isNull() ) {
if ( a.isReference() ) {
// ref, where alt is insertion
encoding = new VCFGenotypeEncoding(Character.toString(previousRefBase));
encoding = new VCFGenotypeEncoding(Character.toString(vcfRefBase));
} else {
// non-ref deletion
encoding = new VCFGenotypeEncoding("D" + Integer.toString(referenceBases.length()));
}
} else if ( a.isReference() ) {
// ref, where alt is deletion
encoding = new VCFGenotypeEncoding(Character.toString(previousRefBase));
encoding = new VCFGenotypeEncoding(Character.toString(vcfRefBase));
} else {
// non-ref insertion
encoding = new VCFGenotypeEncoding("I" + alleleString);
}
} else if ( vc.getType() == VariantContext.Type.NO_VARIATION ) {
// ref
encoding = new VCFGenotypeEncoding(Character.toString(previousRefBase));
encoding = new VCFGenotypeEncoding(Character.toString(vcfRefBase));
} else {
// ref or alt for snp
encoding = new VCFGenotypeEncoding(alleleString);
@ -274,26 +275,54 @@ public class VariantContextAdaptors {
alleleMap.put(a.toString(), encoding);
}
List<String> vcfGenotypeAttributeKeys = new ArrayList<String>(Arrays.asList(VCFGenotypeRecord.GENOTYPE_KEY));
List<String> vcGenotypeKeys = calcVCFGenotypeKeys(vc);
if ( vc.hasGenotypes() ) vcfGenotypeAttributeKeys.addAll(vcGenotypeKeys);
if ( vcfGenotypeAttributeKeys == null ) {
vcfGenotypeAttributeKeys = new ArrayList<String>();
if ( vc.hasGenotypes() ) {
vcfGenotypeAttributeKeys.add(VCFGenotypeRecord.GENOTYPE_KEY);
vcfGenotypeAttributeKeys.addAll(calcVCFGenotypeKeys(vc));
}
}
String genotypeFormatString = Utils.join(VCFRecord.GENOTYPE_FIELD_SEPERATOR, vcfGenotypeAttributeKeys);
List<VCFGenotypeRecord> genotypeObjects = new ArrayList<VCFGenotypeRecord>(vc.getGenotypes().size());
for ( Genotype g : vc.getGenotypesSortedByName() ) {
List<VCFGenotypeEncoding> encodings = new ArrayList<VCFGenotypeEncoding>(g.getPloidy());
for ( Allele a : g.getAlleles() ) {
// TODO -- REMOVE ME ONCE INTEGRATION TESTS PASS!!!
ArrayList<Allele> temporaryList = new ArrayList<Allele>(g.getAlleles());
Collections.sort(temporaryList, new Comparator<Allele>() {
public int compare(Allele a1, Allele a2) {
return a1.toString().charAt(0) - a2.toString().charAt(0);
}
});
for ( Allele a : temporaryList ) {
//for ( Allele a : g.getAlleles() ) {
encodings.add(alleleMap.get(a.toString()));
}
VCFGenotypeRecord.PHASE phasing = g.genotypesArePhased() ? VCFGenotypeRecord.PHASE.PHASED : VCFGenotypeRecord.PHASE.UNPHASED;
VCFGenotypeRecord vcfG = new VCFGenotypeRecord(g.getSampleName(), encodings, phasing);
if ( ! g.isNoCall() ) {
for ( String key : vcGenotypeKeys ) {
String val = key.equals(VCFGenotypeRecord.GENOTYPE_QUALITY_KEY) ? String.format("%.2f", g.getPhredScaledQual()) : g.getAttribute(key).toString();
vcfG.setField(key, val);
for ( String key : vcfGenotypeAttributeKeys ) {
if ( key.equals(VCFGenotypeRecord.GENOTYPE_KEY) )
continue;
Object val = g.getAttribute(key);
// some exceptions
if ( key.equals(VCFGenotypeRecord.GENOTYPE_QUALITY_KEY) ) {
if ( MathUtils.compareDoubles(g.getNegLog10PError(), Genotype.NO_NEG_LOG_10PERROR) == 0 )
val = VCFGenotypeRecord.MISSING_GENOTYPE_QUALITY;
else
val = Math.min(g.getPhredScaledQual(), VCFGenotypeRecord.MAX_QUAL_VALUE);
} else if ( key.equals(VCFGenotypeRecord.DEPTH_KEY) && val == null ) {
ReadBackedPileup pileup = (ReadBackedPileup)g.getAttribute(CalledGenotype.READBACKEDPILEUP_ATTRIBUTE_KEY);
if ( pileup != null )
val = pileup.size();
}
String outputValue = formatVCFField(key, val);
if ( outputValue != null )
vcfG.setField(key, outputValue);
}
genotypeObjects.add(vcfG);
@ -302,13 +331,27 @@ public class VariantContextAdaptors {
Map<String, String> infoFields = new HashMap<String, String>();
for ( Map.Entry<String, Object> elt : vc.getAttributes().entrySet() ) {
String key = elt.getKey();
String val = elt.getValue().toString();
if ( ! key.equals("ID") ) {
infoFields.put(key, val);
}
if ( key.equals("ID") )
continue;
String outputValue = formatVCFField(key, elt.getValue());
if ( outputValue != null )
infoFields.put(key, outputValue);
}
return new VCFRecord(Character.toString(previousRefBase), contig, position, ID, vcfAltAlleles, qual, filters, infoFields, genotypeFormatString, genotypeObjects);
return new VCFRecord(Character.toString(vcfRefBase), contig, position, ID, vcfAltAlleles, qual, filters, infoFields, genotypeFormatString, genotypeObjects);
}
private static String formatVCFField(String key, Object val) {
String result;
if ( val == null )
result = VCFGenotypeRecord.getMissingFieldValue(key);
else if ( val instanceof Double )
result = String.format("%.2f", val);
else
result = val.toString();
return result;
}
private static List<String> calcVCFGenotypeKeys(VariantContext vc) {

6
java/src/org/broadinstitute/sting/gatk/walkers/annotator/Alignability.java
View File

@ -2,11 +2,9 @@ package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.refdata.TabularROD;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import java.util.Map;
@ -17,7 +15,7 @@ public class Alignability implements VariantAnnotation {
public String annotate(RefMetaDataTracker tracker,
ReferenceContext ref,
Map<String, StratifiedAlignmentContext> stratifiedContexts,
Variation variation)
VariantContext vc)
{
TabularROD record = (TabularROD)(tracker.lookup("alignability", null));
int value;

37
java/src/org/broadinstitute/sting/gatk/walkers/annotator/AlleleBalance.java
View File

@ -2,63 +2,52 @@ package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.*;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.genotype.*;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import java.util.List;
import java.util.Map;
public class AlleleBalance extends StandardVariantAnnotation {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation) {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, VariantContext vc) {
if ( !variation.isBiallelic() || !(variation instanceof VariantBackedByGenotype) )
if ( !vc.isBiallelic() || !vc.isSNP() )
return null;
final List<Genotype> genotypes = ((VariantBackedByGenotype)variation).getGenotypes();
final Map<String, Genotype> genotypes = vc.getGenotypes();
if ( genotypes == null || genotypes.size() == 0 )
return null;
double ratio = 0.0;
double totalWeights = 0.0;
for ( Genotype genotype : genotypes ) {
for ( Map.Entry<String, Genotype> genotype : genotypes.entrySet() ) {
// we care only about het calls
if ( !genotype.isHet() )
if ( !genotype.getValue().isHet() )
continue;
if ( !(genotype instanceof SampleBacked) )
continue;
String sample = ((SampleBacked)genotype).getSampleName();
StratifiedAlignmentContext context = stratifiedContexts.get(sample);
StratifiedAlignmentContext context = stratifiedContexts.get(genotype.getKey());
if ( context == null )
continue;
final String genotypeStr = genotype.getBases().toUpperCase();
if ( genotypeStr.length() != 2 )
return null;
final String bases = new String(context.getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).getBasePileup().getBases()).toUpperCase();
if ( bases.length() == 0 )
return null;
char a = genotypeStr.charAt(0);
char b = genotypeStr.charAt(1);
int aCount = Utils.countOccurrences(a, bases);
int bCount = Utils.countOccurrences(b, bases);
char refChr = vc.getReference().toString().charAt(0);
char altChr = vc.getAlternateAllele(0).toString().charAt(0);
int refCount = (a == ref.getBase() ? aCount : bCount);
int altCount = (a == ref.getBase() ? bCount : aCount);
int refCount = Utils.countOccurrences(refChr, bases);
int altCount = Utils.countOccurrences(altChr, bases);
// sanity check
if ( refCount + altCount == 0 )
continue;
// weight the allele balance by genotype quality so that e.g. mis-called homs don't affect the ratio too much
ratio += genotype.getNegLog10PError() * ((double)refCount / (double)(refCount + altCount));
totalWeights += genotype.getNegLog10PError();
ratio += genotype.getValue().getNegLog10PError() * ((double)refCount / (double)(refCount + altCount));
totalWeights += genotype.getValue().getNegLog10PError();
}
// make sure we had a het genotype

4
java/src/org/broadinstitute/sting/gatk/walkers/annotator/DepthOfCoverage.java
View File

@ -2,8 +2,8 @@ package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.vcf.VCFRecord;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
@ -12,7 +12,7 @@ import java.util.Map;
public class DepthOfCoverage extends StandardVariantAnnotation {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation) {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, VariantContext vc) {
int depth = 0;
for ( String sample : stratifiedContexts.keySet() )
depth += stratifiedContexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).size();

22
java/src/org/broadinstitute/sting/gatk/walkers/annotator/HardyWeinberg.java
View File

@ -2,12 +2,12 @@ package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.*;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.genotype.*;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import org.broadinstitute.sting.utils.genotype.HardyWeinbergCalculation;
import org.broadinstitute.sting.utils.QualityUtils;
import java.util.List;
import java.util.Map;
@ -17,31 +17,31 @@ public class HardyWeinberg implements VariantAnnotation {
private static final int MIN_GENOTYPE_QUALITY = 10;
private static final int MIN_NEG_LOG10_PERROR = MIN_GENOTYPE_QUALITY / 10;
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation) {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, VariantContext vc) {
if ( !(variation instanceof VariantBackedByGenotype) )
return null;
final List<Genotype> genotypes = ((VariantBackedByGenotype)variation).getGenotypes();
final Map<String, Genotype> genotypes = vc.getGenotypes();
if ( genotypes == null || genotypes.size() < MIN_SAMPLES )
return null;
int refCount = 0;
int hetCount = 0;
int homCount = 0;
for ( Genotype genotype : genotypes ) {
if ( genotype.isNoCall() )
for ( Map.Entry<String, Genotype> genotype : genotypes.entrySet() ) {
Genotype g = genotype.getValue();
if ( g.isNoCall() )
continue;
// TODO - fix me:
// Right now we just ignore genotypes that are not confident, but this throws off
// our HW ratios. More analysis is needed to determine the right thing to do when
// the genotyper cannot decide whether a given sample is het or hom var.
if ( genotype.getNegLog10PError() < MIN_NEG_LOG10_PERROR )
if ( g.getNegLog10PError() < MIN_NEG_LOG10_PERROR )
continue;
if ( !genotype.isVariant(ref.getBase()) )
if ( g.isHomRef() )
refCount++;
else if ( genotype.isHet() )
else if ( g.isHet() )
hetCount++;
else
homCount++;

8
java/src/org/broadinstitute/sting/gatk/walkers/annotator/HomopolymerRun.java
View File

@ -2,9 +2,9 @@ package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import java.util.Map;
@ -12,12 +12,12 @@ import java.util.Map;
public class HomopolymerRun extends StandardVariantAnnotation {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation) {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, VariantContext vc) {
if ( !variation.isBiallelic() || !variation.isSNP() )
if ( !vc.isBiallelic() || !vc.isSNP() )
return null;
int run = computeHomopolymerRun(variation.getAlternativeBaseForSNP(), ref);
int run = computeHomopolymerRun(vc.getAlternateAllele(0).toString().charAt(0), ref);
return String.format("%d", run);
}

4
java/src/org/broadinstitute/sting/gatk/walkers/annotator/LowMQ.java
View File

@ -2,10 +2,10 @@ package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import java.util.Map;
@ -13,7 +13,7 @@ import java.util.Map;
public class LowMQ implements VariantAnnotation {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation) {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, VariantContext vc) {
double mq0 = 0;
double mq10 = 0;
double total = 0;

4
java/src/org/broadinstitute/sting/gatk/walkers/annotator/MappingQualityZero.java
View File

@ -2,10 +2,10 @@ package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import java.util.Map;
@ -13,7 +13,7 @@ import java.util.Map;
public class MappingQualityZero extends StandardVariantAnnotation {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation) {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, VariantContext vc) {
int mq0 = 0;
for ( String sample : stratifiedContexts.keySet() ) {
ReadBackedPileup pileup = stratifiedContexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).getBasePileup();

42
java/src/org/broadinstitute/sting/gatk/walkers/annotator/QualByDepth.java
View File

@ -2,33 +2,27 @@ package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.*;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.VariantBackedByGenotype;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.genotype.SampleBacked;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import java.util.Map;
import java.util.List;
import java.util.ArrayList;
public class QualByDepth extends StandardVariantAnnotation {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation) {
if ( !(variation instanceof VariantBackedByGenotype) )
return null;
final List<Genotype> genotypes = ((VariantBackedByGenotype)variation).getGenotypes();
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, VariantContext vc) {
final Map<String, Genotype> genotypes = vc.getGenotypes();
if ( genotypes == null || genotypes.size() == 0 )
return null;
//double QbyD = genotypeQualByDepth(ref.getBase(), genotypes, stratifiedContexts);
int qDepth = variationQualByDepth(ref.getBase(), genotypes, stratifiedContexts);
//double QbyD = genotypeQualByDepth(genotypes, stratifiedContexts);
int qDepth = variationQualByDepth(genotypes, stratifiedContexts);
if ( qDepth == 0 )
return null;
double QbyD = 10.0 * variation.getNegLog10PError() / (double)qDepth;
double QbyD = 10.0 * vc.getNegLog10PError() / (double)qDepth;
return String.format("%.2f", QbyD);
}
@ -36,18 +30,15 @@ public class QualByDepth extends StandardVariantAnnotation {
public VCFInfoHeaderLine getDescription() { return new VCFInfoHeaderLine(getKeyName(), 1, VCFInfoHeaderLine.INFO_TYPE.Float, "Variant Confidence/Quality by Depth"); }
private int variationQualByDepth(char ref, final List<Genotype> genotypes, Map<String, StratifiedAlignmentContext> stratifiedContexts) {
private int variationQualByDepth(final Map<String, Genotype> genotypes, Map<String, StratifiedAlignmentContext> stratifiedContexts) {
int depth = 0;
for ( Genotype genotype : genotypes ) {
if ( !(genotype instanceof SampleBacked) )
continue;
for ( Map.Entry<String, Genotype> genotype : genotypes.entrySet() ) {
// we care only about variant calls
if ( !genotype.isVariant(ref) )
if ( genotype.getValue().isHomRef() )
continue;
String sample = ((SampleBacked)genotype).getSampleName();
StratifiedAlignmentContext context = stratifiedContexts.get(sample);
StratifiedAlignmentContext context = stratifiedContexts.get(genotype.getKey());
if ( context != null )
depth += context.getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).size();
}
@ -55,22 +46,19 @@ public class QualByDepth extends StandardVariantAnnotation {
return depth;
}
private double genotypeQualByDepth(char ref, final List<Genotype> genotypes, Map<String, StratifiedAlignmentContext> stratifiedContexts) {
private double genotypeQualByDepth(final Map<String, Genotype> genotypes, Map<String, StratifiedAlignmentContext> stratifiedContexts) {
ArrayList<Double> qualsByDepth = new ArrayList<Double>();
for ( Genotype genotype : genotypes ) {
if ( !(genotype instanceof SampleBacked) )
continue;
for ( Map.Entry<String, Genotype> genotype : genotypes.entrySet() ) {
// we care only about variant calls
if ( !genotype.isVariant(ref) )
if ( genotype.getValue().isHomRef() )
continue;
String sample = ((SampleBacked)genotype).getSampleName();
StratifiedAlignmentContext context = stratifiedContexts.get(sample);
StratifiedAlignmentContext context = stratifiedContexts.get(genotype.getKey());
if ( context == null )
continue;
qualsByDepth.add(genotype.getNegLog10PError() / context.getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).size());
qualsByDepth.add(genotype.getValue().getNegLog10PError() / context.getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).size());
}
double mean = 0.0;

10
java/src/org/broadinstitute/sting/gatk/walkers/annotator/QualityAdjustedSecondBaseLod.java
View File

@ -2,8 +2,8 @@ package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import java.util.Map;
@ -26,12 +26,12 @@ public class QualityAdjustedSecondBaseLod implements VariantAnnotation {
public VCFInfoHeaderLine getDescription() { return new VCFInfoHeaderLine(KEY_NAME, 1, VCFInfoHeaderLine.INFO_TYPE.Float, "Adjusted residual quality based on second-base skew"); }
public String annotate( RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> contexts, Variation variant) {
String chi = skewCalc.annotate(tracker, ref,contexts,variant);
public String annotate( RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> contexts, VariantContext vc) {
String chi = skewCalc.annotate(tracker, ref,contexts,vc);
if ( chi == null )
return null;
double chi_square = Double.valueOf(chi);
double chi_square = Double.valueOf(chi);
double chi_loglik = chi_square <= 0.0 ? 0.0 : Math.max(-(chi_square/2.0)*log10e + log10half,CHI_LOD_MAX); // cap it...
return String.format("%f", 10*(variant.getNegLog10PError() + chi_loglik));
return String.format("%f", 10*(vc.getNegLog10PError() + chi_loglik));
}
}

4
java/src/org/broadinstitute/sting/gatk/walkers/annotator/RMSMappingQuality.java
View File

@ -2,11 +2,11 @@ package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.vcf.VCFRecord;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
@ -16,7 +16,7 @@ import java.util.ArrayList;
public class RMSMappingQuality extends StandardVariantAnnotation {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation) {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, VariantContext vc) {
ArrayList<Integer> qualities = new ArrayList<Integer>();
for ( String sample : stratifiedContexts.keySet() ) {
ReadBackedPileup pileup = stratifiedContexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).getBasePileup();

17
java/src/org/broadinstitute/sting/gatk/walkers/annotator/RankSumTest.java
View File

@ -2,10 +2,10 @@ package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.*;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.genotype.*;
import java.util.List;
import java.util.ArrayList;
@ -16,27 +16,26 @@ public abstract class RankSumTest implements VariantAnnotation {
private final static boolean DEBUG = false;
private static final double minPValue = 1e-10;
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation) {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, VariantContext vc) {
if ( !variation.isBiallelic() || !variation.isSNP() || !(variation instanceof VariantBackedByGenotype) )
if ( !vc.isBiallelic() || !vc.isSNP() )
return null;
final List<Genotype> genotypes = ((VariantBackedByGenotype)variation).getGenotypes();
final Map<String, Genotype> genotypes = vc.getGenotypes();
if ( genotypes == null || genotypes.size() == 0 )
return null;
ArrayList<Integer> refQuals = new ArrayList<Integer>();
ArrayList<Integer> altQuals = new ArrayList<Integer>();
for ( Genotype genotype : genotypes ) {
for ( Map.Entry<String, Genotype> genotype : genotypes.entrySet() ) {
// we care only about het calls
if ( genotype.isHet() ) {
String sample = ((SampleBacked)genotype).getSampleName();
StratifiedAlignmentContext context = stratifiedContexts.get(sample);
if ( genotype.getValue().isHet() ) {
StratifiedAlignmentContext context = stratifiedContexts.get(genotype.getKey());
if ( context == null )
continue;
fillQualsFromPileup(ref.getBase(), variation.getAlternativeBaseForSNP(), context.getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).getBasePileup(), refQuals, altQuals);
fillQualsFromPileup(ref.getBase(), vc.getAlternateAllele(0).toString().charAt(0), context.getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).getBasePileup(), refQuals, altQuals);
}
}

8
java/src/org/broadinstitute/sting/gatk/walkers/annotator/SecondBaseSkew.java
View File

@ -4,10 +4,10 @@ import org.broadinstitute.sting.utils.Pair;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.BaseUtils;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import java.util.Map;
@ -30,11 +30,11 @@ public class SecondBaseSkew implements VariantAnnotation {
public VCFInfoHeaderLine getDescription() { return new VCFInfoHeaderLine(KEY_NAME, 1, VCFInfoHeaderLine.INFO_TYPE.Float, "Chi-square Secondary Base Skew"); }
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation) {
if ( !variation.isBiallelic() || !variation.isSNP() )
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, VariantContext vc) {
if ( !vc.isBiallelic() || !vc.isSNP() )
return null;
char alternate = variation.getAlternativeBaseForSNP();
char alternate = vc.getAlternateAllele(0).toString().charAt(0);
Pair<Integer, Integer> depth = new Pair<Integer, Integer>(0, 0);
for ( String sample : stratifiedContexts.keySet() ) {

4
java/src/org/broadinstitute/sting/gatk/walkers/annotator/SpanningDeletions.java
View File

@ -2,9 +2,9 @@ package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import java.util.Map;
@ -12,7 +12,7 @@ import java.util.Map;
public class SpanningDeletions extends StandardVariantAnnotation {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation) {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, VariantContext vc) {
int deletions = 0;
int depth = 0;
for ( String sample : stratifiedContexts.keySet() ) {

4
java/src/org/broadinstitute/sting/gatk/walkers/annotator/VariantAnnotation.java
View File

@ -2,8 +2,8 @@ package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import java.util.Map;
@ -12,7 +12,7 @@ import java.util.Map;
public interface VariantAnnotation {
// return the annotation for the given variation and context split by sample (return null for no annotation)
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation);
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, VariantContext vc);
// return the INFO key
public String getKeyName();

64
java/src/org/broadinstitute/sting/gatk/walkers/annotator/VariantAnnotator.java
View File

@ -1,16 +1,15 @@
package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.*;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.*;
import org.broadinstitute.sting.gatk.refdata.utils.RODRecordList;
import org.broadinstitute.sting.gatk.walkers.*;
import org.broadinstitute.sting.gatk.datasources.simpleDataSources.ReferenceOrderedDataSource;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.genotype.vcf.*;
import org.broadinstitute.sting.utils.genotype.*;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.cmdLine.Argument;
import org.broadinstitute.sting.playground.gatk.walkers.variantstovcf.VariantsToVCF;
import java.util.*;
import java.io.*;
@ -38,7 +37,6 @@ public class VariantAnnotator extends LocusWalker<Integer, Integer> {
protected Boolean LIST = false;
private VCFWriter vcfWriter;
private VCFHeader vcfHeader;
private HashMap<String, String> nonVCFsampleName = new HashMap<String, String>();
@ -155,7 +153,7 @@ public class VariantAnnotator extends LocusWalker<Integer, Integer> {
hInfo.add(new VCFInfoHeaderLine(VCFRecord.HAPMAP3_KEY,1,VCFInfoHeaderLine.INFO_TYPE.Integer, "Hapmap 3 membership"));
vcfWriter = new VCFWriter(VCF_OUT);
vcfHeader = new VCFHeader(hInfo, samples);
VCFHeader vcfHeader = new VCFHeader(hInfo, samples);
vcfWriter.writeHeader(vcfHeader);
}
@ -192,15 +190,26 @@ public class VariantAnnotator extends LocusWalker<Integer, Integer> {
return 0;
Map<String, String> annotations = new HashMap<String, String>();
VariationRod variant = (VariationRod)rods.get(0);
ReferenceOrderedDatum variant = rods.get(0);
VariantContext vc = VariantContextAdaptors.toVariantContext("variant", variant);
if ( vc == null )
return 0;
// if the reference base is not ambiguous, we can annotate
if ( BaseUtils.simpleBaseToBaseIndex(ref.getBase()) != -1 ) {
Map<String, StratifiedAlignmentContext> stratifiedContexts = StratifiedAlignmentContext.splitContextBySample(context.getBasePileup());
if ( stratifiedContexts != null )
annotations = getAnnotations(tracker, ref, stratifiedContexts, variant, requestedAnnotations, annotateDbsnp, annotateHapmap2, annotateHapmap3);
annotations = getAnnotations(tracker, ref, stratifiedContexts, vc, requestedAnnotations, annotateDbsnp, annotateHapmap2, annotateHapmap3);
}
writeVCF(tracker, ref, context, variant, annotations);
VCFRecord record;
if ( variant instanceof RodVCF )
record = ((RodVCF)variant).mCurrentRecord;
else
record = VariantContextAdaptors.toVCF(vc, ref.getBase());
record.addInfoFields(annotations);
writeVCF(tracker, record);
return 1;
}
@ -240,21 +249,21 @@ public class VariantAnnotator extends LocusWalker<Integer, Integer> {
}
// option #1: don't specify annotations to be used: standard annotations are used by default
public static Map<String, String> getAnnotations(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation, boolean annotateDbsnp, boolean annotateHapmap2, boolean annotateHapmap3) {
public static Map<String, String> getAnnotations(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, VariantContext vc, boolean annotateDbsnp, boolean annotateHapmap2, boolean annotateHapmap3) {
if ( standardAnnotations == null )
determineAllAnnotations();
return getAnnotations(tracker, ref, stratifiedContexts, variation, standardAnnotations.values(), annotateDbsnp, annotateHapmap2, annotateHapmap3);
return getAnnotations(tracker, ref, stratifiedContexts, vc, standardAnnotations.values(), annotateDbsnp, annotateHapmap2, annotateHapmap3);
}
// option #2: specify that all possible annotations be used
public static Map<String, String> getAllAnnotations(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation, boolean annotateDbsnp, boolean annotateHapmap2, boolean annotateHapmap3) {
public static Map<String, String> getAllAnnotations(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, VariantContext vc, boolean annotateDbsnp, boolean annotateHapmap2, boolean annotateHapmap3) {
if ( allAnnotations == null )
determineAllAnnotations();
return getAnnotations(tracker, ref, stratifiedContexts, variation, allAnnotations.values(), annotateDbsnp, annotateHapmap2, annotateHapmap3);
return getAnnotations(tracker, ref, stratifiedContexts, vc, allAnnotations.values(), annotateDbsnp, annotateHapmap2, annotateHapmap3);
}
// option #3: specify the exact annotations to be used
public static Map<String, String> getAnnotations(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation, Collection<VariantAnnotation> annotations, boolean annotateDbsnp, boolean annotateHapmap2, boolean annotateHapmap3) {
public static Map<String, String> getAnnotations(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, VariantContext vc, Collection<VariantAnnotation> annotations, boolean annotateDbsnp, boolean annotateHapmap2, boolean annotateHapmap3) {
HashMap<String, String> results = new HashMap<String, String>();
@ -275,7 +284,7 @@ public class VariantAnnotator extends LocusWalker<Integer, Integer> {
}
for ( VariantAnnotation annotator : annotations) {
String annot = annotator.annotate(tracker, ref, stratifiedContexts, variation);
String annot = annotator.annotate(tracker, ref, stratifiedContexts, vc);
if ( annot != null ) {
results.put(annotator.getKeyName(), annot);
}
@ -284,29 +293,14 @@ public class VariantAnnotator extends LocusWalker<Integer, Integer> {
return results;
}
private void writeVCF(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context, VariationRod variant, Map<String, String> annotations) {
VCFRecord rec = getVCFRecord(tracker, ref, context, variant);
if ( rec != null ) {
rec.addInfoFields(annotations);
// also, annotate dbsnp id if available and not already there
if ( annotateDbsnp && (rec.getID() == null || rec.getID().equals(".")) ) {
rodDbSNP dbsnp = rodDbSNP.getFirstRealSNP(tracker.getTrackData("dbsnp", null));
if ( dbsnp != null )
rec.setID(dbsnp.getRS_ID());
}
vcfWriter.addRecord(rec);
}
}
private VCFRecord getVCFRecord(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context, VariationRod variant) {
if ( variant instanceof RodVCF ) {
return ((RodVCF)variant).mCurrentRecord;
} else {
List<VCFGenotypeRecord> gt = new ArrayList<VCFGenotypeRecord>();
Map<VCFHeader.HEADER_FIELDS, String> map = new HashMap<VCFHeader.HEADER_FIELDS, String>();
return VariantsToVCF.generateVCFRecord(tracker, ref, context, vcfHeader, gt, map, nonVCFsampleName, out, false, false);
private void writeVCF(RefMetaDataTracker tracker, VCFRecord record) {
// annotate dbsnp id if available and not already there
if ( annotateDbsnp && (record.getID() == null || record.getID().equals(VCFRecord.EMPTY_ID_FIELD)) ) {
rodDbSNP dbsnp = rodDbSNP.getFirstRealSNP(tracker.getTrackData("dbsnp", null));
if ( dbsnp != null )
record.setID(dbsnp.getRS_ID());
}
vcfWriter.addRecord(record);
}
/**

3
java/src/org/broadinstitute/sting/gatk/walkers/concordance/CallsetConcordanceWalker.java
View File

@ -4,7 +4,6 @@ import org.broadinstitute.sting.gatk.contexts.*;
import org.broadinstitute.sting.gatk.refdata.*;
import org.broadinstitute.sting.gatk.walkers.*;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.genotype.vcf.*;
import org.broadinstitute.sting.utils.cmdLine.Argument;
@ -140,7 +139,7 @@ public class CallsetConcordanceWalker extends RodWalker<Integer, Integer> {
// pull out all of the individual calls from the rods and insert into a map based on the
// mapping from rod/sample to uniquified name
HashMap<String, Genotype> samplesToRecords = new HashMap<String, Genotype>();
HashMap<String, VCFGenotypeRecord> samplesToRecords = new HashMap<String, VCFGenotypeRecord>();
for ( RodVCF rod : vcfRods ) {
List<VCFGenotypeRecord> records = rod.getVCFGenotypeRecords();
for ( VCFGenotypeRecord vcfRec : records ) {

4
java/src/org/broadinstitute/sting/gatk/walkers/concordance/ConcordanceType.java
View File

@ -1,8 +1,8 @@
package org.broadinstitute.sting.gatk.walkers.concordance;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import org.broadinstitute.sting.utils.genotype.vcf.VCFGenotypeRecord;
import java.util.Map;
import java.util.Set;
@ -10,7 +10,7 @@ import java.util.Set;
public interface ConcordanceType {
public void initialize(Map<String,String> args, Set<String> samples);
public String computeConcordance(Map<String, Genotype> samplesToRecords, ReferenceContext ref);
public String computeConcordance(Map<String, VCFGenotypeRecord> samplesToRecords, ReferenceContext ref);
public String getInfoName();
public VCFInfoHeaderLine getInfoDescription();
}

14
java/src/org/broadinstitute/sting/gatk/walkers/concordance/IndelSubsets.java
View File

@ -3,9 +3,9 @@ package org.broadinstitute.sting.gatk.walkers.concordance;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.StingException;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import org.broadinstitute.sting.utils.genotype.vcf.VCFGenotypeRecord;
import org.broadinstitute.sting.utils.genotype.vcf.VCFRecord;
import java.util.*;
@ -54,10 +54,10 @@ public class IndelSubsets implements ConcordanceType {
}
}
public String computeConcordance(Map<String, Genotype> samplesToRecords, ReferenceContext ref) {
public String computeConcordance(Map<String, VCFGenotypeRecord> samplesToRecords, ReferenceContext ref) {
Genotype indel1 = samplesToRecords.get(sample1);
Genotype indel2 = samplesToRecords.get(sample2);
VCFGenotypeRecord indel1 = samplesToRecords.get(sample1);
VCFGenotypeRecord indel2 = samplesToRecords.get(sample2);
int set1 = ( indel1 != null && !indel1.isPointGenotype() ? 0 : 1 );
int set2 = ( indel2 != null && !indel2.isPointGenotype() ? 0 : 1 );
@ -67,7 +67,7 @@ public class IndelSubsets implements ConcordanceType {
return null;
// only deal with a valid indel
Variation indel = ( indel1 != null ? indel1.toVariation(ref.getBase()) : indel2.toVariation(ref.getBase()) );
VCFRecord indel = ( indel1 != null ? indel1.getRecord() : indel2.getRecord() );
// we only deal with the first allele
int size = ( indel.getAlternateAlleleList().get(0).length() <= sizeCutoff ? 0 : 1 );
@ -76,7 +76,7 @@ public class IndelSubsets implements ConcordanceType {
return tags[set1][set2][size][homopol];
}
private int homopolymerRunSize(ReferenceContext ref, Variation indel) {
private int homopolymerRunSize(ReferenceContext ref, VCFRecord indel) {
char[] bases = ref.getBases();
GenomeLoc window = ref.getWindow();
GenomeLoc locus = ref.getLocus();

6
java/src/org/broadinstitute/sting/gatk/walkers/concordance/NWayVenn.java
View File

@ -1,8 +1,8 @@
package org.broadinstitute.sting.gatk.walkers.concordance;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import org.broadinstitute.sting.utils.genotype.vcf.VCFGenotypeRecord;
import java.util.*;
import java.util.Map.Entry;
@ -19,12 +19,12 @@ public class NWayVenn implements ConcordanceType {
public void initialize(Map<String, String> args, Set<String> samples) { }
public String computeConcordance(Map<String, Genotype> samplesToRecords, ReferenceContext ref) {
public String computeConcordance(Map<String, VCFGenotypeRecord> samplesToRecords, ReferenceContext ref) {
if ( samplesToRecords.size() == 0 )
return null;
TreeSet<String> concordantSamples = new TreeSet<String>();
for ( Entry<String, Genotype> entry : samplesToRecords.entrySet() ) {
for ( Entry<String, VCFGenotypeRecord> entry : samplesToRecords.entrySet() ) {
if ( !entry.getValue().isNoCall() )
concordantSamples.add(entry.getKey());
}

8
java/src/org/broadinstitute/sting/gatk/walkers/concordance/SNPGenotypeConcordance.java
View File

@ -2,8 +2,8 @@ package org.broadinstitute.sting.gatk.walkers.concordance;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.utils.StingException;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import org.broadinstitute.sting.utils.genotype.vcf.VCFGenotypeRecord;
import java.util.*;
@ -30,13 +30,13 @@ public class SNPGenotypeConcordance implements ConcordanceType {
sample2 = iter.next();
}
public String computeConcordance(Map<String, Genotype> samplesToRecords, ReferenceContext ref) {
public String computeConcordance(Map<String, VCFGenotypeRecord> samplesToRecords, ReferenceContext ref) {
char refBase = ref.getBase();
Genotype call1 = samplesToRecords.get(sample1);
VCFGenotypeRecord call1 = samplesToRecords.get(sample1);
if ( call1 != null && call1.isNoCall() )
call1 = null;
Genotype call2 = samplesToRecords.get(sample2);
VCFGenotypeRecord call2 = samplesToRecords.get(sample2);
if ( call2 != null && call2.isNoCall() )
call2 = null;

14
java/src/org/broadinstitute/sting/gatk/walkers/concordance/SimpleVenn.java
View File

@ -1,9 +1,9 @@
package org.broadinstitute.sting.gatk.walkers.concordance;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import org.broadinstitute.sting.utils.genotype.vcf.VCFGenotypeRecord;
import org.broadinstitute.sting.utils.genotype.vcf.VCFRecord;
import org.broadinstitute.sting.utils.StingException;
import java.util.*;
@ -26,12 +26,12 @@ public class SimpleVenn implements ConcordanceType {
sample2 = iter.next();
}
public String computeConcordance(Map<String, Genotype> samplesToRecords, ReferenceContext ref ) {
public String computeConcordance(Map<String, VCFGenotypeRecord> samplesToRecords, ReferenceContext ref ) {
Genotype call1 = samplesToRecords.get(sample1);
VCFGenotypeRecord call1 = samplesToRecords.get(sample1);
if ( call1 != null && call1.isNoCall() )
call1 = null;
Genotype call2 = samplesToRecords.get(sample2);
VCFGenotypeRecord call2 = samplesToRecords.get(sample2);
if ( call2 != null && call2.isNoCall() )
call2 = null;
@ -47,8 +47,8 @@ public class SimpleVenn implements ConcordanceType {
return sample2 + "_only";
// at this point we know that neither is null, so now we need to test for alternate allele concordance
Variation callV1 = call1.toVariation(ref.getBase());
Variation callV2 = call2.toVariation(ref.getBase());
VCFRecord callV1 = call1.getRecord();
VCFRecord callV2 = call2.getRecord();
// we can't really deal with multi-allelic variants
if ( callV1.isBiallelic() && callV2.isBiallelic() ) {

60
java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidGenotypeCalculationModel.java
View File

@ -1,9 +1,12 @@
package org.broadinstitute.sting.gatk.walkers.genotyper;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.genotype.DiploidGenotype;
import org.broadinstitute.sting.utils.genotype.CalledGenotype;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.genotype.*;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.Genotype;
import org.broadinstitute.sting.gatk.contexts.variantcontext.Allele;
import java.util.*;
@ -83,53 +86,44 @@ public class DiploidGenotypeCalculationModel extends JointEstimateGenotypeCalcul
}
}
protected List<Genotype> makeGenotypeCalls(char ref, char alt, int frequency, Map<String, StratifiedAlignmentContext> contexts, GenomeLoc loc) {
ArrayList<Genotype> calls = new ArrayList<Genotype>();
protected Map<String, Genotype> makeGenotypeCalls(char ref, char alt, int frequency, Map<String, StratifiedAlignmentContext> contexts, GenomeLoc loc) {
HashMap<String, Genotype> calls = new HashMap<String, Genotype>();
// set up some variables we'll need in the loop
AlleleFrequencyMatrix matrix = AFMatrixMap.get(alt);
DiploidGenotype refGenotype = DiploidGenotype.createHomGenotype(ref);
DiploidGenotype hetGenotype = DiploidGenotype.createDiploidGenotype(ref, alt);
DiploidGenotype homGenotype = DiploidGenotype.createHomGenotype(alt);
Allele refAllele = new Allele(Character.toString(ref), true);
Allele altAllele = new Allele(Character.toString(alt), false);
for ( String sample : GLs.keySet() ) {
// create the call
GenotypeCall call = GenotypeWriterFactory.createSupportedGenotypeCall(OUTPUT_FORMAT, Character.toString(ref), loc);
// set the genotype and confidence
Pair<Integer, Double> AFbasedGenotype = matrix.getGenotype(frequency, sample);
call.setNegLog10PError(AFbasedGenotype.second);
if ( AFbasedGenotype.first == GenotypeType.REF.ordinal() )
call.setGenotype(refGenotype);
else if ( AFbasedGenotype.first == GenotypeType.HET.ordinal() )
call.setGenotype(hetGenotype);
else // ( AFbasedGenotype.first == GenotypeType.HOM.ordinal() )
call.setGenotype(homGenotype);
if ( call instanceof ReadBacked ) {
ReadBackedPileup pileup = contexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).getBasePileup();
((ReadBacked)call).setPileup(pileup);
}
if ( call instanceof SampleBacked ) {
((SampleBacked)call).setSampleName(sample);
}
if ( call instanceof LikelihoodsBacked ) {
((LikelihoodsBacked)call).setLikelihoods(GLs.get(sample).getLikelihoods());
}
if ( call instanceof PosteriorsBacked ) {
((PosteriorsBacked)call).setPosteriors(GLs.get(sample).getPosteriors());
}
if ( call instanceof AlleleConstrainedGenotype ) {
((AlleleConstrainedGenotype)call).setAlternateAllele(alt);
ArrayList<Allele> myAlleles = new ArrayList<Allele>();
if ( AFbasedGenotype.first == GenotypeType.REF.ordinal() ) {
myAlleles.add(refAllele);
myAlleles.add(refAllele);
} else if ( AFbasedGenotype.first == GenotypeType.HET.ordinal() ) {
myAlleles.add(refAllele);
myAlleles.add(altAllele);
} else { // ( AFbasedGenotype.first == GenotypeType.HOM.ordinal() )
myAlleles.add(altAllele);
myAlleles.add(altAllele);
}
calls.add(call);
CalledGenotype cg = new CalledGenotype(sample, myAlleles, AFbasedGenotype.second);
cg.setLikelihoods(GLs.get(sample).getLikelihoods());
cg.setPosteriors(GLs.get(sample).getPosteriors());
cg.setReadBackedPileup(contexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).getBasePileup());
calls.put(sample, cg);
}
// output to beagle file if requested
if ( beagleWriter != null ) {
DiploidGenotype refGenotype = DiploidGenotype.createHomGenotype(ref);
DiploidGenotype hetGenotype = DiploidGenotype.createDiploidGenotype(ref, alt);
DiploidGenotype homGenotype = DiploidGenotype.createHomGenotype(alt);
for ( String sample : samples ) {
GenotypeLikelihoods gl = GLs.get(sample);
if ( gl == null ) {

106
java/src/org/broadinstitute/sting/gatk/walkers/genotyper/JointEstimateGenotypeCalculationModel.java
View File

@ -2,11 +2,11 @@ package org.broadinstitute.sting.gatk.walkers.genotyper;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.pileup.*;
import org.broadinstitute.sting.utils.genotype.*;
import org.broadinstitute.sting.utils.genotype.Variation.VARIANT_TYPE;
import org.broadinstitute.sting.utils.genotype.vcf.VCFRecord;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.refdata.rodDbSNP;
import org.broadinstitute.sting.gatk.contexts.*;
import org.broadinstitute.sting.gatk.contexts.variantcontext.*;
import java.util.*;
@ -73,7 +73,7 @@ public abstract class JointEstimateGenotypeCalculationModel extends GenotypeCalc
// technically, at this point our confidence in a reference call isn't accurately
// estimated because it didn't take into account samples with no data
if ( vcc.variation == null )
if ( vcc.vc == null )
estimateReferenceConfidence(vcc, contexts, DiploidGenotypePriors.HUMAN_HETEROZYGOSITY, true);
return vcc;
}
@ -306,9 +306,9 @@ public abstract class JointEstimateGenotypeCalculationModel extends GenotypeCalc
verboseWriter.println();
}
protected List<Genotype> makeGenotypeCalls(char ref, char alt, int frequency, Map<String, StratifiedAlignmentContext> contexts, GenomeLoc loc) {
protected Map<String, Genotype> makeGenotypeCalls(char ref, char alt, int frequency, Map<String, StratifiedAlignmentContext> contexts, GenomeLoc loc) {
// by default, we return no genotypes
return new ArrayList<Genotype>();
return new HashMap<String, Genotype>();
}
protected VariantCallContext createCalls(RefMetaDataTracker tracker, char ref, Map<String, StratifiedAlignmentContext> contexts, GenomeLoc loc, int frequencyEstimationPoints) {
@ -345,70 +345,62 @@ public abstract class JointEstimateGenotypeCalculationModel extends GenotypeCalc
}
// populate the sample-specific data (output it to beagle also if requested)
List<Genotype> calls = makeGenotypeCalls(ref, bestAlternateAllele, bestAFguess, contexts, loc);
Map<String, Genotype> genotypes = makeGenotypeCalls(ref, bestAlternateAllele, bestAFguess, contexts, loc);
// close beagle record (if requested)
if ( beagleWriter != null ) {
if ( beagleWriter != null )
beagleWriter.println();
}
// next, the general locus data
// next, the variant context data (alleles, attributes, etc.)
ArrayList<Allele> alleles = new ArrayList<Allele>();
alleles.add(new Allele(Character.toString(ref), true));
if ( bestAFguess != 0 )
alleles.add(new Allele(bestAlternateAllele.toString(), false));
// *** note that calculating strand bias involves overwriting data structures, so we do that last
VariationCall locusdata = GenotypeWriterFactory.createSupportedCall(OUTPUT_FORMAT, ref, loc, bestAFguess == 0 ? VARIANT_TYPE.REFERENCE : VARIANT_TYPE.SNP);
if ( locusdata != null ) {
if ( bestAFguess != 0 ) {
locusdata.addAlternateAllele(bestAlternateAllele.toString());
locusdata.setNonRefAlleleFrequency((double)bestAFguess / (double)(frequencyEstimationPoints-1));
}
if ( locusdata instanceof ConfidenceBacked ) {
((ConfidenceBacked)locusdata).setConfidence(phredScaledConfidence);
}
if ( locusdata instanceof IDBacked ) {
rodDbSNP dbsnp = getDbSNP(tracker);
if ( dbsnp != null )
((IDBacked)locusdata).setID(dbsnp.getRS_ID());
}
if ( locusdata instanceof SLODBacked && REPORT_SLOD ) {
// the overall lod
double overallLog10PofNull = log10AlleleFrequencyPriors[0] + log10PofDgivenAFi[indexOfMax][0];
double overallLog10PofF = log10AlleleFrequencyPriors[bestAFguess] + log10PofDgivenAFi[indexOfMax][bestAFguess];
double lod = overallLog10PofF - overallLog10PofNull;
HashMap<String, Object> attributes = new HashMap<String, Object>();
if ( bestAFguess != 0 )
attributes.put(VCFRecord.ALLELE_FREQUENCY_KEY, new Double((double)bestAFguess / (double)(frequencyEstimationPoints-1)));
// the forward lod
initialize(ref, contexts, StratifiedAlignmentContext.StratifiedContextType.FORWARD);
calculateAlleleFrequencyPosteriors(ref, frequencyEstimationPoints, contexts, StratifiedAlignmentContext.StratifiedContextType.FORWARD);
calculatePofFs(ref, frequencyEstimationPoints);
double forwardLog10PofNull = log10AlleleFrequencyPriors[0] + log10PofDgivenAFi[indexOfMax][0];
double forwardLog10PofF = log10AlleleFrequencyPriors[bestAFguess] + log10PofDgivenAFi[indexOfMax][bestAFguess];
rodDbSNP dbsnp = getDbSNP(tracker);
if ( dbsnp != null )
attributes.put("ID", dbsnp.getRS_ID());
// the reverse lod
initialize(ref, contexts, StratifiedAlignmentContext.StratifiedContextType.REVERSE);
calculateAlleleFrequencyPosteriors(ref, frequencyEstimationPoints, contexts, StratifiedAlignmentContext.StratifiedContextType.REVERSE);
calculatePofFs(ref, frequencyEstimationPoints);
double reverseLog10PofNull = log10AlleleFrequencyPriors[0] + log10PofDgivenAFi[indexOfMax][0];
double reverseLog10PofF = log10AlleleFrequencyPriors[bestAFguess] + log10PofDgivenAFi[indexOfMax][bestAFguess];
if ( REPORT_SLOD ) {
// the overall lod
double overallLog10PofNull = log10AlleleFrequencyPriors[0] + log10PofDgivenAFi[indexOfMax][0];
double overallLog10PofF = log10AlleleFrequencyPriors[bestAFguess] + log10PofDgivenAFi[indexOfMax][bestAFguess];
double lod = overallLog10PofF - overallLog10PofNull;
double forwardLod = forwardLog10PofF + reverseLog10PofNull - overallLog10PofNull;
double reverseLod = reverseLog10PofF + forwardLog10PofNull - overallLog10PofNull;
//logger.debug("forward lod=" + forwardLod + ", reverse lod=" + reverseLod);
// the forward lod
initialize(ref, contexts, StratifiedAlignmentContext.StratifiedContextType.FORWARD);
calculateAlleleFrequencyPosteriors(ref, frequencyEstimationPoints, contexts, StratifiedAlignmentContext.StratifiedContextType.FORWARD);
calculatePofFs(ref, frequencyEstimationPoints);
double forwardLog10PofNull = log10AlleleFrequencyPriors[0] + log10PofDgivenAFi[indexOfMax][0];
double forwardLog10PofF = log10AlleleFrequencyPriors[bestAFguess] + log10PofDgivenAFi[indexOfMax][bestAFguess];
// strand score is max bias between forward and reverse strands
double strandScore = Math.max(forwardLod - lod, reverseLod - lod);
// rescale by a factor of 10
strandScore *= 10.0;
//logger.debug(String.format("SLOD=%f", strandScore));
// the reverse lod
initialize(ref, contexts, StratifiedAlignmentContext.StratifiedContextType.REVERSE);
calculateAlleleFrequencyPosteriors(ref, frequencyEstimationPoints, contexts, StratifiedAlignmentContext.StratifiedContextType.REVERSE);
calculatePofFs(ref, frequencyEstimationPoints);
double reverseLog10PofNull = log10AlleleFrequencyPriors[0] + log10PofDgivenAFi[indexOfMax][0];
double reverseLog10PofF = log10AlleleFrequencyPriors[bestAFguess] + log10PofDgivenAFi[indexOfMax][bestAFguess];
((SLODBacked)locusdata).setSLOD(strandScore);
}
double forwardLod = forwardLog10PofF + reverseLog10PofNull - overallLog10PofNull;
double reverseLod = reverseLog10PofF + forwardLog10PofNull - overallLog10PofNull;
//logger.debug("forward lod=" + forwardLod + ", reverse lod=" + reverseLod);
// strand score is max bias between forward and reverse strands
double strandScore = Math.max(forwardLod - lod, reverseLod - lod);
// rescale by a factor of 10
strandScore *= 10.0;
//logger.debug(String.format("SLOD=%f", strandScore));
attributes.put("SB", new Double(strandScore));
}
// finally, associate the Variation with the Genotypes (if available)
if ( locusdata != null ) {
locusdata.setGenotypeCalls(calls);
for ( Genotype call : calls )
((GenotypeCall)call).setVariation(locusdata);
}
VariantContext vc = new MutableVariantContext("UG_SNP_call", loc, alleles, genotypes, phredScaledConfidence/10.0, null, attributes);
return new VariantCallContext(locusdata, calls, phredScaledConfidence >= CONFIDENCE_THRESHOLD);
return new VariantCallContext(vc, phredScaledConfidence >= CONFIDENCE_THRESHOLD);
}
}

10
java/src/org/broadinstitute/sting/gatk/walkers/genotyper/PooledCalculationModel.java
View File

@ -23,9 +23,9 @@ public class PooledCalculationModel extends JointEstimateGenotypeCalculationMode
/**
*
* @param ref
* @param contexts
* @param contextType
* @param ref reference base
* @param contexts alignment contexts
* @param contextType context type
*/
protected void initialize(char ref, Map<String, StratifiedAlignmentContext> contexts, StratifiedAlignmentContext.StratifiedContextType contextType) {
super.initialize(ref, contexts, contextType);
@ -48,7 +48,7 @@ public class PooledCalculationModel extends JointEstimateGenotypeCalculationMode
protected void calculatelog10PofDgivenAFforAllF(char ref, char alt, int nChromosomes, Map<String, StratifiedAlignmentContext> contexts, StratifiedAlignmentContext.StratifiedContextType contextType) {
StratifiedAlignmentContext context = contexts.get(POOL_SAMPLE_NAME);
ReadBackedPileup pileup = context.getContext(contextType).getPileup();
ReadBackedPileup pileup = context.getContext(contextType).getBasePileup();
int refIndex = BaseUtils.simpleBaseToBaseIndex(ref);
int altIndex = BaseUtils.simpleBaseToBaseIndex(alt);
@ -119,7 +119,7 @@ public class PooledCalculationModel extends JointEstimateGenotypeCalculationMode
}
private double calcPBGivenH(int refIndex, int altIndex, int nAltAlleles, int nChromosomes, char base, byte qual, SAMRecord read, int offset) {
double L = 0.0;
double L;
if ( USE_CACHE ) {
L = getCache(refIndex, altIndex, nAltAlleles, base, qual, read);

20
java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java
View File

@ -193,23 +193,13 @@ public class UnifiedGenotyper extends LocusWalker<VariantCallContext, UnifiedGen
sum.nBasesCalledConfidently += value.confidentlyCalled ? 1 : 0;
// can't make a confident variant call here
if ( value.genotypes == null ||
(UAC.genotypeModel != GenotypeCalculationModel.Model.POOLED && value.genotypes.size() == 0) ) {
if ( value.vc == null )
return sum;
}
// if we have a single-sample call (single sample from PointEstimate model returns no VariationCall data)
if ( value.variation == null || (!writer.supportsMultiSample() && UG_engine.samples.size() <= 1) ) {
writer.addGenotypeCall(value.genotypes.get(0));
}
// use multi-sample mode if we have multiple samples or the output type allows it
else {
try {
writer.addMultiSampleCall(value.genotypes, value.variation);
} catch (IllegalArgumentException e) {
throw new IllegalArgumentException(e.getMessage() + "; this is often caused by using the --assume_single_sample_reads argument with the wrong sample name");
}
try {
writer.addCall(value.vc);
} catch (IllegalArgumentException e) {
throw new IllegalArgumentException(e.getMessage() + "; this is often caused by using the --assume_single_sample_reads argument with the wrong sample name");
}
return sum;

13
java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java
View File

@ -28,6 +28,7 @@ package org.broadinstitute.sting.gatk.walkers.genotyper;
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.gatk.datasources.simpleDataSources.ReferenceOrderedDataSource;
import org.broadinstitute.sting.gatk.contexts.*;
import org.broadinstitute.sting.gatk.contexts.variantcontext.MutableVariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.refdata.ReferenceOrderedData;
import org.broadinstitute.sting.gatk.refdata.rodDbSNP;
@ -108,10 +109,6 @@ public class UnifiedGenotyperEngine {
// in pooled mode we need to check that the format is acceptable
if ( UAC.genotypeModel == GenotypeCalculationModel.Model.POOLED && genotypeWriter != null ) {
// only multi-sample calls use Variations
if ( !genotypeWriter.supportsMultiSample() )
throw new IllegalArgumentException("The POOLED model is not compatible with the specified format; try using VCF instead");
// when using VCF with multiple threads, we need to turn down the validation stringency so that writing temporary files will work
if ( toolkit.getArguments().numberOfThreads > 1 && genotypeWriter instanceof VCFGenotypeWriter )
((VCFGenotypeWriter)genotypeWriter).setValidationStringency(VCFGenotypeWriterAdapter.VALIDATION_STRINGENCY.SILENT);
@ -194,16 +191,16 @@ public class UnifiedGenotyperEngine {
VariantCallContext call = gcm.get().callLocus(tracker, ref, rawContext.getLocation(), stratifiedContexts, priors);
// annotate the call, if possible
if ( call != null && call.variation != null && call.variation instanceof ArbitraryFieldsBacked ) {
if ( call != null && call.vc != null ) {
// first off, we want to use the *unfiltered* context for the annotations
stratifiedContexts = StratifiedAlignmentContext.splitContextBySample(rawContext.getBasePileup());
Map<String, String> annotations;
if ( UAC.ALL_ANNOTATIONS )
annotations = VariantAnnotator.getAllAnnotations(tracker, refContext, stratifiedContexts, call.variation, annotateDbsnp, annotateHapmap2, annotateHapmap3);
annotations = VariantAnnotator.getAllAnnotations(tracker, refContext, stratifiedContexts, call.vc, annotateDbsnp, annotateHapmap2, annotateHapmap3);
else
annotations = VariantAnnotator.getAnnotations(tracker, refContext, stratifiedContexts, call.variation, annotateDbsnp, annotateHapmap2, annotateHapmap3);
((ArbitraryFieldsBacked)call.variation).setFields(annotations);
annotations = VariantAnnotator.getAnnotations(tracker, refContext, stratifiedContexts, call.vc, annotateDbsnp, annotateHapmap2, annotateHapmap3);
((MutableVariantContext)call.vc).putAttributes(annotations);
}
return call;

18
java/src/org/broadinstitute/sting/gatk/walkers/genotyper/VariantCallContext.java
View File

@ -1,31 +1,27 @@
package org.broadinstitute.sting.gatk.walkers.genotyper;
import org.broadinstitute.sting.utils.genotype.VariationCall;
import org.broadinstitute.sting.utils.genotype.Genotype;
import java.util.List;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
/**
* Created by IntelliJ IDEA.
* User: depristo
* User: depristo, ebanks
* Date: Jan 22, 2010
* Time: 2:25:19 PM
*
* Useful helper class to communicate the results of calculateGenotype to framework
*/
public class VariantCallContext {
public VariationCall variation = null;
public List<Genotype> genotypes = null;
public VariantContext vc = null;
/** Was the site called confidently, either reference or variant? */
// Was the site called confidently, either reference or variant?
public boolean confidentlyCalled = false;
VariantCallContext(VariationCall variation, List<Genotype> genotypes, boolean confidentlyCalledP) {
this.variation = variation;
this.genotypes = genotypes;
VariantCallContext(VariantContext vc, boolean confidentlyCalledP) {
this.vc = vc;
this.confidentlyCalled = confidentlyCalledP;
}
/** blank variation and genotypes => we're a ref site */
// blank variant context => we're a ref site
VariantCallContext(boolean confidentlyCalledP) {
this.confidentlyCalled = confidentlyCalledP;
}

3
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/CallableBasesAnalysis.java
View File

@ -61,7 +61,8 @@ public class CallableBasesAnalysis extends BasicVariantAnalysis implements Genot
// we actually have a record here
if (!(eval instanceof VariantBackedByGenotype)) { // evaluation record isn't a genotype, die!
throw new RuntimeException("Evaluation track isn't backed by a Genotype!");
return null;
//throw new RuntimeException("Evaluation track isn't backed by a Genotype!");
}
all_calls++;

3
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/ChipConcordance.java
View File

@ -94,6 +94,9 @@ public abstract class ChipConcordance extends BasicVariantAnalysis {
}
public String inc(Map<String, Genotype> chips, Variation eval, char ref) {
// TODO -- needed to add this for now while we're moving over to VE2
if ( !(eval instanceof VariantBackedByGenotype) )
return null;
// each implementing class can decide whether the Variation is valid
assertVariationIsValid(eval);

5
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/GenotypeConcordance.java
View File

@ -50,8 +50,9 @@ public class GenotypeConcordance extends ChipConcordance implements GenotypeAnal
for ( Genotype eval : evals ) {
if ( providedChipName != null )
hash.put(providedChipName, eval);
else if ( eval instanceof SampleBacked )
hash.put(((SampleBacked)eval).getSampleName(), eval);
// TODO -- fix me in VE2
//else if ( eval instanceof SampleBacked )
// hash.put(((SampleBacked)eval).getSampleName(), eval);
else if ( rodNames.size() == 1 )
hash.put(rodNames.get(0), eval);
else

2
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/VariantCounter.java
View File

@ -33,7 +33,7 @@ public class VariantCounter extends BasicVariantAnalysis implements GenotypeAnal
// TODO -- break the het check out to a different module used only for single samples
if ( this.getMaster().evalContainsGenotypes && eval != null && eval.isBiallelic() && eval.isSNP() ) {
if ( this.getMaster().evalContainsGenotypes && eval != null && eval.isBiallelic() && eval.isSNP() && eval instanceof VariantBackedByGenotype ) {
List<Genotype> genotypes = ((VariantBackedByGenotype)eval).getGenotypes();
if ( genotypes.size() == 1 && genotypes.get(0).isHet() ) {
nHets++;

9
java/src/org/broadinstitute/sting/oneoffprojects/refdata/HapmapVCFROD.java
View File

@ -70,7 +70,8 @@ public class HapmapVCFROD extends BasicReferenceOrderedDatum implements Variatio
}
public List<Genotype> getGenotypes() {
return rod.getGenotypes();
return null;
//return rod.getGenotypes();
}
public String getReference() {
@ -130,7 +131,8 @@ public class HapmapVCFROD extends BasicReferenceOrderedDatum implements Variatio
}
public Genotype getCalledGenotype() {
return rod.getCalledGenotype();
return null;
//return rod.getCalledGenotype();
}
public char getReferenceForSNP() {
@ -138,7 +140,8 @@ public class HapmapVCFROD extends BasicReferenceOrderedDatum implements Variatio
}
public boolean hasGenotype(DiploidGenotype g) {
return rod.hasGenotype(g);
return false;
//return rod.hasGenotype(g);
}
public VCFHeader getHeader() {

6
java/src/org/broadinstitute/sting/oneoffprojects/walkers/BaseTransitionTableCalculatorJavaWalker.java
View File

@ -74,7 +74,7 @@ public class BaseTransitionTableCalculatorJavaWalker extends LocusWalker<Set<Bas
}
public Set<BaseTransitionTable> map( RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context ) {
ReadBackedPileup pileup = context.getPileup();
ReadBackedPileup pileup = context.getBasePileup();
Set<BaseTransitionTable> newCounts = null;
//System.out.println(pileup.getBases());
if ( baseIsUsable(tracker, ref, pileup, context) ) {
@ -360,9 +360,9 @@ public class BaseTransitionTableCalculatorJavaWalker extends LocusWalker<Set<Bas
if ( !BaseUtils.isRegularBase(ref.getBase()) )
return false;
VariantCallContext calls = ug.runGenotyper(tracker,ref,context);
if ( calls == null || calls.genotypes == null)
if ( calls == null || calls.vc == null)
return false;
return ( calls.genotypes.size() > 0 && !calls.genotypes.get(0).isVariant(ref.getBase()) && calls.genotypes.get(0).getNegLog10PError() > confidentRefThreshold );
return ( calls.vc.getNSamples() > 0 && calls.vc.getGenotype(0).isHomRef() && calls.vc.getGenotype(0).getNegLog10PError() > confidentRefThreshold );
}

11
java/src/org/broadinstitute/sting/oneoffprojects/walkers/annotator/ProportionOfNonrefBasesSupportingSNP.java
View File

@ -3,13 +3,12 @@ package org.broadinstitute.sting.oneoffprojects.walkers.annotator;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.walkers.annotator.VariantAnnotation;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.Pair;
import org.broadinstitute.sting.utils.BaseUtils;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.walkers.annotator.*;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import java.util.Map;
@ -30,14 +29,14 @@ public class ProportionOfNonrefBasesSupportingSNP implements VariantAnnotation {
1,VCFInfoHeaderLine.INFO_TYPE.Float,"Simple proportion of non-reference bases that are the SNP base");
}
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> context, Variation var) {
if ( ! var.isSNP() || ! var.isBiallelic() ) {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> context, VariantContext vc) {
if ( ! vc.isSNP() || ! vc.isBiallelic() ) {
return null;
} else {
Pair<Integer,Integer> totalNonref_totalSNP = new Pair<Integer,Integer>(0,0);
for ( String sample : context.keySet() ) {
ReadBackedPileup pileup = context.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).getPileup();
totalNonref_totalSNP = getNonrefAndSNP(pileup, ref.getBase(), var.getAlternativeBaseForSNP(), totalNonref_totalSNP);
ReadBackedPileup pileup = context.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).getBasePileup();
totalNonref_totalSNP = getNonrefAndSNP(pileup, ref.getBase(), vc.getAlternateAllele(0).toString().charAt(0), totalNonref_totalSNP);
}
if ( totalNonref_totalSNP.equals(new Pair<Integer,Integer>(0,0)) )

19
java/src/org/broadinstitute/sting/oneoffprojects/walkers/annotator/ProportionOfRefSecondBasesSupportingSNP.java
View File

@ -3,23 +3,16 @@ package org.broadinstitute.sting.oneoffprojects.walkers.annotator;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.walkers.annotator.VariantAnnotation;
import org.broadinstitute.sting.utils.BaseUtils;
import org.broadinstitute.sting.utils.StingException;
import org.broadinstitute.sting.utils.Pair;
import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import java.io.IOException;
import java.io.PrintWriter;
import java.io.FileOutputStream;
import java.util.List;
import java.util.Map;
import org.broadinstitute.sting.gatk.walkers.annotator.*;
/**
* Created by IntelliJ IDEA.
* User: chartl
@ -33,14 +26,14 @@ public class ProportionOfRefSecondBasesSupportingSNP implements VariantAnnotatio
public String getKeyName() { return KEY_NAME; }
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> context, Variation var) {
if ( ! var.isSNP() || ! var.isBiallelic() ) {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> context, VariantContext vc) {
if ( ! vc.isSNP() || ! vc.isBiallelic() ) {
return null;
} else {
Pair<Integer,Integer> totalAndSNPSupporting = new Pair<Integer,Integer>(0,0);
for ( String sample : context.keySet() ) {
ReadBackedPileup pileup = context.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).getPileup();
totalAndSNPSupporting = getTotalRefAndSNPSupportCounts(pileup, ref.getBase(), var.getAlternativeBaseForSNP(), totalAndSNPSupporting);
ReadBackedPileup pileup = context.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).getBasePileup();
totalAndSNPSupporting = getTotalRefAndSNPSupportCounts(pileup, ref.getBase(), vc.getAlternateAllele(0).toString().charAt(0), totalAndSNPSupporting);
}
if ( totalAndSNPSupporting.equals(new Pair<Integer,Integer>(0,0)) )

12
java/src/org/broadinstitute/sting/oneoffprojects/walkers/annotator/ProportionOfSNPSecondBasesSupportingRef.java
View File

@ -2,9 +2,9 @@ package org.broadinstitute.sting.oneoffprojects.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.walkers.annotator.VariantAnnotation;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import org.broadinstitute.sting.utils.Pair;
import org.broadinstitute.sting.utils.BaseUtils;
@ -12,7 +12,7 @@ import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import java.util.Map;
import org.broadinstitute.sting.gatk.walkers.annotator.*;
/**
* Created by IntelliJ IDEA.
* User: chartl
@ -29,14 +29,14 @@ public class ProportionOfSNPSecondBasesSupportingRef implements VariantAnnotatio
public boolean useZeroQualityReads() { return USE_MAPQ0_READS; }
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> context, Variation var) {
if ( ! var.isSNP() || ! var.isBiallelic() ) {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> context, VariantContext vc) {
if ( ! vc.isSNP() || ! vc.isBiallelic() ) {
return null;
} else {
Pair<Integer,Integer> totalAndSNPSupporting = new Pair<Integer,Integer>(0,0);
for ( String sample : context.keySet() ) {
ReadBackedPileup pileup = context.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).getPileup();
totalAndSNPSupporting = getTotalSNPandRefSupporting(pileup, ref.getBase(), var.getAlternativeBaseForSNP(), totalAndSNPSupporting);
ReadBackedPileup pileup = context.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).getBasePileup();
totalAndSNPSupporting = getTotalSNPandRefSupporting(pileup, ref.getBase(), vc.getAlternateAllele(0).toString().charAt(0), totalAndSNPSupporting);
}
if ( totalAndSNPSupporting.equals(new Pair<Integer,Integer>(0,0)) )

8
java/src/org/broadinstitute/sting/oneoffprojects/walkers/annotator/ThousandGenomesAnnotator.java
View File

@ -2,11 +2,11 @@ package org.broadinstitute.sting.oneoffprojects.walkers.annotator;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.refdata.utils.RODRecordList;
import org.broadinstitute.sting.gatk.walkers.annotator.VariantAnnotation;
import org.broadinstitute.sting.oneoffprojects.refdata.HapmapVCFROD;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
import java.util.Map;
@ -14,8 +14,8 @@ import java.util.Map;
/**
* IF THERE IS NO JAVADOC RIGHT HERE, YELL AT chartl
*
* @Author chartl
* @Date Feb 1, 2010
* @author chartl
* @date Feb 1, 2010
*/
public class ThousandGenomesAnnotator implements VariantAnnotation {
@ -28,7 +28,7 @@ public class ThousandGenomesAnnotator implements VariantAnnotation {
1,VCFInfoHeaderLine.INFO_TYPE.String,"Is this site seen in Pilot1 or Pilot2 of 1KG?");
}
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> context, Variation variation) {
public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> context, VariantContext vc) {
if ( tracker == null ) {
return null;
}

17
java/src/org/broadinstitute/sting/playground/gatk/walkers/DeNovoSNPWalker.java
View File

@ -3,14 +3,13 @@ package org.broadinstitute.sting.playground.gatk.walkers;
import net.sf.samtools.SAMRecord;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.Genotype;
import org.broadinstitute.sting.gatk.filters.ZeroMappingQualityReadFilter;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.refdata.VariationRod;
import org.broadinstitute.sting.gatk.walkers.*;
import org.broadinstitute.sting.gatk.walkers.genotyper.UnifiedGenotyper;
import org.broadinstitute.sting.gatk.walkers.genotyper.VariantCallContext;
import org.broadinstitute.sting.utils.genotype.*;
import org.broadinstitute.sting.utils.Pair;
import java.util.ArrayList;
import java.util.List;
@ -45,8 +44,8 @@ public class DeNovoSNPWalker extends RefWalker<String, Integer>{
}
public String map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
Variation child = (Variation)tracker.lookup("child", null);
Variation dbsnp = (Variation)tracker.lookup("dbSNP", null);
VariationRod child = (VariationRod)tracker.lookup("child", null);
VariationRod dbsnp = (VariationRod)tracker.lookup("dbSNP", null);
if (child != null) {
if (child.isSNP() && child.getNegLog10PError() > 5) { // BTR > 5
@ -79,13 +78,13 @@ public class DeNovoSNPWalker extends RefWalker<String, Integer>{
AlignmentContext parent2_subContext = new AlignmentContext(context.getLocation(), parent2_reads, parent2_offsets);
VariantCallContext parent2 = UG.map(tracker, ref, parent2_subContext);
if ( parent1 != null && parent1.genotypes != null && parent2 != null && parent2.genotypes != null ) {
Genotype parent1call = parent1.genotypes.get(0);
Genotype parent2call = parent2.genotypes.get(0);
if ( parent1 != null && parent1.vc != null && parent2 != null && parent2.vc != null ) {
Genotype parent1call = parent1.vc.getGenotype(0);
Genotype parent2call = parent2.vc.getGenotype(0);
if (!parent1call.isVariant(parent1call.getReference().charAt(0)) &&
if (parent1call.isHomRef() &&
parent1call.getNegLog10PError() > 5 &&
!parent2call.isVariant(parent2call.getReference().charAt(0)) &&
!parent2call.isHomRef() &&
parent2call.getNegLog10PError() > 5) {
double sumConfidences = 0.5 * (0.5 * child.getNegLog10PError() +

17
java/src/org/broadinstitute/sting/playground/gatk/walkers/LocusMismatchWalker.java
View File

@ -4,13 +4,14 @@ import org.broadinstitute.sting.gatk.walkers.*;
import org.broadinstitute.sting.gatk.walkers.genotyper.*;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.refdata.ReferenceOrderedDatum;
import org.broadinstitute.sting.gatk.refdata.VariationRod;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.Genotype;
import org.broadinstitute.sting.utils.cmdLine.Argument;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.genotype.*;
/**
@ -56,7 +57,7 @@ public class LocusMismatchWalker extends LocusWalker<String,Integer> implements
ReadBackedPileup pileup = context.getBasePileup();
if ( locusIsUsable(tracker, ref, pileup, context) ) {
Genotype g = getGenotype(tracker, ref, context);
if ( g != null && g.isPointGenotype() )
if ( g != null )
result = errorCounts( ref, pileup, g );
}
@ -105,15 +106,15 @@ public class LocusMismatchWalker extends LocusWalker<String,Integer> implements
}
return String.format("%s %c %10s %5.2f %d %d %d %s",
pileup.getLocation(), ref.getBase(),
getGenotypeClass(g, ref.getBase()), 10 * g.getNegLog10PError(),
getGenotypeClass(g), 10 * g.getNegLog10PError(),
usableDepth, nMismatches, qSumMismatches, baseCountString);
}
return null;
}
private String getGenotypeClass(Genotype g, char ref) {
if ( ! g.isVariant(ref) ) return "HOM-REF";
private String getGenotypeClass(Genotype g) {
if ( g.isHomRef() ) return "HOM-REF";
else if ( g.isHet() ) return "HET";
else if ( g.isHom() ) return "HOM-NONREF";
else throw new StingException("Unexpected genotype in getGenotypeClass " + g);
@ -142,7 +143,7 @@ public class LocusMismatchWalker extends LocusWalker<String,Integer> implements
private boolean notCoveredByVariations( RefMetaDataTracker tracker ) {
for ( ReferenceOrderedDatum datum : tracker.getAllRods() ) {
if ( datum instanceof Variation || datum instanceof Genotype ) {
if ( datum instanceof VariationRod || datum instanceof Genotype ) {
//System.out.printf("Ignoring site because of %s%n", datum);
return false;
}
@ -163,10 +164,10 @@ public class LocusMismatchWalker extends LocusWalker<String,Integer> implements
private Genotype getGenotype( RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context ) {
VariantCallContext calls = ug.runGenotyper(tracker,ref,context);
if ( calls == null || calls.variation == null || calls.genotypes == null )
if ( calls == null || calls.vc == null || calls.vc.getNSamples() == 0 || !calls.vc.isSNP() )
return null;
else {
return calls.genotypes.get(0);
return calls.vc.getGenotype(0);
}
}

66
java/src/org/broadinstitute/sting/playground/gatk/walkers/SnpCallRateByCoverageWalker.java
View File

@ -3,21 +3,17 @@ package org.broadinstitute.sting.playground.gatk.walkers;
import net.sf.samtools.SAMRecord;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.Genotype;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.refdata.RodGenotypeChipAsGFF;
import org.broadinstitute.sting.gatk.walkers.LocusWalker;
import org.broadinstitute.sting.gatk.walkers.genotyper.UnifiedGenotyper;
import org.broadinstitute.sting.gatk.walkers.genotyper.VariantCallContext;
import org.broadinstitute.sting.utils.genotype.DiploidGenotype;
import org.broadinstitute.sting.utils.BaseUtils;
import org.broadinstitute.sting.utils.ListUtils;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.Pair;
import org.broadinstitute.sting.utils.cmdLine.Argument;
import org.broadinstitute.sting.utils.genotype.*;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
/**
@ -45,7 +41,7 @@ public class SnpCallRateByCoverageWalker extends LocusWalker<List<String>, Strin
}
public boolean filter(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
return (BaseUtils.simpleBaseToBaseIndex(ref.getBase()) != -1 && context.getPileup().size() != 0);
return (BaseUtils.simpleBaseToBaseIndex(ref.getBase()) != -1 && context.getBasePileup().size() != 0);
}
public List<String> map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
@ -79,10 +75,11 @@ public class SnpCallRateByCoverageWalker extends LocusWalker<List<String>, Strin
AlignmentContext subContext = new AlignmentContext(context.getLocation(), sub_reads, sub_offsets);
VariantCallContext calls = UG.map(tracker, ref, subContext);
if (calls != null && calls.genotypes != null && calls.genotypes.size() > 0) {
Genotype call = calls.genotypes.get(0);
String callType = (call.isVariant(call.getReference().charAt(0))) ? ((call.isHom()) ? "HomozygousSNP" : "HeterozygousSNP") : "HomozygousReference";
GenotypeCalls.add(coverage+"\t"+coverage_available+"\t"+hc_genotype+"\t"+callType+"\t"+toGeliString(call));
if (calls != null && calls.vc != null && calls.vc.getNSamples() > 0) {
Genotype call = calls.vc.getGenotype(0);
String callType = (!call.isHomRef()) ? ((call.isHom()) ? "HomozygousSNP" : "HeterozygousSNP") : "HomozygousReference";
// TODO -- fixme: the old way of doing things isn't being supported anymore
GenotypeCalls.add(coverage+"\t"+coverage_available+"\t"+hc_genotype+"\t"+callType+"\t"+call.toString());
}
}
}
@ -105,55 +102,6 @@ public class SnpCallRateByCoverageWalker extends LocusWalker<List<String>, Strin
return "";
}
// a method to support getting the geli string, since the AlleleFrequencyEstimate is going away
public String toGeliString (Genotype locus) {
double posteriors[];
int readDepth = -1;
double nextVrsBest = 0;
double nextVrsRef = 0;
char ref = locus.getReference().charAt(0);
if (locus instanceof ReadBacked) {
readDepth = ((ReadBacked)locus).getReadCount();
}
if (!(locus instanceof PosteriorsBacked)) {
posteriors = new double[10];
Arrays.fill(posteriors, 0.0);
} else {
posteriors = ((PosteriorsBacked) locus).getPosteriors();
double[] lks;
lks = Arrays.copyOf(posteriors,posteriors.length);
Arrays.sort(lks);
nextVrsBest = lks[9] - lks[8];
if (ref != 'X') {
int index = (DiploidGenotype.valueOf(Utils.dupString(ref,2)).ordinal());
nextVrsRef = lks[9] - posteriors[index];
}
}
// we have to calcuate our own
return new String(String.format("%s %16d %c %8d %d %s %.6f %.6f %6.6f %6.6f %6.6f %6.6f %6.6f %6.6f %6.6f %6.6f %6.6f %6.6f",
locus.getLocation().getContig(),
locus.getLocation().getStart(),
ref,
readDepth,
-1,
locus.getBases(),
nextVrsRef,
nextVrsBest,
posteriors[0],
posteriors[1],
posteriors[2],
posteriors[3],
posteriors[4],
posteriors[5],
posteriors[6],
posteriors[7],
posteriors[8],
posteriors[9]));
}
}

4
java/src/org/broadinstitute/sting/playground/gatk/walkers/contamination/FindContaminatingReadGroupsWalker.java
View File

@ -79,8 +79,8 @@ public class FindContaminatingReadGroupsWalker extends LocusWalker<Integer, Inte
if (altBalance > 0.70) {
VariantCallContext ugResult = ug.runGenotyper(tracker, ref, context);
if (ugResult != null && ugResult.genotypes != null && ugResult.genotypes.size() > 0) {
return ugResult.genotypes.get(0).isHet();
if (ugResult != null && ugResult.vc != null && ugResult.vc.getNSamples() > 0) {
return ugResult.vc.getGenotype(0).isHet();
}
}

10
java/src/org/broadinstitute/sting/playground/gatk/walkers/secondaryBases/SecondaryBaseTransitionTableWalker.java
View File

@ -2,6 +2,7 @@ package org.broadinstitute.sting.playground.gatk.walkers.secondaryBases;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.Genotype;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.walkers.LocusWalker;
import org.broadinstitute.sting.gatk.walkers.Reference;
@ -9,7 +10,6 @@ import org.broadinstitute.sting.gatk.walkers.Window;
import org.broadinstitute.sting.gatk.walkers.genotyper.*;
import org.broadinstitute.sting.playground.utils.NamedTable;
import org.broadinstitute.sting.utils.BaseUtils;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
@ -50,12 +50,12 @@ public class SecondaryBaseTransitionTableWalker extends LocusWalker<Integer, Int
if (contextBases.length == 3 && refBase != 'N' && pileup.getBases() != null && pileup.getSecondaryBases() != null) {
VariantCallContext ugResult = ug.runGenotyper(tracker,ref,context);
if (ugResult != null && ugResult.variation != null) {
Genotype res = ugResult.genotypes.get(0);
String call = res.getBases();
if (ugResult != null && ugResult.vc != null) {
Genotype res = ugResult.vc.getGenotype(0);
String call = res.getGenotypeString();
String type;
String alleleBalance = "N/A";
if (!res.isVariant(refBase)) {
if (res.isHomRef()) {
type = "homref";
}
else if (!res.isHet()) {type = "homvar";}

5
java/src/org/broadinstitute/sting/playground/gatk/walkers/vcftools/VCFSubsetWalker.java
View File

@ -80,8 +80,7 @@ public class VCFSubsetWalker extends RodWalker<ArrayList<VCFRecord>, VCFWriter>
private VCFRecord subsetRecord(VCFRecord record) {
ArrayList<VCFGenotypeRecord> genotypeRecords = new ArrayList<VCFGenotypeRecord>();
for (int i = 0; i < record.getGenotypes().size(); i++) {
VCFGenotypeRecord gr = (VCFGenotypeRecord)record.getGenotypes().get(i);
for ( VCFGenotypeRecord gr : record.getVCFGenotypeRecords() ) {
//if (gr.getSampleName().equalsIgnoreCase(SAMPLE)) {
if (SAMPLES.contains(gr.getSampleName())) {
@ -109,7 +108,7 @@ public class VCFSubsetWalker extends RodWalker<ArrayList<VCFRecord>, VCFWriter>
VCFRecord subset = subsetRecord(record);
boolean isVariant = false;
for (VCFGenotypeEncoding ge : ((VCFGenotypeRecord)subset.getGenotypes().get(0)).getAlleles()) {
for ( VCFGenotypeEncoding ge : subset.getVCFGenotypeRecords().get(0).getAlleles() ) {
if (!record.getReference().equals(ge.getBases())) {
isVariant = true;
}

17
java/src/org/broadinstitute/sting/utils/SampleUtils.java
View File

@ -5,8 +5,6 @@ import net.sf.samtools.SAMFileHeader;
import java.util.*;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.genotype.SampleBacked;
import org.broadinstitute.sting.utils.genotype.vcf.VCFReader;
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.gatk.refdata.ReferenceOrderedData;
@ -42,21 +40,6 @@ public class SampleUtils {
return samples;
}
/**
* get the samples names from genotype objects if they are backed by samples
*
* @param genotypes the genotype list
* @return list of strings representing the sample names
*/
public static List<String> getGenotypeSamples(List<Genotype> genotypes) {
List<String> samples = new ArrayList<String>();
for ( Genotype genotype : genotypes ) {
if ( genotype instanceof SampleBacked )
samples.add(((SampleBacked)genotype).getSampleName());
}
return samples;
}
/**
* Gets all of the unique sample names from all VCF rods input by the user
*

55
java/src/org/broadinstitute/sting/utils/genotype/AlleleConstrainedGenotype.java
View File

@ -1,55 +0,0 @@
package org.broadinstitute.sting.utils.genotype;
/**
* @author ebanks
* <p/>
* Class AlleleConstrainedGenotype
* <p/>
* A genotype that can have only one of 3 genotypes AA,AB,BB
*/
public abstract class AlleleConstrainedGenotype implements Genotype {
protected final static char NO_CONSTRAINT = 'N';
private char ref = NO_CONSTRAINT;
private char alt = NO_CONSTRAINT;
public AlleleConstrainedGenotype(String ref) {
this.ref = ref.charAt(0);
}
/**
* set the allowed alternate allele
*
* @param alt the alternate allele
*/
public void setAlternateAllele(char alt) {
this.alt = alt;
}
/**
* @return returns the allowed alternate allele
*/
public char getAlternateAllele() {
return alt;
}
/**
*
* @return returns the best genotype
*/
protected abstract DiploidGenotype getBestGenotype();
/**
* get the bases that represent this
*
* @return the bases, as a string
*/
public String getBases() {
DiploidGenotype g = getBestGenotype();
if ( alt != NO_CONSTRAINT && ((g.base1 != ref && g.base1 != alt) || (g.base2 != ref && g.base2 != alt)) )
throw new IllegalStateException("The best genotype " + g + " is composed of an allele that is not " + ref + " or " + alt);
return g.toString();
}
}

25
java/src/org/broadinstitute/sting/utils/genotype/ArbitraryFieldsBacked.java
View File

@ -1,25 +0,0 @@
package org.broadinstitute.sting.utils.genotype;
import java.util.Map;
/**
* @author ebanks
* Interface AlleleBalanceBacked
*
* this interface indicates that the genotype is
* backed up by allele balance information.
*/
public interface ArbitraryFieldsBacked {
/**
*
* @return returns the list of arbitrary fields to emit
*/
public Map<String, String> getFields();
/**
*
* @param fields a list of arbitrary fields to emit
*/
public void setFields(Map<String, String> fields);
}

2
java/src/org/broadinstitute/sting/utils/genotype/BasicGenotypeBackedVariation.java
View File

@ -11,7 +11,7 @@ import java.util.*;
* <p/>
* represents a genotype-backed Variation.
*/
public class BasicGenotypeBackedVariation implements Variation, VariantBackedByGenotype, ConfidenceBacked {
public class BasicGenotypeBackedVariation implements Variation, VariantBackedByGenotype {
// the discovery lod score
private double mConfidence = 0.0;

45
java/src/org/broadinstitute/sting/utils/genotype/CalledGenotype.java 100755
View File

@ -0,0 +1,45 @@
package org.broadinstitute.sting.utils.genotype;
import org.broadinstitute.sting.gatk.contexts.variantcontext.*;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import java.util.*;
/**
* This class emcompasses all the basic information about a genotype. It is immutable.
*
* @author ebanks
*/
public class CalledGenotype extends MutableGenotype {
public static final String LIKELIHOODS_ATTRIBUTE_KEY = "Likelihoods";
public static final String POSTERIORS_ATTRIBUTE_KEY = "Posteriors";
public static final String READBACKEDPILEUP_ATTRIBUTE_KEY = "ReadBackedPileup";
public CalledGenotype(String sampleName, List<Allele> alleles, double negLog10PError, Set<String> filters, Map<String, ?> attributes, boolean genotypesArePhased) {
super(sampleName, alleles, negLog10PError, filters, attributes, genotypesArePhased);
}
public CalledGenotype(String sampleName, List<Allele> alleles, double negLog10PError) {
super(sampleName, alleles, negLog10PError);
}
public CalledGenotype(String sampleName, List<Allele> alleles) {
super(sampleName, alleles, NO_NEG_LOG_10PERROR, null, null, false);
}
// ---------------------------------------------------------------------------------------------------------
//
// routines to modify useful attribute fields
//
// ---------------------------------------------------------------------------------------------------------
public void setLikelihoods(double likelihoods[]) { putAttribute(LIKELIHOODS_ATTRIBUTE_KEY, likelihoods); }
public void setPosteriors(double posteriors[]) { putAttribute(POSTERIORS_ATTRIBUTE_KEY, posteriors); }
public void setReadBackedPileup(ReadBackedPileup pileup) { putAttribute(READBACKEDPILEUP_ATTRIBUTE_KEY, pileup); }
public double[] getLikelihoods() { return (double[])getAttribute(LIKELIHOODS_ATTRIBUTE_KEY); }
public double[] getPosteriors() { return (double[])getAttribute(POSTERIORS_ATTRIBUTE_KEY); }
public ReadBackedPileup getReadBackedPileup() { return (ReadBackedPileup)getAttribute(READBACKEDPILEUP_ATTRIBUTE_KEY); }
}

24
java/src/org/broadinstitute/sting/utils/genotype/ConfidenceBacked.java
View File

@ -1,24 +0,0 @@
package org.broadinstitute.sting.utils.genotype;
/**
* @author ebanks
* Interface ConfidenceBacked
*
* this interface indicates that the genotype is
* backed up by sample information.
*/
public interface ConfidenceBacked {
/**
*
* @return returns the confidence for this genotype
*/
public double getConfidence();
/**
*
* @param confidence the confidence for this genotype
*/
public void setConfidence(double confidence);
}

30
java/src/org/broadinstitute/sting/utils/genotype/GenotypeWriter.java
View File

@ -1,6 +1,7 @@
package org.broadinstitute.sting.utils.genotype;
import java.util.List;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
/*
* Copyright (c) 2009 The Broad Institute
@ -28,7 +29,7 @@ import java.util.List;
*/
/**
* @author aaron
* @author aaron, ebanks
* <p/>
* Class GenotypeWriter
* <p/>
@ -36,31 +37,12 @@ import java.util.List;
*/
public interface GenotypeWriter {
/**
* Add a genotype, given a genotype locus
* @param call the locus to add
* Add a genotype, given a variant context
* @param vc the variant context representing the call to add
*/
public void addGenotypeCall(Genotype call);
/**
* add a no call to the genotype file, if supported.
*
* @param position the position to add the no call at
*/
public void addNoCall(int position);
public void addCall(VariantContext vc);
/** finish writing, closing any open files. */
public void close();
/**
* add a multi-sample call if we support it
* @param genotypes the list of genotypes, that are backed by sample information
* @param variation the variation
*/
public void addMultiSampleCall(List<Genotype> genotypes, VariationCall variation);
/**
* @return true if we support multisample, false otherwise
*/
public boolean supportsMultiSample();
}

44
java/src/org/broadinstitute/sting/utils/genotype/GenotypeWriterFactory.java
View File

@ -1,7 +1,6 @@
package org.broadinstitute.sting.utils.genotype;
import net.sf.samtools.SAMFileHeader;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.StingException;
import org.broadinstitute.sting.utils.genotype.geli.*;
import org.broadinstitute.sting.utils.genotype.glf.*;
@ -77,47 +76,4 @@ public class GenotypeWriterFactory {
}
// nothing to do for GELI TEXT
}
/**
* create a genotype call
* @param format the format
* @param ref the reference base
* @param loc the location
* @return an unpopulated genotype call object
*/
public static GenotypeCall createSupportedGenotypeCall(GENOTYPE_FORMAT format, String ref, GenomeLoc loc) {
switch (format) {
case VCF:
return new VCFGenotypeCall(ref, loc);
case GELI:
case GELI_BINARY:
return new GeliGenotypeCall(ref, loc);
case GLF:
return new GLFGenotypeCall(ref, loc);
default:
throw new StingException("Genotype format " + format.toString() + " is not implemented");
}
}
/**
* create a genotype locus data object
* @param format the format
* @param ref the reference base
* @param loc the location
* @param type the variant type
* @return an unpopulated genotype locus data object
*/
public static VariationCall createSupportedCall(GENOTYPE_FORMAT format, char ref, GenomeLoc loc, Variation.VARIANT_TYPE type) {
switch (format) {
case VCF:
return new VCFVariationCall(ref, loc, type);
case GELI:
case GELI_BINARY:
return null;
case GLF:
return null;
default:
throw new StingException("Genotype format " + format.toString() + " is not implemented");
}
}
}

24
java/src/org/broadinstitute/sting/utils/genotype/IDBacked.java
View File

@ -1,24 +0,0 @@
package org.broadinstitute.sting.utils.genotype;
/**
* @author ebanks
* Interface IDBacked
*
* this interface indicates that the genotype is
* backed up by sample information.
*/
public interface IDBacked {
/**
*
* @return returns the ID for this genotype
*/
public String getID();
/**
*
* @param id the ID for this genotype
*/
public void setID(String id);
}

24
java/src/org/broadinstitute/sting/utils/genotype/LikelihoodsBacked.java
View File

@ -1,24 +0,0 @@
package org.broadinstitute.sting.utils.genotype;
/**
* @author aaron
* Interface LikelihoodsBacked
*
* this interface indicates that the genotype is
* backed up by genotype likelihood information.
*/
public interface LikelihoodsBacked {
/**
*
* @return the likelihood information for this genotype
*/
public double[] getLikelihoods();
/**
*
* @param likelihoods the likelihoods for this genotype
*/
public void setLikelihoods(double[] likelihoods);
}

25
java/src/org/broadinstitute/sting/utils/genotype/PosteriorsBacked.java
View File

@ -1,25 +0,0 @@
package org.broadinstitute.sting.utils.genotype;
/**
* @author aaron
* Interface PosteriorsBacked
*
* This interface indicates that the genotype is backed up by *diploid* genotype posterior information.
* The posteriors array is expected to have 10 entries (one for each of the possible diploid genotypes).
*/
public interface PosteriorsBacked {
/**
*
* @return the likelihood information for this genotype
*/
public double[] getPosteriors();
/**
*
* @param posteriors the posteriors for this genotype
*/
public void setPosteriors(double[] posteriors);
}

24
java/src/org/broadinstitute/sting/utils/genotype/SLODBacked.java
View File

@ -1,24 +0,0 @@
package org.broadinstitute.sting.utils.genotype;
/**
* @author ebanks
* Interface SLODBacked
*
* this interface indicates that the genotype is
* backed up by strand lod information.
*/
public interface SLODBacked {
/**
*
* @return returns the strand lod for this genotype or null if not set
*/
public Double getSLOD();
/**
*
* @param slod the strand lod for this genotype
*/
public void setSLOD(double slod);
}

24
java/src/org/broadinstitute/sting/utils/genotype/SampleBacked.java
View File

@ -1,24 +0,0 @@
package org.broadinstitute.sting.utils.genotype;
/**
* @author aaron
* Interface SampleBacked
*
* this interface indicates that the genotype is
* backed up by sample information.
*/
public interface SampleBacked {
/**
*
* @return returns the sample name for this genotype
*/
public String getSampleName();
/**
*
* @param name the sample name for this genotype
*/
public void setSampleName(String name);
}

115
java/src/org/broadinstitute/sting/utils/genotype/geli/GeliAdapter.java
View File

@ -3,13 +3,15 @@ package org.broadinstitute.sting.utils.genotype.geli;
import edu.mit.broad.picard.genotype.geli.GeliFileWriter;
import edu.mit.broad.picard.genotype.geli.GenotypeLikelihoods;
import net.sf.samtools.SAMFileHeader;
import net.sf.samtools.SAMRecord;
import net.sf.samtools.SAMSequenceRecord;
import org.broadinstitute.sting.utils.GenomeLocParser;
import org.broadinstitute.sting.utils.genotype.*;
import org.broadinstitute.sting.utils.genotype.LikelihoodObject;
import org.broadinstitute.sting.utils.genotype.CalledGenotype;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.gatk.contexts.variantcontext.*;
import java.io.File;
import java.util.List;
/*
@ -38,7 +40,7 @@ import java.util.List;
*/
/**
* @author aaron
* @author aaron, ebanks
* @version 1.0
* <p/>
* Class GeliAdapter
@ -81,12 +83,12 @@ public class GeliAdapter implements GeliGenotypeWriter {
* @param readCount the read count
* @param likelihoods the likelihoods of each of the possible alleles
*/
private void addGenotypeCall(SAMSequenceRecord contig,
int position,
char referenceBase,
double maxMappingQuality,
int readCount,
LikelihoodObject likelihoods) {
private void addCall(SAMSequenceRecord contig,
int position,
char referenceBase,
double maxMappingQuality,
int readCount,
LikelihoodObject likelihoods) {
GenotypeLikelihoods lk = likelihoods.convertToGenotypeLikelihoods(writer.getFileHeader(), contig.getSequenceIndex(), position, (byte) referenceBase);
lk.setNumReads(readCount);
@ -94,22 +96,6 @@ public class GeliAdapter implements GeliGenotypeWriter {
writer.addGenotypeLikelihoods(lk);
}
/**
* add a variable length call to the genotyper
*
* @param contig the contig you're calling in
* @param position the position on the genome
* @param rmsMapQuals the root mean square of the mapping qualities
* @param readDepth the read depth
* @param refBase the reference base
* @param firstHomZyg the first homozygous indel
* @param secondHomZyg the second homozygous indel (if present, null if not)
* @param hetLikelihood the heterozygous likelihood
*/
public void addVariableLengthCall(SAMSequenceRecord contig, int position, float rmsMapQuals, int readDepth, char refBase, IndelLikelihood firstHomZyg, IndelLikelihood secondHomZyg, byte hetLikelihood) {
throw new UnsupportedOperationException("Geli format does not support variable length allele calls");
}
public void addGenotypeLikelihoods(GenotypeLikelihoods gl) {
if ( writer == null )
throw new IllegalStateException("The Geli Header must be written before records can be added");
@ -118,45 +104,52 @@ public class GeliAdapter implements GeliGenotypeWriter {
}
/**
* Add a genotype, given a genotype call
* Add a genotype, given a variant context
*
* @param call the call to add
* @param vc the variant context representing the call to add
*/
public void addGenotypeCall(Genotype call) {
public void addCall(VariantContext vc) {
if ( writer == null )
throw new IllegalStateException("The Geli Header must be written before calls can be added");
if ( !(call instanceof GeliGenotypeCall) )
throw new IllegalArgumentException("Only GeliGenotypeCalls should be passed in to the Geli writers");
GeliGenotypeCall gCall = (GeliGenotypeCall)call;
char ref = vc.getReference().toString().charAt(0);
if ( vc.getNSamples() != 1 )
throw new IllegalArgumentException("The Geli format does not support multi-sample or no-calls");
char ref = gCall.getReference().charAt(0);
int readCount = gCall.getReadCount();
Genotype genotype = vc.getGenotypes().values().iterator().next();
if ( genotype.isNoCall() )
throw new IllegalArgumentException("The Geli format does not support no-calls");
ReadBackedPileup pileup;
double[] posteriors;
if ( genotype instanceof CalledGenotype ) {
pileup = ((CalledGenotype)genotype).getReadBackedPileup();
posteriors = ((CalledGenotype)genotype).getPosteriors();
} else {
pileup = (ReadBackedPileup)genotype.getAttribute(CalledGenotype.READBACKEDPILEUP_ATTRIBUTE_KEY);
posteriors = (double[])genotype.getAttribute(CalledGenotype.POSTERIORS_ATTRIBUTE_KEY);
}
if ( posteriors == null )
throw new IllegalArgumentException("The Geli format requires posteriors");
int readCount = 0;
double maxMappingQual = 0;
if ( gCall.getPileup() != null ) {
List<SAMRecord> recs = gCall.getPileup().getReads();
for (SAMRecord rec : recs) {
if (maxMappingQual < rec.getMappingQuality()) maxMappingQual = rec.getMappingQuality();
if ( pileup != null ) {
readCount = pileup.size();
for (PileupElement p : pileup ) {
if ( maxMappingQual < p.getMappingQual() )
maxMappingQual = p.getMappingQual();
}
}
double[] posteriors = gCall.getPosteriors();
LikelihoodObject obj = new LikelihoodObject(posteriors, LikelihoodObject.LIKELIHOOD_TYPE.LOG);
this.addGenotypeCall(GenomeLocParser.getContigInfo(gCall.getLocation().getContig()),
(int)gCall.getLocation().getStart(),
ref,
maxMappingQual,
readCount,
obj);
}
/**
* add a no call to the genotype file, if supported.
*
* @param position the position
*/
public void addNoCall(int position) {
throw new UnsupportedOperationException("Geli format does not support no-calls");
addCall(GenomeLocParser.getContigInfo(vc.getLocation().getContig()),
(int)vc.getLocation().getStart(),
ref,
maxMappingQual,
readCount,
obj);
}
/** finish writing, closing any open files. */
@ -165,18 +158,4 @@ public class GeliAdapter implements GeliGenotypeWriter {
this.writer.close();
}
}
/**
* add a multi-sample call if we support it
*
* @param genotypes the list of genotypes
*/
public void addMultiSampleCall( List<Genotype> genotypes, VariationCall metadata) {
throw new UnsupportedOperationException("Geli binary doesn't support multisample calls");
}
/** @return true if we support multisample, false otherwise */
public boolean supportsMultiSample() {
return false;
}
}

300
java/src/org/broadinstitute/sting/utils/genotype/geli/GeliGenotypeCall.java
View File

@ -1,300 +0,0 @@
package org.broadinstitute.sting.utils.genotype.geli;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.genotype.*;
import java.util.Arrays;
/**
* @author ebanks
* <p/>
* Class GeliGenotypeCall
* <p/>
* The implementation of the genotype interface, specific to Geli
*/
public class GeliGenotypeCall extends AlleleConstrainedGenotype implements GenotypeCall, ReadBacked, PosteriorsBacked {
private final char mRefBase;
private final GenomeLoc mLocation;
private ReadBackedPileup mPileup = null;
private double[] mPosteriors;
private Variation mVariation = null;
// the reference genotype, the best genotype, and the next best genotype, lazy loaded
private DiploidGenotype mRefGenotype = null;
private DiploidGenotype mBestGenotype = null;
private DiploidGenotype mNextGenotype = null;
/**
* Generate a single sample genotype object
*
* @param ref the ref character
* @param loc the genome loc
*/
public GeliGenotypeCall(String ref, GenomeLoc loc) {
super(ref);
mRefBase = ref.charAt(0);
mLocation = loc;
// fill in empty data
mPosteriors = new double[10];
Arrays.fill(mPosteriors, Double.MIN_VALUE);
}
public GeliGenotypeCall(String ref, GenomeLoc loc, String genotype, double negLog10PError) {
super(ref);
mRefBase = ref.charAt(0);
mLocation = loc;
mBestGenotype = DiploidGenotype.valueOf(genotype);
mRefGenotype = DiploidGenotype.createHomGenotype(mRefBase);
mNextGenotype = mRefGenotype;
// set general posteriors to min double value
mPosteriors = new double[10];
Arrays.fill(mPosteriors, Double.MIN_VALUE);
// set the ref to PError
mPosteriors[mRefGenotype.ordinal()] = -1.0 * negLog10PError;
// set the best genotype to zero (need to do this after ref in case ref=best)
mPosteriors[mBestGenotype.ordinal()] = 0.0;
// choose a smart next best genotype and set it to PError
if ( mBestGenotype == mRefGenotype )
mNextGenotype = DiploidGenotype.valueOf(BaseUtils.simpleComplement(genotype));
else
mNextGenotype = mRefGenotype;
mPosteriors[mNextGenotype.ordinal()] = -1.0 * negLog10PError;
}
public void setPileup(ReadBackedPileup pileup) {
mPileup = pileup;
}
public void setPosteriors(double[] posteriors) {
mPosteriors = posteriors;
}
public void setVariation(Variation variation) {
this.mVariation = variation;
}
public void setGenotype(DiploidGenotype genotype) {
; // do nothing: geli uses diploid posteriors to calculate the genotype
}
public void setNegLog10PError(double value) {
; // do nothing: geli uses diploid posteriors to calculate the P(error)
}
@Override
public boolean equals(Object other) {
lazyEval();
if (other == null)
return false;
if (other instanceof GeliGenotypeCall) {
GeliGenotypeCall otherCall = (GeliGenotypeCall) other;
if (!this.mBestGenotype.equals(otherCall.mBestGenotype))
return false;
return (this.mRefBase == otherCall.mRefBase);
}
return false;
}
public String toString() {
lazyEval();
return String.format("%s best=%s cmp=%s ref=%s depth=%d negLog10PError=%.2f",
getLocation(), mBestGenotype, mRefGenotype, mRefBase, getReadCount(), getNegLog10PError());
}
private void lazyEval() {
if (mBestGenotype == null) {
char ref = this.getReference().charAt(0);
char alt = this.getAlternateAllele();
mRefGenotype = DiploidGenotype.createHomGenotype(ref);
// if we are constrained to a single alternate allele, use only that one
if ( alt != AlleleConstrainedGenotype.NO_CONSTRAINT ) {
DiploidGenotype hetGenotype = ref < alt ? DiploidGenotype.valueOf(String.valueOf(ref) + String.valueOf(alt)) : DiploidGenotype.valueOf(String.valueOf(alt) + String.valueOf(ref));
DiploidGenotype homGenotype = DiploidGenotype.createHomGenotype(alt);
boolean hetOverHom = mPosteriors[hetGenotype.ordinal()] > mPosteriors[homGenotype.ordinal()];
boolean hetOverRef = mPosteriors[hetGenotype.ordinal()] > mPosteriors[mRefGenotype.ordinal()];
boolean homOverRef = mPosteriors[homGenotype.ordinal()] > mPosteriors[mRefGenotype.ordinal()];
if ( hetOverHom ) {
mBestGenotype = (hetOverRef ? hetGenotype : mRefGenotype);
mNextGenotype = (!hetOverRef ? hetGenotype : (homOverRef ? homGenotype : mRefGenotype));
} else {
mBestGenotype = (homOverRef ? homGenotype : mRefGenotype);
mNextGenotype = (!homOverRef ? homGenotype : (hetOverRef ? hetGenotype : mRefGenotype));
}
} else {
Integer sorted[] = Utils.SortPermutation(mPosteriors);
mBestGenotype = DiploidGenotype.values()[sorted[DiploidGenotype.values().length - 1]];
mNextGenotype = DiploidGenotype.values()[sorted[DiploidGenotype.values().length - 2]];
}
}
}
/**
* get the confidence we have
*
* @return get the one minus error value
*/
public double getNegLog10PError() {
return Math.abs(mPosteriors[getBestGenotype().ordinal()] - mPosteriors[getNextBest().ordinal()]);
}
// get the best genotype
protected DiploidGenotype getBestGenotype() {
lazyEval();
return mBestGenotype;
}
// get the alternate genotype
private DiploidGenotype getNextBest() {
lazyEval();
return mNextGenotype;
}
// get the ref genotype
private DiploidGenotype getRefGenotype() {
lazyEval();
return mRefGenotype;
}
/**
* get the bases that represent this
*
* @return the bases, as a string
*/
public String getBases() {
return getBestGenotype().toString();
}
/**
* get the ploidy
*
* @return the ploidy value
*/
public int getPloidy() {
return 2;
}
/**
* Returns true if both observed alleles are the same (regardless of whether they are ref or alt)
*
* @return true if we're homozygous, false otherwise
*/
public boolean isHom() {
return getBestGenotype().isHom();
}
/**
* Returns true if observed alleles differ (regardless of whether they are ref or alt)
*
* @return true if we're het, false otherwise
*/
public boolean isHet() {
return !isHom();
}
public boolean isNoCall() { return false; }
/**
* Location of this genotype on the reference (on the forward strand). If the allele is insertion/deletion, the first inserted/deleted base
* is located right <i>after</i> the specified location
*
* @return position on the genome wrapped in GenomeLoc object
*/
public GenomeLoc getLocation() {
return this.mLocation;
}
/**
* returns true if the genotype is a point genotype, false if it's a indel / deletion
*
* @return true is a SNP
*/
public boolean isPointGenotype() {
return true;
}
/**
* given the reference, are we a variant? (non-ref)
*
* @param ref the reference base or bases
*
* @return true if we're a variant
*/
public boolean isVariant(char ref) {
return !Utils.dupString(ref, 2).equals(getBestGenotype().toString());
}
/**
* are we a variant? (non-ref)
*
* @return true if we're a variant
*/
public boolean isVariant() {
return isVariant(mRefBase);
}
/**
*
* @return return this genotype as a variant
*/
public Variation toVariation(char ref) {
if ( mVariation == null ) {
BasicGenotypeBackedVariation var = new BasicGenotypeBackedVariation(ref, mLocation, isVariant() ? Variation.VARIANT_TYPE.SNP : Variation.VARIANT_TYPE.REFERENCE);
double confidence = Math.abs(mPosteriors[getBestGenotype().ordinal()] - mPosteriors[getRefGenotype().ordinal()]);
var.setConfidence(confidence);
if ( isVariant() )
var.addAlternateAllele(Character.toString(mBestGenotype.base1 != ref ? mBestGenotype.base1 : mBestGenotype.base2));
mVariation = var;
}
return mVariation;
}
/**
* get the pileup that backs this genotype call
*
* @return a pileup
*/
public ReadBackedPileup getPileup() {
return mPileup;
}
/**
* get the count of reads
*
* @return the number of reads we're backed by
*/
public int getReadCount() {
return (mPileup != null ? mPileup.getReads().size() : 0);
}
/**
* get the reference character
*
* @return the reference character
*/
public String getReference() {
return Character.toString(mRefBase);
}
/**
* get the posteriors
*
* @return the posteriors
*/
public double[] getPosteriors() {
return mPosteriors;
}
}

89
java/src/org/broadinstitute/sting/utils/genotype/geli/GeliTextWriter.java
View File

@ -1,16 +1,19 @@
package org.broadinstitute.sting.utils.genotype.geli;
import net.sf.samtools.SAMRecord;
import net.sf.samtools.SAMFileHeader;
import org.broadinstitute.sting.utils.StingException;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.genotype.*;
import org.broadinstitute.sting.gatk.contexts.variantcontext.*;
import java.io.File;
import java.io.FileNotFoundException;
import java.io.PrintStream;
import java.io.PrintWriter;
import java.util.Arrays;
import java.util.List;
import java.util.ArrayList;
import java.util.Collections;
import edu.mit.broad.picard.genotype.geli.GenotypeLikelihoods;
@ -56,18 +59,34 @@ public class GeliTextWriter implements GeliGenotypeWriter {
}
/**
* Add a genotype, given a call
* Add a genotype, given a variant context
*
* @param call the call to add
* @param vc the variant context representing the call to add
*/
public void addGenotypeCall(Genotype call) {
if ( !(call instanceof GeliGenotypeCall) )
throw new IllegalArgumentException("Only GeliGenotypeCalls should be passed in to the Geli writers");
GeliGenotypeCall gCall = (GeliGenotypeCall)call;
public void addCall(VariantContext vc) {
char ref = gCall.getReference().charAt(0);
char ref = vc.getReference().toString().charAt(0);
if ( vc.getNSamples() != 1 )
throw new IllegalArgumentException("The Geli format does not support multi-sample or no-calls");
org.broadinstitute.sting.gatk.contexts.variantcontext.Genotype genotype = vc.getGenotypes().values().iterator().next();
if ( genotype.isNoCall() )
throw new IllegalArgumentException("The Geli format does not support no-calls");
ReadBackedPileup pileup;
double[] posteriors;
if ( genotype instanceof CalledGenotype ) {
pileup = ((CalledGenotype)genotype).getReadBackedPileup();
posteriors = ((CalledGenotype)genotype).getPosteriors();
} else {
pileup = (ReadBackedPileup)genotype.getAttribute(CalledGenotype.READBACKEDPILEUP_ATTRIBUTE_KEY);
posteriors = (double[])genotype.getAttribute(CalledGenotype.POSTERIORS_ATTRIBUTE_KEY);
}
if ( posteriors == null )
throw new IllegalArgumentException("The Geli format requires posteriors");
double[] posteriors = gCall.getPosteriors();
double[] lks;
lks = Arrays.copyOf(posteriors, posteriors.length);
Arrays.sort(lks);
@ -77,20 +96,31 @@ public class GeliTextWriter implements GeliGenotypeWriter {
if (ref != 'X')
nextVrsRef = lks[9] - posteriors[DiploidGenotype.createHomGenotype(ref).ordinal()];
int readCount = 0;
double maxMappingQual = 0;
List<SAMRecord> recs = gCall.getPileup().getReads();
int readDepth = recs.size();
for (SAMRecord rec : recs) {
if (maxMappingQual < rec.getMappingQuality()) maxMappingQual = rec.getMappingQuality();
if ( pileup != null ) {
readCount = pileup.size();
for (PileupElement p : pileup ) {
if ( maxMappingQual < p.getMappingQual() )
maxMappingQual = p.getMappingQual();
}
}
ArrayList<Character> alleles = new ArrayList<Character>();
for ( Allele a : genotype.getAlleles() )
alleles.add(a.toString().charAt(0));
Collections.sort(alleles);
StringBuffer sb = new StringBuffer();
for ( Character base : alleles )
sb.append(base);
mWriter.println(String.format("%s %16d %c %8d %.0f %s %.6f %.6f %6.6f %6.6f %6.6f %6.6f %6.6f %6.6f %6.6f %6.6f %6.6f %6.6f",
gCall.getLocation().getContig(),
gCall.getLocation().getStart(),
vc.getLocation().getContig(),
vc.getLocation().getStart(),
ref,
readDepth,
readCount,
maxMappingQual,
gCall.getBases(),
sb.toString(),
nextVrsRef,
nextVrsBest,
posteriors[0],
@ -111,32 +141,9 @@ public class GeliTextWriter implements GeliGenotypeWriter {
mWriter.flush(); // necessary so that writing to an output stream will work
}
/**
* add a no call to the genotype file, if supported.
*
* @param position the position to add the no call at
*/
public void addNoCall(int position) {
throw new UnsupportedOperationException("Geli text format doesn't support a no-call call.");
}
/** finish writing, closing any open files. */
public void close() {
mWriter.flush();
mWriter.close();
}
/**
* add a multi-sample call if we support it
*
* @param genotypes the list of genotypes
*/
public void addMultiSampleCall(List<Genotype> genotypes, VariationCall metadata) {
throw new UnsupportedOperationException("Geli text doesn't support multisample calls");
}
/** @return true if we support multisample, false otherwise */
public boolean supportsMultiSample() {
return false;
}
}

208
java/src/org/broadinstitute/sting/utils/genotype/glf/GLFGenotypeCall.java
View File

@ -1,208 +0,0 @@
package org.broadinstitute.sting.utils.genotype.glf;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.genotype.*;
import java.util.Arrays;
/**
* @author ebanks
* <p/>
* Class GLFGenotypeCall
* <p/>
* The implementation of the genotype interface, specific to GLF
*/
public class GLFGenotypeCall implements GenotypeCall, ReadBacked, LikelihoodsBacked {
private final char mRefBase;
private final GenomeLoc mLocation;
private ReadBackedPileup mPileup;
private double[] mLikelihoods;
private double mNegLog10PError;
private String mGenotype;
private Variation mVariation = null;
/**
* Generate a single sample genotype object, containing everything we need to represent calls out of a genotyper object
*
* @param ref the ref character
* @param loc the genome loc
*/
public GLFGenotypeCall(String ref, GenomeLoc loc) {
mRefBase = ref.charAt(0);
mGenotype = Utils.dupString(mRefBase, 2);
// fill in empty data
mLocation = loc;
mLikelihoods = new double[10];
Arrays.fill(mLikelihoods, Double.MIN_VALUE);
mPileup = null;
mNegLog10PError = Double.MIN_VALUE;
}
public void setPileup(ReadBackedPileup pileup) {
mPileup = pileup;
}
public void setLikelihoods(double[] likelihoods) {
mLikelihoods = likelihoods;
}
public void setNegLog10PError(double negLog10PError) {
mNegLog10PError = negLog10PError;
}
public void setVariation(Variation variation) {
this.mVariation = variation;
}
public void setGenotype(String genotype) {
mGenotype = genotype;
}
public void setGenotype(DiploidGenotype genotype) {
setGenotype(genotype.toString());
}
@Override
public boolean equals(Object other) {
if (other == null || !(other instanceof GLFGenotypeCall))
return false;
return (this.mRefBase == ((GLFGenotypeCall)other).mRefBase);
}
public String toString() {
return String.format("%s ref=%s depth=%d negLog10PError=%.2f",
getLocation(), mRefBase, getReadCount(), getNegLog10PError());
}
/**
* get the confidence we have
*
* @return get the one minus error value
*/
public double getNegLog10PError() {
return mNegLog10PError;
}
/**
* get the bases that represent this
*
* @return the bases, as a string
*/
public String getBases() {
return Character.toString(mRefBase);
}
/**
* get the ploidy
*
* @return the ploidy value
*/
public int getPloidy() {
return 2;
}
/**
* Returns true if both observed alleles are the same (regardless of whether they are ref or alt)
*
* @return true if we're homozygous, false otherwise
*/
public boolean isHom() {
return true;
}
/**
* Returns true if observed alleles differ (regardless of whether they are ref or alt)
*
* @return true if we're het, false otherwise
*/
public boolean isHet() {
return !isHom();
}
public boolean isNoCall() { return false; }
/**
*
* @return return this genotype as a variant
*/
public Variation toVariation(char ref) {
if ( mVariation == null ) {
mVariation = new BasicGenotypeBackedVariation(ref, mLocation, Variation.VARIANT_TYPE.REFERENCE);
}
return mVariation;
}
/**
* Location of this genotype on the reference (on the forward strand). If the allele is insertion/deletion, the first inserted/deleted base
* is located right <i>after</i> the specified location
*
* @return position on the genome wrapped in GenomeLoc object
*/
public GenomeLoc getLocation() {
return this.mLocation;
}
/**
* returns true if the genotype is a point genotype, false if it's a indel / deletion
*
* @return true is a SNP
*/
public boolean isPointGenotype() {
return true;
}
/**
* given the reference, are we a variant? (non-ref)
*
* @param ref the reference base or bases
*
* @return true if we're a variant
*/
public boolean isVariant(char ref) {
return !Utils.dupString(mRefBase, 2).equals(mGenotype);
}
/**
* get the pileup that backs this genotype call
*
* @return a pileup
*/
public ReadBackedPileup getPileup() {
return mPileup;
}
/**
* get the count of reads
*
* @return the number of reads we're backed by
*/
public int getReadCount() {
return (mPileup != null ? mPileup.getReads().size() : 0);
}
/**
* get the reference character
*
* @return the reference character
*/
public String getReference() {
return Character.toString(mRefBase);
}
/**
* get the posteriors
*
* @return the posteriors
*/
public double[] getLikelihoods() {
return mLikelihoods;
}
}

99
java/src/org/broadinstitute/sting/utils/genotype/glf/GLFWriter.java
View File

@ -1,17 +1,20 @@
package org.broadinstitute.sting.utils.genotype.glf;
import net.sf.samtools.SAMRecord;
import net.sf.samtools.SAMSequenceRecord;
import net.sf.samtools.util.BinaryCodec;
import net.sf.samtools.util.BlockCompressedOutputStream;
import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.GenomeLocParser;
import org.broadinstitute.sting.utils.genotype.*;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.genotype.LikelihoodObject;
import org.broadinstitute.sting.utils.genotype.IndelLikelihood;
import org.broadinstitute.sting.utils.genotype.CalledGenotype;
import org.broadinstitute.sting.gatk.contexts.variantcontext.*;
import java.io.DataOutputStream;
import java.io.File;
import java.io.OutputStream;
import java.util.List;
/*
* Copyright (c) 2009 The Broad Institute
*
@ -113,12 +116,12 @@ public class GLFWriter implements GLFGenotypeWriter {
* @param rmsMapQ the root mean square of the mapping quality
* @param lhValues the GenotypeLikelihoods object, representing the genotype likelyhoods
*/
public void addGenotypeCall(SAMSequenceRecord contig,
int genomicLoc,
float rmsMapQ,
char refBase,
int readDepth,
LikelihoodObject lhValues) {
public void addCall(SAMSequenceRecord contig,
int genomicLoc,
float rmsMapQ,
char refBase,
int readDepth,
LikelihoodObject lhValues) {
if ( headerText == null )
throw new IllegalStateException("The GLF Header must be written before calls can be added");
@ -137,45 +140,61 @@ public class GLFWriter implements GLFGenotypeWriter {
}
/**
* Add a genotype, given a genotype call
* Add a genotype, given a variant context
*
* @param call the genotype call
* @param vc the variant context representing the call to add
*/
public void addGenotypeCall(Genotype call) {
public void addCall(VariantContext vc) {
if ( headerText == null )
throw new IllegalStateException("The GLF Header must be written before calls can be added");
if ( !(call instanceof GLFGenotypeCall) )
throw new IllegalArgumentException("Only GLFGenotypeCall should be passed in to the GLF writers");
GLFGenotypeCall gCall = (GLFGenotypeCall) call;
char ref = gCall.getReference().charAt(0);
char ref = vc.getReference().toString().charAt(0);
if ( vc.getNSamples() != 1 )
throw new IllegalArgumentException("The GLF format does not support multi-sample or no-calls");
// get likelihood information if available
LikelihoodObject obj = new LikelihoodObject(gCall.getLikelihoods(), LikelihoodObject.LIKELIHOOD_TYPE.LOG);
org.broadinstitute.sting.gatk.contexts.variantcontext.Genotype genotype = vc.getGenotypes().values().iterator().next();
if ( genotype.isNoCall() )
throw new IllegalArgumentException("The GLF format does not support no-calls");
ReadBackedPileup pileup;
double[] likelihoods;
if ( genotype instanceof CalledGenotype) {
pileup = ((CalledGenotype)genotype).getReadBackedPileup();
likelihoods = ((CalledGenotype)genotype).getLikelihoods();
} else {
pileup = (ReadBackedPileup)genotype.getAttribute(CalledGenotype.READBACKEDPILEUP_ATTRIBUTE_KEY);
likelihoods = (double[])genotype.getAttribute(CalledGenotype.LIKELIHOODS_ATTRIBUTE_KEY);
}
if ( likelihoods == null )
throw new IllegalArgumentException("The GLF format requires likelihoods");
LikelihoodObject obj = new LikelihoodObject(likelihoods, LikelihoodObject.LIKELIHOOD_TYPE.LOG);
obj.setLikelihoodType(LikelihoodObject.LIKELIHOOD_TYPE.NEGATIVE_LOG); // transform! ... to negitive log likelihoods
// calculate the RMS mapping qualities and the read depth
double rms = 0.0;
if (gCall.getPileup() != null)
rms = calculateRMS(gCall.getPileup().getReads());
int readCount = gCall.getReadCount();
this.addGenotypeCall(GenomeLocParser.getContigInfo(gCall.getLocation().getContig()), (int) gCall.getLocation().getStart(), (float) rms, ref, readCount, obj);
int readCount = 0;
if ( pileup != null ) {
rms = calculateRMS(pileup);
readCount = pileup.size();
}
addCall(GenomeLocParser.getContigInfo(vc.getLocation().getContig()), (int)vc.getLocation().getStart(), (float) rms, ref, readCount, obj);
}
/**
* calculate the rms , given the read pileup
*
* @param reads the read array
* @param pileup the pileup
*
* @return the rms of the read mapping qualities
*/
private double calculateRMS(List<SAMRecord> reads) {
int[] qualities = new int[reads.size()];
for (int i = 0; i < reads.size(); i++) {
qualities[i] = reads.get(i).getMappingQuality();
}
private double calculateRMS(ReadBackedPileup pileup) {
int[] qualities = new int[pileup.size()];
int index = 0;
for (PileupElement p : pileup )
qualities[index++] = p.getMappingQual();
return MathUtils.rms(qualities);
}
@ -224,16 +243,6 @@ public class GLFWriter implements GLFGenotypeWriter {
mLastRecord = call;
}
/**
* add a no call to the genotype file, if supported.
*
* @param position the position
*/
public void addNoCall(int position) {
// glf doesn't support this operation
throw new UnsupportedOperationException("GLF doesn't support a 'no call' call.");
}
/**
* add a GLF record to the output file
*
@ -295,20 +304,6 @@ public class GLFWriter implements GLFGenotypeWriter {
writeEndRecord();
outputBinaryCodec.close();
}
/**
* add a multi-sample call if we support it
*
* @param genotypes the list of genotypes
*/
public void addMultiSampleCall(List<Genotype> genotypes, VariationCall metadata) {
throw new UnsupportedOperationException("GLF writer doesn't support multisample calls");
}
/** @return true if we support multisample, false otherwise */
public boolean supportsMultiSample() {
return false;
}
}

80
java/src/org/broadinstitute/sting/utils/genotype/tabular/TabularLFWriter.java
View File

@ -1,13 +1,14 @@
package org.broadinstitute.sting.utils.genotype.tabular;
import org.broadinstitute.sting.utils.StingException;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.genotype.*;
import org.broadinstitute.sting.gatk.contexts.variantcontext.*;
import java.io.File;
import java.io.FileNotFoundException;
import java.io.PrintStream;
import java.util.Arrays;
import java.util.List;
/**
@ -26,7 +27,7 @@ public class TabularLFWriter implements GenotypeWriter {
/**
* construct, writing to a specified file
* @param writeTo
* @param writeTo file to write to
*/
public TabularLFWriter(File writeTo) {
try {
@ -39,52 +40,56 @@ public class TabularLFWriter implements GenotypeWriter {
}
/**
* Add a genotype, given a genotype locus
* Add a genotype, given a variant context
*
* @param locus the locus to add
* @param vc the variant context representing the call to add
*/
public void addGenotypeCall(Genotype locus) {
double likelihoods[];
int readDepth = -1;
double nextVrsBest = 0;
double nextVrsRef = 0;
if (!(locus instanceof LikelihoodsBacked)) {
public void addCall(VariantContext vc) {
if ( vc.getNSamples() != 1 )
throw new IllegalArgumentException("The tabular LF format does not support multi-sample or no-calls");
org.broadinstitute.sting.gatk.contexts.variantcontext.Genotype genotype = vc.getGenotypes().values().iterator().next();
if ( genotype.isNoCall() )
throw new IllegalArgumentException("The tabular LF format does not support no-calls");
ReadBackedPileup pileup;
double[] likelihoods;
if ( genotype instanceof CalledGenotype) {
pileup = ((CalledGenotype)genotype).getReadBackedPileup();
likelihoods = ((CalledGenotype)genotype).getLikelihoods();
} else {
pileup = (ReadBackedPileup)genotype.getAttribute(CalledGenotype.READBACKEDPILEUP_ATTRIBUTE_KEY);
likelihoods = (double[])genotype.getAttribute(CalledGenotype.LIKELIHOODS_ATTRIBUTE_KEY);
}
if ( likelihoods == null ) {
likelihoods = new double[10];
Arrays.fill(likelihoods, Double.MIN_VALUE);
} else {
likelihoods = ((LikelihoodsBacked) locus).getLikelihoods();
}
char ref = locus.getReference().charAt(0);
if (locus instanceof ReadBacked) {
readDepth = ((ReadBacked)locus).getReadCount();
}
int readDepth = pileup == null ? -1 : pileup.size();
double nextVrsBest = 0;
double nextVrsRef = 0;
char ref = vc.getReference().toString().charAt(0);
/**
* This output is not correct, but I don't we even use this format anymore. If we do, someone
* should change this code
*/
outStream.println(String.format("%s %s %c %s %s %f %f %f %f %d %s",
locus.getLocation().toString(),
vc.getLocation().toString(),
"NOT OUTPUTED",
ref,
locus.getBases(),
locus.getBases(),
genotype.getGenotypeString(),
genotype.getGenotypeString(),
-1,
-1,
nextVrsRef,
nextVrsBest,
readDepth,
locus.getBases()));
}
/**
* add a no call to the genotype file, if supported.
*
* @param position
*/
public void addNoCall(int position) {
throw new StingException("TabularLFWriter doesn't support no-calls");
genotype.getGenotypeString()));
}
/** finish writing, closing any open files. */
@ -93,19 +98,4 @@ public class TabularLFWriter implements GenotypeWriter {
outStream.close();
}
}
/**
* add a multi-sample call if we support it
*
* @param genotypes the list of genotypes, that are backed by sample information
*/
public void addMultiSampleCall(List<Genotype> genotypes, VariationCall metadata) {
throw new UnsupportedOperationException("Tabular LF doesn't support multisample calls");
}
/** @return true if we support multisample, false otherwise */
public boolean supportsMultiSample() {
return false;
}
}

226
java/src/org/broadinstitute/sting/utils/genotype/vcf/VCFGenotypeCall.java
View File

@ -1,226 +0,0 @@
package org.broadinstitute.sting.utils.genotype.vcf;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.genotype.*;
/**
* @author ebanks
* <p/>
* Class VCFGenotypeCall
* <p/>
* The implementation of the genotype interface, specific to VCF
*/
public class VCFGenotypeCall extends AlleleConstrainedGenotype implements GenotypeCall, ReadBacked, SampleBacked, Cloneable {
private final String mRef;
private final GenomeLoc mLocation;
private ReadBackedPileup mPileup = null;
private int mCoverage = 0;
private double mNegLog10PError = -1;
private Variation mVariation = null;
// the best genotype
private DiploidGenotype mGenotype = null;
// the sample name, used to propulate the SampleBacked interface
private String mSampleName;
public VCFGenotypeCall(String ref, GenomeLoc loc) {
super(ref);
mRef = ref;
mLocation = loc;
// fill in empty data
mGenotype = DiploidGenotype.createHomGenotype(ref.charAt(0));
mSampleName = "";
}
public VCFGenotypeCall(String ref, GenomeLoc loc, DiploidGenotype genotype, double negLog10PError, int coverage, String sample) {
super(ref);
mRef = ref;
mLocation = loc;
mGenotype = genotype;
mNegLog10PError = negLog10PError;
mCoverage = coverage;
mSampleName = sample;
}
public void setPileup(ReadBackedPileup pileup) {
mPileup = pileup;
}
public void setGenotype(DiploidGenotype genotype) {
mGenotype = genotype;
}
public void setNegLog10PError(double negLog10PError) {
mNegLog10PError = negLog10PError;
}
public void setVariation(Variation variation) {
this.mVariation = variation;
}
public void setSampleName(String name) {
mSampleName = name;
}
@Override
public boolean equals(Object other) {
if ( other == null || !(other instanceof VCFGenotypeCall) )
return false;
VCFGenotypeCall otherCall = (VCFGenotypeCall) other;
return mGenotype.equals(otherCall.mGenotype) &&
mNegLog10PError == otherCall.mNegLog10PError &&
mLocation.equals(otherCall.mLocation) &&
mRef.equals(otherCall.mRef);
}
public String toString() {
return String.format("%s best=%s ref=%s depth=%d negLog10PError=%.2f",
getLocation(), mGenotype, mRef, getReadCount(), getNegLog10PError());
}
/**
* get the confidence we have
*
* @return get the one minus error value
*/
public double getNegLog10PError() {
return mNegLog10PError;
}
// get the best genotype
protected DiploidGenotype getBestGenotype() {
return mGenotype;
}
/**
* get the ploidy
*
* @return the ploidy value
*/
public int getPloidy() {
return 2;
}
/**
* Returns true if both observed alleles are the same (regardless of whether they are ref or alt)
*
* @return true if we're homozygous, false otherwise
*/
public boolean isHom() {
return getBestGenotype().isHom();
}
/**
* Returns true if observed alleles differ (regardless of whether they are ref or alt)
*
* @return true if we're het, false otherwise
*/
public boolean isHet() {
return !isHom();
}
// You can't make a 'no call' genotype call
public boolean isNoCall() { return false; }
/**
* Location of this genotype on the reference (on the forward strand). If the allele is insertion/deletion, the first inserted/deleted base
* is located right <i>after</i> the specified location
*
* @return position on the genome wrapped in GenomeLoc object
*/
public GenomeLoc getLocation() {
return this.mLocation;
}
/**
* returns true if the genotype is a point genotype, false if it's a indel / deletion
*
* @return true is a SNP
*/
public boolean isPointGenotype() {
return true;
}
/**
* given the reference, are we a variant? (non-ref)
*
* @param ref the reference base or bases
*
* @return true if we're a variant
*/
public boolean isVariant(char ref) {
return !getBestGenotype().isHomRef(ref);
}
/**
* are we a variant? (non-ref)
*
* @return true if we're a variant
*/
public boolean isVariant() {
return isVariant(mRef.charAt(0));
}
/**
*
* @return return this genotype as a variant
*/
public Variation toVariation(char ref) {
if ( mVariation == null ) {
VCFVariationCall var = new VCFVariationCall(ref, mLocation, isVariant() ? Variation.VARIANT_TYPE.SNP : Variation.VARIANT_TYPE.REFERENCE);
var.setConfidence(10 * mNegLog10PError);
if ( !mGenotype.isHomRef(ref) ) {
if ( mGenotype.base1 != ref )
var.addAlternateAllele(Character.toString(mGenotype.base1));
if ( mGenotype.isHet() && mGenotype.base2 != ref )
var.addAlternateAllele(Character.toString(mGenotype.base2));
}
mVariation = var;
}
return mVariation;
}
/**
* get the pileup that backs this genotype call
*
* @return a pileup
*/
public ReadBackedPileup getPileup() {
return mPileup;
}
/**
* get the count of reads
*
* @return the number of reads we're backed by
*/
public int getReadCount() {
return (mCoverage > 0 ? mCoverage : (mPileup != null ? mPileup.size() : 0));
}
/**
* get the reference string
*
* @return the reference string
*/
public String getReference() {
return mRef;
}
/**
* @return returns the sample name for this genotype
*/
public String getSampleName() {
return mSampleName;
}
}

11
java/src/org/broadinstitute/sting/utils/genotype/vcf/VCFGenotypeRecord.java
View File

@ -1,7 +1,6 @@
package org.broadinstitute.sting.utils.genotype.vcf;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.genotype.*;
import java.util.*;
@ -12,7 +11,7 @@ import java.util.*;
* Class VCFGenotypeRecord
* <p/>
*/
public class VCFGenotypeRecord implements Genotype, SampleBacked {
public class VCFGenotypeRecord {
// key names
public static final String GENOTYPE_KEY = "GT";
@ -143,10 +142,6 @@ public class VCFGenotypeRecord implements Genotype, SampleBacked {
return mRecord != null ? mRecord.getReference() : "N";
}
public Variation toVariation(char ref) {
return mRecord != null ? mRecord : null;
}
public String getBases() {
String genotype = "";
for ( VCFGenotypeEncoding encoding : mGenotypeAlleles )
@ -201,6 +196,10 @@ public class VCFGenotypeRecord implements Genotype, SampleBacked {
return 2;
}
public VCFRecord getRecord() {
return mRecord;
}
private String toGenotypeString(List<VCFGenotypeEncoding> altAlleles) {
String str = "";
boolean first = true;

191
java/src/org/broadinstitute/sting/utils/genotype/vcf/VCFGenotypeWriterAdapter.java
View File

@ -1,6 +1,8 @@
package org.broadinstitute.sting.utils.genotype.vcf;
import org.broadinstitute.sting.utils.genotype.*;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.Allele;
import org.broadinstitute.sting.gatk.refdata.VariantContextAdaptors;
import org.apache.log4j.Logger;
import java.io.File;
@ -28,6 +30,10 @@ public class VCFGenotypeWriterAdapter implements VCFGenotypeWriter {
// validation stringency
private VALIDATION_STRINGENCY validationStringency = VALIDATION_STRINGENCY.STRICT;
// standard genotype format strings
private static String[] standardGenotypeFormatStrings = { VCFGenotypeRecord.GENOTYPE_KEY,
VCFGenotypeRecord.DEPTH_KEY,
VCFGenotypeRecord.GENOTYPE_QUALITY_KEY };
public VCFGenotypeWriterAdapter(File writeTo) {
if (writeTo == null) throw new RuntimeException("VCF output file must not be null");
@ -58,188 +64,47 @@ public class VCFGenotypeWriterAdapter implements VCFGenotypeWriter {
mWriter.writeHeader(mHeader);
}
/**
* Add a genotype, given a genotype locus
*
* @param call the locus to add
*/
public void addGenotypeCall(Genotype call) {
throw new UnsupportedOperationException("VCF calls require locusdata; use the addMultiSampleCall method instead");
}
/**
* add a no call to the genotype file, if supported.
*
* @param position the position to add the no call at
*/
public void addNoCall(int position) {
throw new UnsupportedOperationException("We don't currently support no-calls in VCF");
}
/** finish writing, closing any open files. */
public void close() {
mWriter.close();
}
/**
* add a multi-sample call if we support it
* Add a genotype, given a variant context
*
* @param genotypes the list of genotypes
* @param vc the variant context representing the call to add
*/
public void addMultiSampleCall(List<Genotype> genotypes, VariationCall locusdata) {
public void addCall(VariantContext vc) {
if ( mHeader == null )
throw new IllegalStateException("The VCF Header must be written before records can be added");
if ( locusdata != null && !(locusdata instanceof VCFVariationCall) )
throw new IllegalArgumentException("Only VCFVariationCall objects should be passed in to the VCF writers");
List<String> formatStrings;
if ( vc.getChromosomeCount() > 0 )
formatStrings = Arrays.asList(standardGenotypeFormatStrings);
else
formatStrings = new ArrayList<String>();
VCFRecord call = VariantContextAdaptors.toVCF(vc, vc.getReference().toString().charAt(0), formatStrings, false);
VCFParameters params = new VCFParameters();
if ( genotypes.size() > 0 )
params.addFormatItem(VCFGenotypeRecord.GENOTYPE_KEY);
// get the location and reference
if ( genotypes.size() == 0 ) {
if ( locusdata == null )
throw new IllegalArgumentException("Unable to parse out the current location: genotype array must contain at least one entry or have variation data");
params.setLocations(locusdata.getLocation(), locusdata.getReference());
// if there is no genotype data, we'll also need to set an alternate allele
if ( locusdata.isBiallelic() && locusdata.isSNP() )
params.addAlternateBase(new VCFGenotypeEncoding(locusdata.getAlternateAlleleList().get(0)));
} else {
params.setLocations(genotypes.get(0).getLocation(), genotypes.get(0).getReference());
Set<Allele> altAlleles = vc.getAlternateAlleles();
StringBuffer altAlleleCountString = new StringBuffer();
for ( Allele allele : altAlleles ) {
if ( altAlleleCountString.length() > 0 )
altAlleleCountString.append(",");
altAlleleCountString.append(vc.getChromosomeCount(allele));
}
if ( vc.getChromosomeCount() > 0 ) {
call.addInfoField(VCFRecord.ALLELE_NUMBER_KEY, String.format("%d", vc.getChromosomeCount()));
if ( altAlleleCountString.length() > 0 )
call.addInfoField(VCFRecord.ALLELE_COUNT_KEY, altAlleleCountString.toString());
}
Map<String, VCFGenotypeCall> genotypeMap = genotypeListToSampleNameMap(genotypes);
int totalAlleles = 0;
for (String name : mHeader.getGenotypeSamples()) {
if (genotypeMap.containsKey(name)) {
Genotype gtype = genotypeMap.get(name);
VCFGenotypeRecord record = VCFUtils.createVCFGenotypeRecord(params, (VCFGenotypeCall)gtype);
params.addGenotypeRecord(record);
totalAlleles += record.getAlleles().size();
genotypeMap.remove(name);
} else {
VCFGenotypeRecord record = createNoCallRecord(name);
params.addGenotypeRecord(record);
}
}
if ( validationStringency == VALIDATION_STRINGENCY.STRICT && genotypeMap.size() > 0 ) {
for (String name : genotypeMap.keySet())
logger.fatal("Genotype " + name + " was present in the VCFHeader");
throw new IllegalArgumentException("Genotype array passed to VCFGenotypeWriterAdapter contained Genotypes not in the VCF header");
}
// info fields
Map<String, String> infoFields = getInfoFields((VCFVariationCall)locusdata);
if ( totalAlleles > 0 ) {
infoFields.put(VCFRecord.ALLELE_NUMBER_KEY, String.format("%d", totalAlleles));
// count up the alternate counts
List<Integer> altAlleleCounts = params.getAlleleCounts();
if ( altAlleleCounts.size() > 0 ) {
StringBuffer sb = new StringBuffer();
sb.append(altAlleleCounts.get(0));
for (int i = 1; i < altAlleleCounts.size(); i++ ) {
sb.append(",");
sb.append(altAlleleCounts.get(i));
}
infoFields.put(VCFRecord.ALLELE_COUNT_KEY, sb.toString());
}
}
// q-score
double qual = (locusdata == null) ? 0 : ((VCFVariationCall)locusdata).getConfidence();
// min Q-score is zero
qual = Math.max(qual, 0);
// dbsnp id
String dbSnpID = null;
if ( locusdata != null )
dbSnpID = ((VCFVariationCall)locusdata).getID();
VCFRecord vcfRecord = new VCFRecord(params.getReferenceBases(),
params.getContig(),
params.getPosition(),
(dbSnpID == null ? VCFRecord.EMPTY_ID_FIELD : dbSnpID),
params.getAlternateBases(),
qual,
VCFRecord.UNFILTERED,
infoFields,
params.getFormatString(),
params.getGenotypeRecords());
mWriter.addRecord(vcfRecord, validationStringency);
mWriter.addRecord(call, validationStringency);
}
public void addRecord(VCFRecord vcfRecord) {
mWriter.addRecord(vcfRecord, validationStringency);
}
/**
* get the information fields of the VCF record, given the meta data and parameters
*
* @param locusdata the metadata associated with this multi sample call
*
* @return a mapping of info field to value
*/
private static Map<String, String> getInfoFields(VCFVariationCall locusdata) {
Map<String, String> infoFields = new HashMap<String, String>();
if ( locusdata != null ) {
if ( locusdata.getSLOD() != null )
infoFields.put(VCFRecord.STRAND_BIAS_KEY, String.format("%.2f", locusdata.getSLOD()));
if ( locusdata.hasNonRefAlleleFrequency() )
infoFields.put(VCFRecord.ALLELE_FREQUENCY_KEY, String.format("%.2f", locusdata.getNonRefAlleleFrequency()));
Map<String, String> otherFields = locusdata.getFields();
if ( otherFields != null ) {
infoFields.putAll(otherFields);
}
}
return infoFields;
}
/**
* create a no call record
*
* @param sampleName the sample name
*
* @return a VCFGenotypeRecord for the no call situation
*/
private VCFGenotypeRecord createNoCallRecord(String sampleName) {
List<VCFGenotypeEncoding> alleles = new ArrayList<VCFGenotypeEncoding>();
alleles.add(new VCFGenotypeEncoding(VCFGenotypeRecord.EMPTY_ALLELE));
alleles.add(new VCFGenotypeEncoding(VCFGenotypeRecord.EMPTY_ALLELE));
return new VCFGenotypeRecord(sampleName, alleles, VCFGenotypeRecord.PHASE.UNPHASED);
}
/** @return true if we support multisample, false otherwise */
public boolean supportsMultiSample() {
return true;
}
/**
* create a genotype mapping from a list and their sample names
* while we're at it, checks that all genotypes are VCF-based
*
* @param list a list of genotype samples
*
* @return a mapping of the sample name to genotype fields
*/
private static Map<String, VCFGenotypeCall> genotypeListToSampleNameMap(List<Genotype> list) {
Map<String, VCFGenotypeCall> map = new HashMap<String, VCFGenotypeCall>();
for (Genotype rec : list) {
if ( !(rec instanceof VCFGenotypeCall) )
throw new IllegalArgumentException("Only VCFGenotypeCalls should be passed in to the VCF writers");
map.put(((VCFGenotypeCall) rec).getSampleName(), (VCFGenotypeCall) rec);
}
return map;
}
/**
* set the validation stringency
* @param value validation stringency value

62
java/src/org/broadinstitute/sting/utils/genotype/vcf/VCFRecord.java
View File

@ -9,7 +9,7 @@ import java.util.*;
/**
* the basic VCF record type
*/
public class VCFRecord implements Variation, VariantBackedByGenotype {
public class VCFRecord {
// standard info field keys
public static final String ANCESTRAL_ALLELE_KEY = "AA";
@ -58,7 +58,7 @@ public class VCFRecord implements Variation, VariantBackedByGenotype {
private String mGenotypeFormatString;
// our genotype sample fields
private final List<Genotype> mGenotypeRecords = new ArrayList<Genotype>();
private final List<VCFGenotypeRecord> mGenotypeRecords = new ArrayList<VCFGenotypeRecord>();
/**
* given a reference base, a location, and the format string, create a VCF record.
@ -274,9 +274,9 @@ public class VCFRecord implements Variation, VariantBackedByGenotype {
return 0.0;
}
public VARIANT_TYPE getType() {
public Variation.VARIANT_TYPE getType() {
if ( ! hasAlternateAllele() )
return VARIANT_TYPE.REFERENCE;
return Variation.VARIANT_TYPE.REFERENCE;
VCFGenotypeEncoding.TYPE type = mAlts.get(0).getType();
for (int i = 1; i < mAlts.size(); i++) {
@ -286,25 +286,25 @@ public class VCFRecord implements Variation, VariantBackedByGenotype {
switch ( type ) {
case SINGLE_BASE:
return VARIANT_TYPE.SNP;
return Variation.VARIANT_TYPE.SNP;
case UNCALLED:
// If there are no alt alleles, all of the genotypes are reference or no calls, so we're a reference site
return VARIANT_TYPE.REFERENCE;
return Variation.VARIANT_TYPE.REFERENCE;
case DELETION:
return VARIANT_TYPE.DELETION;
return Variation.VARIANT_TYPE.DELETION;
case INSERTION:
return VARIANT_TYPE.INSERTION;
return Variation.VARIANT_TYPE.INSERTION;
}
throw new IllegalStateException("The record contains unknown genotype encodings");
}
public boolean isDeletion() {
return getType() == VARIANT_TYPE.DELETION;
return getType() == Variation.VARIANT_TYPE.DELETION;
}
public boolean isInsertion() {
return getType() == VARIANT_TYPE.INSERTION;
return getType() == Variation.VARIANT_TYPE.INSERTION;
}
public boolean isIndel() {
@ -312,7 +312,7 @@ public class VCFRecord implements Variation, VariantBackedByGenotype {
}
public boolean isSNP() {
return getType() == VARIANT_TYPE.SNP;
return getType() == Variation.VARIANT_TYPE.SNP;
}
public boolean isNovel() {
@ -385,7 +385,7 @@ public class VCFRecord implements Variation, VariantBackedByGenotype {
*
* @return a map, of the info key-value pairs
*/
public Map<String, String> getInfoValues() {
public final Map<String, String> getInfoValues() {
if (mInfoFields.size() < 1) {
Map<String, String> map = new HashMap<String, String>();
map.put(".", "");
@ -395,43 +395,16 @@ public class VCFRecord implements Variation, VariantBackedByGenotype {
}
public List<VCFGenotypeRecord> getVCFGenotypeRecords() {
ArrayList<VCFGenotypeRecord> list = new ArrayList<VCFGenotypeRecord>();
for ( Genotype g : mGenotypeRecords )
list.add((VCFGenotypeRecord)g);
return list;
}
public List<Genotype> getGenotypes() {
return mGenotypeRecords;
}
public Genotype getCalledGenotype() {
if ( mGenotypeRecords == null || mGenotypeRecords.size() != 1 )
throw new IllegalStateException("There is not one and only one genotype associated with this record");
VCFGenotypeRecord record = (VCFGenotypeRecord)mGenotypeRecords.get(0);
if ( record.isEmptyGenotype() )
return null;
return record;
}
public boolean hasGenotype(DiploidGenotype x) {
if ( mGenotypeRecords == null )
return false;
for ( Genotype g : mGenotypeRecords ) {
if ( DiploidGenotype.valueOf(g.getBases()).equals(x) )
return true;
}
return false;
}
/**
* @return a List of the sample names
*/
public String[] getSampleNames() {
String names[] = new String[mGenotypeRecords.size()];
for (int i = 0; i < mGenotypeRecords.size(); i++) {
VCFGenotypeRecord rec = (VCFGenotypeRecord)mGenotypeRecords.get(i);
names[i] = rec.getSampleName();
names[i] = mGenotypeRecords.get(i).getSampleName();
}
return names;
}
@ -611,6 +584,7 @@ public class VCFRecord implements Variation, VariantBackedByGenotype {
protected String createInfoString() {
String info = "";
for (String str : getInfoValues().keySet()) {
if (str.equals(EMPTY_INFO_FIELD))
return EMPTY_INFO_FIELD;
else
@ -626,10 +600,10 @@ public class VCFRecord implements Variation, VariantBackedByGenotype {
* @param header the header object
*/
private void addGenotypeData(StringBuilder builder, VCFHeader header) {
Map<String, VCFGenotypeRecord> gMap = genotypeListToMap(getGenotypes());
Map<String, VCFGenotypeRecord> gMap = genotypeListToMap(getVCFGenotypeRecords());
StringBuffer tempStr = new StringBuffer();
if ( header.getGenotypeSamples().size() < getGenotypes().size() ) {
if ( header.getGenotypeSamples().size() < getVCFGenotypeRecords().size() ) {
for ( String sample : gMap.keySet() ) {
if ( !header.getGenotypeSamples().contains(sample) )
System.err.println("Sample " + sample + " is a duplicate or is otherwise not present in the header");
@ -687,10 +661,10 @@ public class VCFRecord implements Variation, VariantBackedByGenotype {
* @param list a list of genotype samples
* @return a mapping of the sample name to VCF genotype record
*/
private static Map<String, VCFGenotypeRecord> genotypeListToMap(List<Genotype> list) {
private static Map<String, VCFGenotypeRecord> genotypeListToMap(List<VCFGenotypeRecord> list) {
Map<String, VCFGenotypeRecord> map = new HashMap<String, VCFGenotypeRecord>();
for (int i = 0; i < list.size(); i++) {
VCFGenotypeRecord rec = (VCFGenotypeRecord)list.get(i);
VCFGenotypeRecord rec = list.get(i);
map.put(rec.getSampleName(), rec);
}
return map;

33
java/src/org/broadinstitute/sting/utils/genotype/vcf/VCFUtils.java
View File

@ -6,7 +6,6 @@ import org.broadinstitute.sting.gatk.refdata.RodVCF;
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.utils.Pair;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.genotype.*;
import java.util.*;
@ -168,36 +167,6 @@ public class VCFUtils {
return record;
}
/**
* create the VCF genotype record
*
* @param params the VCF parameters object
* @param gtype the genotype
*
* @return a VCFGenotypeRecord
*/
public static VCFGenotypeRecord createVCFGenotypeRecord(VCFParameters params, VCFGenotypeCall gtype) {
List<VCFGenotypeEncoding> alleles = createAlleleArray(gtype);
for (VCFGenotypeEncoding allele : alleles) {
params.addAlternateBase(allele);
}
VCFGenotypeRecord record = new VCFGenotypeRecord(gtype.getSampleName(), alleles, VCFGenotypeRecord.PHASE.UNPHASED);
// calculate the RMS mapping qualities and the read depth
record.setField(VCFGenotypeRecord.DEPTH_KEY, String.valueOf(gtype.getReadCount()));
params.addFormatItem(VCFGenotypeRecord.DEPTH_KEY);
double qual = Math.min(10.0 * gtype.getNegLog10PError(), VCFGenotypeRecord.MAX_QUAL_VALUE);
if ( qual >= 0 )
record.setField(VCFGenotypeRecord.GENOTYPE_QUALITY_KEY, String.format("%.2f", qual));
else
record.setField(VCFGenotypeRecord.GENOTYPE_QUALITY_KEY, String.format("%d", VCFGenotypeRecord.MISSING_GENOTYPE_QUALITY));
params.addFormatItem(VCFGenotypeRecord.GENOTYPE_QUALITY_KEY);
return record;
}
/**
* create the allele array?
*
@ -205,7 +174,7 @@ public class VCFUtils {
*
* @return a list of string representing the string array of alleles
*/
private static List<VCFGenotypeEncoding> createAlleleArray(Genotype gtype) {
private static List<VCFGenotypeEncoding> createAlleleArray(VCFGenotypeRecord gtype) {
List<VCFGenotypeEncoding> alleles = new ArrayList<VCFGenotypeEncoding>();
for (char allele : gtype.getBases().toCharArray()) {
alleles.add(new VCFGenotypeEncoding(String.valueOf(allele)));

250
java/src/org/broadinstitute/sting/utils/genotype/vcf/VCFVariationCall.java
View File

@ -1,250 +0,0 @@
package org.broadinstitute.sting.utils.genotype.vcf;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.genotype.*;
import java.util.*;
/**
* @author ebanks
* <p/>
* Class VCFVariationCall
* <p/>
* represents a VCF Variation
*/
public class VCFVariationCall implements VariationCall, VariantBackedByGenotype, ConfidenceBacked, SLODBacked, IDBacked, ArbitraryFieldsBacked {
// the discovery lod score
private double mConfidence = 0.0;
// the strand score lod
private Double mSLOD = null;
// the allele frequency field
private double mAlleleFrequency = -1.0;
// the location
private GenomeLoc mLoc;
// the ref base and alt bases
private char mRefBase;
private List<String> mAltBases = new ArrayList<String>();
// the variant type
private VARIANT_TYPE mType = VARIANT_TYPE.SNP;
// the id
private String mID;
// the genotypes
private List<Genotype> mGenotypes = null;
// the various info field values
private Map<String, String> mInfoFields;
/**
* create a VCF Variation object, given the following fields
*
* @param ref the reference base
* @param loc the locus
* @param type the variant type
*/
public VCFVariationCall(char ref, GenomeLoc loc, VARIANT_TYPE type) {
mRefBase = ref;
mLoc = loc;
mType = type;
}
/**
* get the reference base.
* @return a character, representing the reference base
*/
public String getReference() {
return String.valueOf(mRefBase);
}
/**
* get the genotype's location
*
* @return a GenomeLoc representing the location
*/
public GenomeLoc getLocation() {
return mLoc;
}
public boolean isBiallelic() {
return mAltBases.size() == 1;
}
public boolean isSNP() {
return mType == VARIANT_TYPE.SNP;
}
public boolean isInsertion() {
return mType == VARIANT_TYPE.INSERTION;
}
public boolean isIndel() {
return mType == VARIANT_TYPE.INSERTION || mType == VARIANT_TYPE.DELETION;
}
public boolean isDeletion() {
return mType == VARIANT_TYPE.DELETION;
}
public boolean isReference() {
return mType == VARIANT_TYPE.REFERENCE;
}
public VARIANT_TYPE getType() {
return mType;
}
public boolean hasNonRefAlleleFrequency() {
return mAlleleFrequency >= 0.0;
}
public double getNonRefAlleleFrequency() {
return mAlleleFrequency;
}
public double getNegLog10PError() {
return mConfidence / 10.0;
}
public List<String> getAlternateAlleleList() {
return mAltBases;
}
public void addAlternateAllele(String alt) {
mAltBases.add(alt);
}
public List<String> getAlleleList() {
LinkedList<String> alleles = new LinkedList<String>(mAltBases);
alleles.addFirst(getReference());
return alleles;
}
public char getAlternativeBaseForSNP() {
if ( !isSNP() )
throw new IllegalStateException("This variant is not a SNP");
if ( mAltBases.size() == 0 )
throw new IllegalStateException("No alternate alleles have been set");
return mAltBases.get(0).charAt(0);
}
public char getReferenceForSNP() {
if ( !isSNP() )
throw new IllegalStateException("This variant is not a SNP");
return mRefBase;
}
/**
* get the confidence
*
* @return the confidence
*/
public double getConfidence() {
return mConfidence;
}
/**
*
* @param confidence the confidence for this genotype
*/
public void setConfidence(double confidence) {
mConfidence = confidence;
}
/**
* get the strand lod
*
* @return the strand lod
*/
public Double getSLOD() {
return mSLOD;
}
/**
*
* @param slod the strand lod for this genotype
*/
public void setSLOD(double slod) {
mSLOD = slod;
}
/**
*
* @param frequency the allele frequency for this genotype
*/
public void setNonRefAlleleFrequency(double frequency) {
mAlleleFrequency = frequency;
}
/**
* @return returns the dbsnp id for this genotype
*/
public String getID() {
return mID;
}
public void setID(String id) {
mID = id;
}
/**
*
* @param calls the GenotypeCalls for this variation
*/
public void setGenotypeCalls(List<Genotype> calls) {
mGenotypes = calls;
}
/**
* @return a specific genotype that represents the called genotype
*/
public Genotype getCalledGenotype() {
if ( mGenotypes == null || mGenotypes.size() != 1 )
throw new IllegalStateException("There is not one and only one Genotype associated with this Variation");
return mGenotypes.get(0);
}
/**
* @return a list of all the genotypes
*/
public List<Genotype> getGenotypes() {
return mGenotypes;
}
/**
* do we have the specified genotype? not all backedByGenotypes
* have all the genotype data.
*
* @param x the genotype
*
* @return true if available, false otherwise
*/
public boolean hasGenotype(DiploidGenotype x) {
if ( mGenotypes == null )
return false;
for ( Genotype g : mGenotypes ) {
if ( DiploidGenotype.valueOf(g.getBases()).equals(x) )
return true;
}
return false;
}
/**
* @return returns te arbitrary info fields
*/
public Map<String, String> getFields() {
return mInfoFields;
}
public void setFields(Map<String, String> fields) {
mInfoFields = new HashMap<String, String>(fields);
}
}

22
java/test/org/broadinstitute/sting/gatk/refdata/RodVCFTest.java
View File

@ -5,11 +5,10 @@ import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.GenomeLocParser;
import org.broadinstitute.sting.utils.StingException;
import org.broadinstitute.sting.utils.fasta.IndexedFastaSequenceFile;
import org.broadinstitute.sting.utils.genotype.DiploidGenotype;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.genotype.vcf.VCFHeader;
import org.broadinstitute.sting.utils.genotype.vcf.VCFReader;
import org.broadinstitute.sting.utils.genotype.vcf.VCFRecord;
import org.broadinstitute.sting.utils.genotype.vcf.VCFGenotypeRecord;
import org.junit.Assert;
import static org.junit.Assert.fail;
import org.junit.BeforeClass;
@ -31,11 +30,11 @@ import java.util.List;
*/
public class RodVCFTest extends BaseTest {
private static IndexedFastaSequenceFile seq;
private static File vcfFile = new File(validationDataLocation + "vcfexample.vcf");
private VCFHeader mHeader;
@BeforeClass
public static void beforeTests() {
IndexedFastaSequenceFile seq;
try {
seq = new IndexedFastaSequenceFile(new File(oneKGLocation + "reference/human_b36_both.fasta"));
} catch (FileNotFoundException e) {
@ -147,13 +146,13 @@ public class RodVCFTest extends BaseTest {
Iterator<RodVCF> iter = vcf.createIterator("VCF", vcfFile);
RodVCF rec = iter.next();
// we should get back a ref 'G' and an alt 'A'
List<Genotype> list = rec.getGenotypes();
List<VCFGenotypeRecord> list = rec.getVCFGenotypeRecords();
List<String> knowngenotypes = new ArrayList<String>();
knowngenotypes.add("GG");
knowngenotypes.add("AG");
knowngenotypes.add("AA");
Assert.assertEquals(3, list.size());
for (Genotype g: list) {
for (VCFGenotypeRecord g: list) {
Assert.assertTrue(knowngenotypes.contains(g.getBases()));
}
}
@ -163,22 +162,11 @@ public class RodVCFTest extends BaseTest {
Iterator<RodVCF> iter = vcf.createIterator("VCF", vcfFile);
RodVCF rec = iter.next();
// we should get back a ref 'G' and an alt 'A'
List<Genotype> list = rec.getGenotypes();
List<VCFGenotypeRecord> list = rec.getVCFGenotypeRecords();
Assert.assertEquals(4.8,list.get(0).getNegLog10PError(),0.0001);
Assert.assertEquals(4.8,list.get(1).getNegLog10PError(),0.0001);
Assert.assertEquals(4.3,list.get(2).getNegLog10PError(),0.0001);
}
@Test
public void testHasGenotypes() {
RodVCF vcf = getVCFObject();
Iterator<RodVCF> iter = vcf.createIterator("VCF", vcfFile);
RodVCF rec = iter.next();
Assert.assertTrue(rec.hasGenotype(DiploidGenotype.valueOf("GG")));
Assert.assertTrue(rec.hasGenotype(DiploidGenotype.valueOf("AG")));
Assert.assertTrue(rec.hasGenotype(DiploidGenotype.valueOf("AA")));
Assert.assertTrue(!rec.hasGenotype(DiploidGenotype.valueOf("TT")));
}
@Test
public void testGetLocation() {
RodVCF vcf = getVCFObject();

8
java/test/org/broadinstitute/sting/gatk/walkers/annotator/SecondBaseSkewIntegrationTest.java
View File

@ -14,7 +14,9 @@ import java.util.Arrays;
*/
public class SecondBaseSkewIntegrationTest extends WalkerTest {
@Test
// TODO -- reinstitute these integration tests when GELI -> VariantContext is supported
//@Test
public void secondBaseSkewTest() {
for ( int test = 1; test <= 2; test ++ ) {
String bamFilePath = VariantAnnotatorIntegrationTest.validationDataPath()+VariantAnnotatorIntegrationTest.secondBaseTestFile(test)+".a2b.recal.annotation_subset.bam";
@ -25,7 +27,7 @@ public class SecondBaseSkewIntegrationTest extends WalkerTest {
}
}
@Test
//@Test
public void testOnE2File() {
String test_args = "-T VariantAnnotator -A SecondBaseSkew "
+"-R " + seqLocation + "references/Homo_sapiens_assembly18/v0/Homo_sapiens_assembly18.fasta "
@ -41,7 +43,7 @@ public class SecondBaseSkewIntegrationTest extends WalkerTest {
}
@Test
//@Test
public void testOnUnannotatedFile() {
String test_args = "-T VariantAnnotator -A SecondBaseSkew "
+"-R " + seqLocation + "references/Homo_sapiens_assembly18/v0/Homo_sapiens_assembly18.fasta "

48
java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java
View File

@ -99,6 +99,30 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
executeTest("testConfidence", spec);
}
// --------------------------------------------------------------------------------------------------------------
//
// testing beagle output
//
// --------------------------------------------------------------------------------------------------------------
@Test
public void testOtherOutput() {
String[] md5s = {"8c7dd53a402b727753002ebcd76168ac", "8cba0b8752f18fc620b4697840bc7291"};
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T UnifiedGenotyper" +
" -R " + oneKGLocation + "reference/human_b36_both.fasta" +
" -I " + validationDataLocation + "low_coverage_CEU.chr1.10k-11k.bam" +
" -varout %s" +
" -beagle %s" +
" -L 1:10,023,400-10,024,000" +
" -bm empirical" +
" -gm JOINT_ESTIMATE" +
" -vf VCF",
2,
Arrays.asList(md5s));
executeTest(String.format("testOtherOutput"), spec);
}
// --------------------------------------------------------------------------------------------------------------
//
// testing other output formats
@ -183,28 +207,4 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
executeTest(String.format("testMultiTechnologies"), spec);
}
// --------------------------------------------------------------------------------------------------------------
//
// testing beagle output
//
// --------------------------------------------------------------------------------------------------------------
@Test
public void testOtherOutput() {
String[] md5s = {"78482125d51f9eb2ee850a6b25921e84", "8cba0b8752f18fc620b4697840bc7291"};
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T UnifiedGenotyper" +
" -R " + oneKGLocation + "reference/human_b36_both.fasta" +
" -I " + validationDataLocation + "low_coverage_CEU.chr1.10k-11k.bam" +
" -varout %s" +
" -beagle %s" +
" -L 1:10,023,400-10,024,000" +
" -bm empirical" +
" -gm JOINT_ESTIMATE" +
" -vf GELI",
2,
Arrays.asList(md5s));
executeTest(String.format("testOtherOutput"), spec);
}
}

14
java/test/org/broadinstitute/sting/gatk/walkers/varianteval/VariantEvalWalkerIntegrationTest.java
View File

@ -116,7 +116,7 @@ public class VariantEvalWalkerIntegrationTest extends WalkerTest {
@Test
public void testEvalGenotypeROD() {
List<String> md5 = new ArrayList<String>();
md5.add("cbe45debc69e83d2b90988ee72042074");
md5.add("6ed44fd586c89dafd40cb8e0194dc456");
/**
* the above MD5 was calculated after running the following command:
*
@ -150,7 +150,7 @@ public class VariantEvalWalkerIntegrationTest extends WalkerTest {
@Test
public void testEvalMarksGenotypingExample() {
List<String> md5 = new ArrayList<String>();
md5.add("2694cf3eb73a9fecdda8cf5b76d0135d");
md5.add("c0396cfe89a63948aebbbae0a0e06678");
/**
* Run with the following commands:
*
@ -193,11 +193,11 @@ public class VariantEvalWalkerIntegrationTest extends WalkerTest {
@Test
public void testVCFVariantEvals() {
HashMap<String, String> md5 = new HashMap<String, String>();
md5.put("", "a1c5c8e9b185108969528e8c9fbef15e");
md5.put("-A", "2ae1456de2375502689c4af8f40b8693");
md5.put("-A --includeFilteredRecords", "9b2f446f0f42ab10d9d27f5221748f5e");
md5.put("-A --sampleName NA12878", "38deda8ab3816f083f0aeb97ccf8c347");
md5.put("-A -vcfInfoSelector AF=0.50", "715ea2da250f58aa35004386a890f5ba");
md5.put("", "ee6b096169d6c5e2ce49d394fbec799b");
md5.put("-A", "a443193c0810363f85278b1cfaed2fff");
md5.put("-A --includeFilteredRecords", "812d7f2ecac28b1be7e7028af17df9c0");
md5.put("-A --sampleName NA12878", "a443193c0810363f85278b1cfaed2fff");
md5.put("-A -vcfInfoSelector AF=0.50", "afed4bf0c9f11b86f6e5356012f9cf2d");
for ( Map.Entry<String, String> e : md5.entrySet() ) {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(

4
java/test/org/broadinstitute/sting/playground/gatk/walkers/variantstovcf/VariantsToVCFIntegrationTest.java
View File

@ -18,7 +18,7 @@ import java.io.File;
public class VariantsToVCFIntegrationTest extends WalkerTest {
@Test
//@Test
public void testVariantsToVCFUsingGeliInput() {
List<String> md5 = new ArrayList<String>();
md5.add("a94c15f2e8905fd3e98301375cf0f42a");
@ -47,7 +47,7 @@ public class VariantsToVCFIntegrationTest extends WalkerTest {
List<File> result = executeTest("testVariantsToVCFUsingGeliInput", spec).getFirst();
}
@Test
//@Test
public void testGenotypesToVCFUsingGeliInput() {
List<String> md5 = new ArrayList<String>();
md5.add("6b18f33e25edbd2154c17a949656644b");

128
java/test/org/broadinstitute/sting/utils/genotype/glf/GLFWriterTest.java
View File

@ -5,8 +5,7 @@ import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.GenomeLocParser;
import org.broadinstitute.sting.utils.StingException;
import org.broadinstitute.sting.utils.fasta.IndexedFastaSequenceFile;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.genotype.GenotypeWriter;
import org.broadinstitute.sting.utils.genotype.LikelihoodObject;
import org.junit.Assert;
import org.junit.Before;
import org.junit.BeforeClass;
@ -14,8 +13,8 @@ import org.junit.Test;
import java.io.File;
import java.io.FileNotFoundException;
import java.util.ArrayList;
import java.util.List;
import net.sf.samtools.SAMSequenceRecord;
/*
@ -55,9 +54,8 @@ public class GLFWriterTest extends BaseTest {
/** some made up values that we use to generate the GLF */
private final String header = "";
private static final int GENOTYPE_COUNT = 10;
private GenotypeWriter rec;
private GLFWriter rec;
protected static final String[] genotypes = {"AA", "AC", "AG", "AT", "CC", "CG", "CT", "GG", "GT", "TT"};
private static IndexedFastaSequenceFile seq;
protected final static double SIGNIFICANCE = 5.1;
@Before
@ -67,6 +65,7 @@ public class GLFWriterTest extends BaseTest {
@BeforeClass
public static void beforeTests() {
IndexedFastaSequenceFile seq;
try {
seq = new IndexedFastaSequenceFile(new File(oneKGLocation + "reference/human_b36_both.fasta"));
} catch (FileNotFoundException e) {
@ -77,37 +76,31 @@ public class GLFWriterTest extends BaseTest {
}
/**
* create a fake genotype
* create a fake genotype likehoods set
* @param bestGenotype the best genotype, as an index into the array of values
* @param location the location we're creating the genotype at
* @param ref the reference base
* @return a FakeGenotype (a fake genotype)
* @return fake genotype likelihoods
*/
private FakeGenotype createGenotype(int bestGenotype, GenomeLoc location, char ref) {
private LikelihoodObject createLikelihoods(int bestGenotype) {
double lk[] = new double[GENOTYPE_COUNT];
for (int x = 0; x < GENOTYPE_COUNT; x++) {
lk[x] = -15.0 - (double) x; // they'll all be unique like a snowflake
}
lk[bestGenotype] = -10.0; // lets make the best way better
return new FakeGenotype(location, genotypes[bestGenotype], ref, SIGNIFICANCE, lk);
return new LikelihoodObject(lk, LikelihoodObject.LIKELIHOOD_TYPE.NEGATIVE_LOG);
}
/**
* create a fake genotype with a minimum likelihood greater than 255
* create a fake genotype likelihhods set with a minimum likelihood greater than 255
* @param bestGenotype the best genotype, as an index into the array of values
* @param location the location we're creating the genotype at
* @param ref the reference base
* @return a FakeGenotype (a fake genotype)
* @return fake genotype likelihoods
*/
private FakeGenotype createGreaterThan255MinimumGenotype(int bestGenotype, GenomeLoc location, char ref) {
private LikelihoodObject createGreaterThan255MinimumGenotype(int bestGenotype) {
double lk[] = new double[GENOTYPE_COUNT];
for (int x = 0; x < GENOTYPE_COUNT; x++) {
lk[x] = -355.0 - (double) x; // they'll all be unique like a snowflake
}
lk[bestGenotype] = -256.0; // lets make the best way better
return new FakeGenotype(location, genotypes[bestGenotype], ref, SIGNIFICANCE, lk);
return new LikelihoodObject(lk, LikelihoodObject.LIKELIHOOD_TYPE.NEGATIVE_LOG);
}
@ -120,11 +113,10 @@ public class GLFWriterTest extends BaseTest {
writeTo.deleteOnExit();
rec = new GLFWriter(writeTo);
((GLFWriter)rec).writeHeader(header);
rec.writeHeader(header);
for (int x = 0; x < 100; x++) {
GenomeLoc loc = GenomeLocParser.createGenomeLoc(1, x + 1);
Genotype type = createGenotype(x % 10, loc, 'A');
rec.addGenotypeCall(type);
rec.addCall(new SAMSequenceRecord("test", 0), (int)loc.getStart(), 10, 'A', 9, createLikelihoods(x % 10));
}
rec.close();
@ -139,11 +131,10 @@ public class GLFWriterTest extends BaseTest {
writeTo.deleteOnExit();
rec = new GLFWriter(writeTo);
((GLFWriter)rec).writeHeader(header);
rec.writeHeader(header);
for (int x = 0; x < 5; x++) {
GenomeLoc loc = GenomeLocParser.createGenomeLoc(1, x + 1);
Genotype type = createGreaterThan255MinimumGenotype(x % 10, loc, 'A');
rec.addGenotypeCall(type);
rec.addCall(new SAMSequenceRecord("test", 0), (int)loc.getStart(), 10, 'A', 9, createGreaterThan255MinimumGenotype(x % 10));
}
rec.close();
@ -158,96 +149,19 @@ public class GLFWriterTest extends BaseTest {
public void basicWriteThenRead() {
File writeTo = new File("testGLF2.glf");
writeTo.deleteOnExit();
List<FakeGenotype> types = new ArrayList<FakeGenotype>();
rec = new GLFWriter(writeTo);
((GLFWriter)rec).writeHeader(header);
rec.writeHeader(header);
for (int x = 0; x < 100; x++) {
GenomeLoc loc = GenomeLocParser.createGenomeLoc(1, x + 1);
FakeGenotype type = createGenotype(x % 10, loc, 'A');
types.add(type);
rec.addGenotypeCall(type);
rec.addCall(new SAMSequenceRecord("test", 0), (int)loc.getStart(), 10, 'A', 9, createLikelihoods(x % 10));
}
rec.close();
GLFReader reader = new GLFReader(writeTo);
int count = 0;
while (reader.hasNext()) {
GLFRecord rec = reader.next();
Assert.assertTrue(types.get(count).compareTo(FakeGenotype.toFakeGenotype((GLFSingleCall) rec, rec.getContig(), (int)rec.getPosition())) == 0);
reader.next();
count++;
}
Assert.assertEquals(count, 100);
}
}
class FakeGenotype extends GLFGenotypeCall implements Comparable<FakeGenotype> {
private double[] likelihoods;
/**
* create a basic genotype, given the following fields
*
* @param location the genomic location
* @param genotype the genotype, as a string, where ploidy = string.length
* @param ref the reference base as a char
* @param negLog10PError the confidence score
*/
public FakeGenotype(GenomeLoc location, String genotype, char ref, double negLog10PError, double likelihoods[]) {
super(Character.toString(ref), location);
setLikelihoods(likelihoods);
setGenotype(genotype);
setNegLog10PError(negLog10PError);
}
/**
* get the likelihood information for this
*
* @return
*/
@Override
public double[] getLikelihoods() {
return likelihoods;
}
public void setLikelihoods(double[] likelihoods) {
this.likelihoods = likelihoods;
}
public int compareTo(FakeGenotype that) {
if (this.getLocation().compareTo(that.getLocation()) != 0) {
System.err.println("Location's aren't equal; this = " + this.getLocation() + " that = " + that.getLocation());
return this.getLocation().compareTo(that.getLocation());
}
if (!this.getBases().equals(that.getBases())) {
System.err.println("getBases's aren't equal; this = " + this.getBases() + " that = " + that.getBases());
return -1;
}
for (int x = 0; x < this.likelihoods.length; x++) {
if (this.likelihoods[x] != that.getLikelihoods()[x]) {
System.err.println("likelihoods' aren't equal; this = " + this.likelihoods[x] + " that = " + that.getLikelihoods()[x]);
return -1;
}
}
return 0;
}
public static FakeGenotype toFakeGenotype(GLFSingleCall record, String contig, int postition) {
double likelihoods[] = record.getLikelihoods();
char ref = record.getRefBase().toChar();
double significance = GLFWriterTest.SIGNIFICANCE;
int minIndex = 0;
for (int i = 0; i < likelihoods.length; i++) {
if (likelihoods[i] < likelihoods[minIndex]) minIndex = i;
}
for (int i = 0; i < likelihoods.length; i++) {
likelihoods[i] = likelihoods[i] * -1;
}
String genotype = GLFWriterTest.genotypes[minIndex];
GenomeLoc loc = GenomeLocParser.createGenomeLoc(contig, postition);
return new FakeGenotype(loc, genotype, ref, significance, likelihoods);
}
}

9
java/test/org/broadinstitute/sting/utils/genotype/vcf/VCFReaderTest.java
View File

@ -15,6 +15,7 @@ import org.broadinstitute.sting.BaseTest;
import org.broadinstitute.sting.gatk.refdata.RodVCF;
import org.broadinstitute.sting.utils.GenomeLocParser;
import org.broadinstitute.sting.utils.StingException;
import org.broadinstitute.sting.utils.genotype.Variation;
import org.broadinstitute.sting.utils.fasta.IndexedFastaSequenceFile;
import org.junit.Assert;
import org.junit.BeforeClass;
@ -291,7 +292,7 @@ public class VCFReaderTest extends BaseTest {
for ( VCFGenotypeRecord grec : grecords ) {
if ( !grec.isEmptyGenotype() ) {
Assert.assertTrue(grec.isVariant(rec.getReference().charAt(0)));
Assert.assertEquals(rec, grec.toVariation(rec.getReference().charAt(0)));
Assert.assertEquals(rec, grec.getRecord());
}
}
@ -318,19 +319,19 @@ public class VCFReaderTest extends BaseTest {
rec = reader.next();
Assert.assertTrue(!rec.isFiltered());
Assert.assertTrue(rec.getFilterString().equals("."));
Assert.assertEquals(VCFRecord.VARIANT_TYPE.SNP, rec.getType());
Assert.assertEquals(Variation.VARIANT_TYPE.SNP, rec.getType());
// record #9: deletion
if (!reader.hasNext()) Assert.fail("The reader should have a record");
rec = reader.next();
Assert.assertEquals(VCFRecord.VARIANT_TYPE.DELETION, rec.getType());
Assert.assertEquals(Variation.VARIANT_TYPE.DELETION, rec.getType());
Assert.assertEquals(1, rec.getAlternateAlleleList().size());
Assert.assertTrue(rec.getAlternateAlleleList().get(0).equals(""));
// record #10: insertion
if (!reader.hasNext()) Assert.fail("The reader should have a record");
rec = reader.next();
Assert.assertEquals(VCFRecord.VARIANT_TYPE.INSERTION, rec.getType());
Assert.assertEquals(Variation.VARIANT_TYPE.INSERTION, rec.getType());
Assert.assertEquals(rec.getAlternateAlleleList().size(), 1);
Assert.assertTrue(rec.getAlternateAlleleList().get(0).equals("CAT"));