Bug Bulletin: The recent 3.2 release fixes many issues. If you run into a problem, please try the latest version before posting a bug report, as your problem may already have been solved.

Using the GATK API to annotate my VCF (2)

lindenblindenb Posts: 5Member

So I started to play with the GATK API to annotate my VCF. :-)

My first simple aim is to translate mt VCFBigWig to the GATK . Here is my code so far (see below)

I compiled it with:

javac -cp /home/lindenb/package/GenomeAnalysisTK-2.4-3-g2a7af43/GenomeAnalysisTK.jar:/home/lindenb/package/bigwig-read-only/dist/BigWig.jar:/home/lindenb/package/bigwig-read-only/lib/log4j-1.2.15.jar -sourcepath src/main/java \
./src/main/java/com/github/lindenb/gatk/tools/vcfbigwig/VCFBigWig.java

(it is compiled but there is a bunch of warnings like "/home/lindenb/package/GenomeAnalysisTK-2.4-3-g2a7af43/GenomeAnalysisTK.jar(org/broadinstitute/sting/gatk/GenomeAnalysisEngine.class): warning: Cannot find annotation method 'value()' in type 'Ensures': class file for com.google.java.contract.Ensures not found"

And I tried to run my code using various invocations of java, like

java  -cp /home/lindenb/package/GenomeAnalysisTK-2.4-3-g2a7af43/GenomeAnalysisTK.jar:/home/lindenb/package/bigwig-read-only/dist/BigWig.jar:/home/lindenb/package/bigwig-read-only/lib/log4j-1.2.15.jar:. org.broadinstitute.sting.gatk.CommandLineGATK \
-T VariantAnnotator \
-A com.github.lindenb.gatk.tools.vcfbigwig.VCFBigWig \
-R /test.fa test.vcf.gz

but it raises an exeception

java.lang.ExceptionInInitializerError
at org.broadinstitute.sting.gatk.walkers.annotator.VariantAnnotatorEngine.initializeAnnotations(VariantAnnotatorEngine.java:116)
(...)
Caused by: java.lang.NullPointerException
at org.broadinstitute.sting.utils.classloader.JVMUtils.isAnonymous(JVMUtils.java:91)
at org.broadinstitute.sting.utils.classloader.PluginManager.<init>(PluginManager.java:158)
at org.broadinstitute.sting.utils.classloader.PluginManager.<init>(PluginManager.java:107)
at org.broadinstitute.sting.gatk.walkers.annotator.interfaces.AnnotationInterfaceManager.<clinit>(AnnotationInterfaceManager.java:34)
... 9 more

Questions:

  • What the proper to declare my class VCFBigWig , how can I make it available for the VariantAnnotatorEngine ?
  • I used the @Argument annotation, is it the right way to catch an argument from the command line?
  • Can I use the same Annotator for one GATK invocation but with different arguments ?
  • there is a "initialize" method; Is there a "dispose" method to release the resources ?
  • what is the interface "RodRequiringAnnotation" ?
  • is "UserException" the right way to send an error ?
  • Should an InfoFieldAnnotation or a Walker be thread free (like a HttpServlet) ? I mean can I safely use class members ?
  • I saw in VariantAnnotatorEngine that an annotator can either change the ID or the INFO. How should I design my class in order to alter both fields ? For example my tool VCFTabix with use another VCF file (e.g: the VCF for 1KGenomes) to update the ID and the INFO fields.

Thank you for your help,

Pierre

.

package com.github.lindenb.gatk.tools.vcfbigwig;

import java.io.IOException;
import java.util.Arrays;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Set;

import org.broad.igv.bbfile.BBFileReader;
import org.broad.igv.bbfile.BigWigIterator;
import org.broad.igv.bbfile.WigItem;
import org.broadinstitute.sting.commandline.Argument;
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.AnnotatorCompatible;
import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.InfoFieldAnnotation;
import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.RodRequiringAnnotation;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap;
import org.broadinstitute.variant.variantcontext.VariantContext;
import org.broadinstitute.variant.vcf.VCFHeaderLine;
import org.broadinstitute.variant.vcf.VCFHeaderLineType;
import org.broadinstitute.variant.vcf.VCFInfoHeaderLine;


public class VCFBigWig
extends InfoFieldAnnotation
implements RodRequiringAnnotation
{
/** PATH to the bigwig file */
@Argument(fullName="bigWigFile",shortName="bw",doc="BigWig to process",required=true)
private String bigWigFile=null;
/** the ID in the INFO field */
@Argument(fullName = "bigwigid", shortName = "bwid", doc="bigwig ID ", required=true)
private String bigWigId=null;
/** BBFileReader to the bigwig file */
private BBFileReader bbFileReader=null;

public VCFBigWig()
{
}

@Override
public void initialize(AnnotatorCompatible walker,
GenomeAnalysisEngine toolkit, Set<VCFHeaderLine> headerLines) {
super.initialize(walker, toolkit, headerLines);
/* open the BIGFile, check it is a bigwig */
try {
this.bbFileReader=new BBFileReader(bigWigFile);
}
catch (IOException e) {
throw new UserException("Cannot read "+bigWigFile,e);
}
if(!this.bbFileReader.isBigWigFile())
{
throw new UserException("The following file is not a BigWig File : "+bigWigFile);
}
}


@Override
public Map<String, Object> annotate(
RefMetaDataTracker tracker,
AnnotatorCompatible walker,
ReferenceContext ref,
Map<String, AlignmentContext> stratifiedContexts,
VariantContext vc,
Map<String, PerReadAlleleLikelihoodMap> stratifiedPerReadAlleleLikelihoodMap)
{
return annotate(tracker, walker,ref,stratifiedContexts,vc);
}



@Override
public Map<String, Object> annotate(
RefMetaDataTracker tracker,
AnnotatorCompatible walker,
ReferenceContext ref,
Map<String, AlignmentContext> stratifiedContexts,
VariantContext vc)
{
final boolean contained=true;
double sum=0;
int count=0;
BigWigIterator iter=this.bbFileReader.getBigWigIterator(
vc.getChr(), vc.getStart()-1,
vc.getChr(), vc.getEnd(),
contained
);

while(iter.hasNext())
{
WigItem item=iter.next();
float v=item.getWigValue();
sum+=v;
count++;
}

if(count==0)
{
return Collections.emptyMap();
}
Map<String,Object> m=new HashMap<String, Object>(1);
m.put(this.bigWigId, (float)(sum/count));
return m ;
}
@Override
public List<VCFInfoHeaderLine> getDescriptions()
{
return Arrays.asList(new VCFInfoHeaderLine[]{
new VCFInfoHeaderLine(bigWigId,1, VCFHeaderLineType.Float,"Mean values for bigwig file "+this.bigWigFile)
});
}
@Override
public List<String> getKeyNames() {
return Arrays.asList(new String[]{
bigWigId
});
}

}

Best Answer

Sign In or Register to comment.