YOKOFAKUN: Creando un GATK Walker personalizado (GATK 3.6) : mi cuaderno

Este es mi cuaderno para crear un motor personalizado en GATK.

Quiero leer un archivo VCF y obtener una tabla de categoría/recuento. Algo como esto:

creo una clase Categoría describiendo cada fila de la tabla. Es solo una lista de cadenas

static class Category
        implements Comparable<Category>
        {
        private final List<String> labels ;
        Category(final List<String> labels) 

como iban a usar Categoría en una matriz asociativa / Mapa que necesitamos implementar hashCode y equals:

static class Category
        implements Comparable<Category>
        

# El Caminante Principal:

El motor principal se llama CountPredictions extiende un RodWalker. Por qué ? Realmente no sé, he copiado esto de otro motor. Vamos a hacer un mapa/reducir un ‘Mapa’ para que mi clase se declare como

public class CountPredictions 
  extends RodWalker< Map<CountPredictions.Category,Long>, Map<CountPredictions.Category,Long>>
  

# Documentando al Caminante

Se agrega una anotación ‘DocumentedGATKFeature’ para que GATKEngine pueda encontrar nuestro caminante

  (...)
  @DocumentedGATKFeature(
        summary="Count Predictions",
        groupName = HelpConstants.DOCS_CAT_VARMANIP,
        extraDocs =  )
public class CountPredictions  extends RodWalker
  (....)

La entrada es un archivo VCF:

@Input(fullName="variant", shortName = "V", doc="Input VCF file", required=true)
    public RodBinding<VariantContext> variants;

La salida es un PrintStream (el valor predeterminado es stdout) donde escribiremos la tabla:

@Output(doc="File to which result should be written")
    protected PrintStream out = System.out;

Los otros argumentos también se decoran con anotaciones de Java. Estos son los interruptores para agregar algunas columnas a la tabla final:

(...)
    @Argument(fullName="chrom",shortName="chrom",required=false,doc="Group by Chromosome/Contig")
    public boolean bychrom = false;
    @Argument(fullName="ID",shortName="ID",required=false,doc="Group by having/not-having ID")
    public boolean byID = false;
    @Argument(fullName="variantType",shortName="variantType",required=false,doc="Group by VariantType")
    (...)

este método se llama después de analizar los argumentos. Como queremos poder analizar la anotación VEP donde vamos a obtener el encabezado VCF y extraer los componentes del ANA atributo para obtener los índices de ‘Annotation_Impact’ y ‘Transcript_BioType’

@Override
    public void initialize()  bybiotype) 

        super.initialize();
    

según tengo entendido, ‘reduceInit()’ se usa para crear un primer elemento durante el proceso de mapa/reducción. Así que estamos creando un mapa asociativo vacío:

@Override
    public Map<CountPredictions.Category,Long> reduceInit() 

Este método reduce dos procesos de mapeo, por lo que estamos creando un mapa que combina ambos recuentos:

@Override
    public Map<CountPredictions.Category,Long> reduce(Map<CountPredictions.Category,Long> value, Map<CountPredictions.Category,Long> sum) 

Este es el caballo de batalla del motor. Por lo que puedo ver, se llama para cada variante. El método ‘tracker.getValues()’ devuelve una matriz de todas las variantes (aquí, solo una porque solo tenemos una entrada) en la posición actual. Al final del método, debemos devolver un recuento de Categoría para esas variantes…

@Override
    public Map<CountPredictions.Category,Long> map(final RefMetaDataTracker tracker,final ReferenceContext ref, final AlignmentContext context) 

para cada variante, se llena una nueva Categoría:

(...)
    final List<String> labels=new ArrayList<>();
    if(bychrom) labels.add(ctx.getContig());
    if(byID) labels.add(ctx.hasID()?"Y":".");
    (...)
    final Category cat=new Category(labels);
    Long n=count.get(cat);
    count.put(cat, n==null?1L:n+1);
    (...)

Cuando se llama a onTraversalDone, la tabla se imprime utilizando la clase GATKReport:

@Override
    public void onTraversalDone(final Map<CountPredictions.Category,Long> counts) 

el código fuente completo:

package mygatk;

import java.io.PrintStream;
import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.regex.Pattern;

import org.broadinstitute.gatk.engine.CommandLineGATK;
import org.broadinstitute.gatk.engine.GATKVCFUtils;
import org.broadinstitute.gatk.engine.walkers.RodWalker;
import org.broadinstitute.gatk.utils.commandline.Argument;
import org.broadinstitute.gatk.utils.commandline.Input;
import org.broadinstitute.gatk.utils.commandline.Output;
import org.broadinstitute.gatk.utils.commandline.RodBinding;
import org.broadinstitute.gatk.utils.contexts.AlignmentContext;
import org.broadinstitute.gatk.utils.contexts.ReferenceContext;
import org.broadinstitute.gatk.utils.help.DocumentedGATKFeature;
import org.broadinstitute.gatk.utils.help.HelpConstants;
import org.broadinstitute.gatk.utils.refdata.RefMetaDataTracker;
import org.broadinstitute.gatk.utils.report.GATKReport;
import org.broadinstitute.gatk.utils.report.GATKReportTable;

import htsjdk.variant.variantcontext.Genotype;
import htsjdk.variant.variantcontext.VariantContext;
import htsjdk.variant.vcf.VCFHeader;
import htsjdk.variant.vcf.VCFInfoHeaderLine;

/**
 * Test Documentation
 */
@DocumentedGATKFeature(
        summary="Count Predictions",
        groupName = HelpConstants.DOCS_CAT_VARMANIP,
        extraDocs =  )
public class CountPredictions  extends RodWalker<Map<CountPredictions.Category,Long>, Map<CountPredictions.Category,Long>> ]");

    static class Category
        implements Comparable<Category>
        
    private int ann_impact_column=-1;/* eg: MODIFIER / LOW */
    private int ann_transcript_biotype_column=-1;/* eg: transcript /intergenic_region */

    @Override
    public void initialize() 

    @Override
    public Map<CountPredictions.Category,Long> map(final RefMetaDataTracker tracker,final ReferenceContext ref, final AlignmentContext context) 

    @Override
    public Map<CountPredictions.Category,Long> reduce(Map<CountPredictions.Category,Long> value, Map<CountPredictions.Category,Long> sum) 

    @Override
    public Map<CountPredictions.Category,Long> reduceInit() 


    @Override
    public void onTraversalDone(final Map<CountPredictions.Category,Long> counts) 

    

gatk.jar=/path/to/GenomeAnalysisTK.jar
VCF=input.vcf.gz
test: dist/mygatk.jar
        java -cp $:dist/mygatk.jar 
                org.broadinstitute.gatk.engine.CommandLineGATK 
                -R /path/to/ref.fasta 
                -T CountPredictions -V $ -ID

dist/mygatk.jar : $ ./mygatk/CountPredictions.java
        mkdir -p tmp dist
        javac -d tmp -cp .:$  ./mygatk/CountPredictions.java
        jar cfv [email protected] -C tmp .
        rm -rf tmp

INFO  16:47:16,545 HelpFormatter - ---------------------------------------------------------------------------------- 
INFO  16:47:16,547 HelpFormatter - The Genome Analysis Toolkit (GATK) v3.6-0-g89b7209, Compiled 2016/06/01 22:27:29 
INFO  16:47:16,547 HelpFormatter - Copyright (c) 2010-2016 The Broad Institute 
INFO  16:47:16,547 HelpFormatter - For support and documentation go to https://www.broadinstitute.org/gatk 
INFO  16:47:16,714 HelpFormatter - [Fri Jan 13 16:47:16 CET 2017] Executing on Linux 3.10.0-327.36.3.el7.x86_64 amd64 
INFO  16:47:16,714 HelpFormatter - Java HotSpot(TM) 64-Bit Server VM 1.8.0_102-b14 JdkDeflater 
INFO  16:47:16,717 HelpFormatter - Program Args: -R ref.fasta -T CountPredictions -V input.vcf.gz -ID 
INFO  16:47:16,724 HelpFormatter - Date/Time: 2017/01/13 16:47:16 
INFO  16:47:16,725 HelpFormatter - ---------------------------------------------------------------------------------- 
INFO  16:47:16,725 HelpFormatter - ---------------------------------------------------------------------------------- 
INFO  16:47:16,740 GenomeAnalysisEngine - Strictness is SILENT 
INFO  16:47:16,829 GenomeAnalysisEngine - Downsampling Settings: Method: BY_SAMPLE, Target Coverage: 1000 
WARN  16:47:16,881 IndexDictionaryUtils - Track variant doesn't have a sequence dictionary built in, skipping dictionary validation 
INFO  16:47:16,936 GenomeAnalysisEngine - Preparing for traversal 
INFO  16:47:16,940 GenomeAnalysisEngine - Done preparing for traversal 
INFO  16:47:16,940 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING] 
INFO  16:47:16,941 ProgressMeter -                 | processed |    time |    per 1M |           |   total | remaining 
INFO  16:47:16,941 ProgressMeter -        Location |     sites | elapsed |     sites | completed | runtime |   runtime 
INFO  16:47:46,946 ProgressMeter -  chr22:26185095    167649.0    30.0 s       3.0 m       91.0%    32.0 s       2.0 s 
#:GATKReport.v1.1:1
#:GATKTable:2:1:%s:%s:;
#:GATKTable:Variants:Variants input.vcf.gz
IN_DBSNP  COUNT 
.         139984

INFO  16:47:52,501 CountPredictions - TraversalDone 
INFO  16:47:52,502 ProgressMeter -            done    202931.0    35.0 s       2.9 m       91.1%    38.0 s       3.0 s 
INFO  16:47:52,502 ProgressMeter - Total runtime 35.56 secs, 0.59 min, 0.01 hours 
------------------------------------------------------------------------------------------
Done. There were 1 WARN messages, the first 1 are repeated below.
WARN  16:47:16,881 IndexDictionaryUtils - Track variant doesn't have a sequence dictionary built in, skipping dictionary validation 
------------------------------------------------------------------------------------------