Escribiendo un ReadFilter personalizado para el GATK, mi cuaderno.

El GATK contiene un conjunto de filtros de lectura predefinidos que «filtrar o transferir archivos de datos SAM/BAM entrantes»:

Con la ayuda de la arquitectura modular del GATK, es posible escribir un Filtro de lectura personalizado. En este post voy a escribir un LeerFiltro que elimina las lecturas si contienen una secuencia específica en un segmento de recorte suave.
Los filtros de lectura amplían la clase. org.broadinstitute.gatk.engine.filters.ReadFilter:

public abstract class ReadFilter implements SamRecordFilter 

Así, la implementación de nuestro filtro, llamado Mi se define a continuación:
La clase se compila y empaqueta utilizando el propio gatk como biblioteca:

mkdir -p  org/broadinstitute/gatk/engine/filters
cp MyFilter.java org/broadinstitute/gatk/engine/filters
javac -cp /path/to/GenomeAnalysisTK.jar org/broadinstitute/gatk/engine/filters/MyFilter.java
jar cvf myfilter.jar org
rm -rf org

La clase principal de GATK es org.broadinstitute.gatk.engine.CommandLineGATK, ahora podemos verificar que hay un filtro llamado Mi en la lista de filtros de lectura.

$ java -cp myfilter.jar:/path/to/GenomeAnalysisTK.jar  org.broadinstitute.gatk.engine.CommandLineGATK  
  -T ReadLengthDistribution 
  --read_filter generate_and_error_please
(...)
##### ERROR MESSAGE: Invalid command line: Malformed read filter: Read filter generate_and_error_please not found. Available read filters:
##### ERROR 
##### ERROR                                 FilterName        Documentation
##### ERROR                           MissingReadGroup        https://www.broadinstitute.org/gatk/guide/tooldocs/org_broadinstitute_gatk_engine_filters_MissingReadGroupFilter.php
##### ERROR                                         My        https://www.broadinstitute.org/gatk/guide/tooldocs/org_broadinstitute_gatk_engine_filters_MyFilter.php
##### ERROR                    NoOriginalQualityScores        https://www.broadinstitute.org/gatk/guide/tooldocs/org_broadinstitute_gatk_engine_filters_NoOriginalQualityScoresFilter.php

Ahora podemos usar las herramientas GATK con o sin el filtro ‘Mi’.

java -cp myfilter.jar:/path/to/GenomeAnalysisTK.jar org.broadinstitute.gatk.engine.CommandLineGATK 
  -T ReadLengthDistribution 
  -R /path/to/ref.fa 
  -I /path/to/test.bam 
  -L 22 -o dist.ATGATA.tbl 
  --read_filter My --clippattern ATGATA 
(...)
INFO  10:53:32,412 MicroScheduler - 32 reads were filtered out during the traversal out of approximately 12434 total reads (0.26%) 
INFO  10:53:32,412 MicroScheduler -   -> 0 reads (0.00% of total) failing BadCigarFilter 
INFO  10:53:32,412 MicroScheduler -   -> 0 reads (0.00% of total) failing MalformedReadFilter 
INFO  10:53:32,413 MicroScheduler -   -> 32 reads (0.26% of total) failing MyFilter

$ cat dist.ATGATA.tbl
#:GATKReport.v1.1:1
#:GATKTable:2:130:%s:%s:;
#:GATKTable:ReadLengthDistribution:Table of read length distributions
readLength  Sample
        19        6
        20        4
        21        4
        22       13
        23        4
        24       12


mientras que sin ‘Mi’ la salida es:

java -cp myfilter.jar:/path/to/GenomeAnalysisTK.jar org.broadinstitute.gatk.engine.CommandLineGATK 
  -T ReadLengthDistribution 
  -R /path/to/ref.fa 
  -I /path/to/test.bam 
  -L 22 -o dist.all.tbl
(...)
INFO  10:53:35,713 MicroScheduler - 0 reads were filtered out during the traversal out of approximately 12434 total reads (0.00%) 
INFO  10:53:35,713 MicroScheduler -   -> 0 reads (0.00% of total) failing BadCigarFilter 
INFO  10:53:35,714 MicroScheduler -   -> 0 reads (0.00% of total) failing MalformedReadFilter 
------------------------------------------------------------------------------------------
Done. There were no warn messages.
------------------------------------------------------------------------------------------

$ cat  dist.all.tbl
#:GATKReport.v1.1:1
#:GATKTable:2:130:%s:%s:;
#:GATKTable:ReadLengthDistribution:Table of read length distributions
readLength  Sample
        19        6
        20        4
        21        4
        22       13
        23        4
        24       12

El archivo MAKE

Eso es todo,
Pedro

Fuente del artículo

Deja un comentario