Saltar a contenido

Parte 2: Llamado de variantes por muestra

Traducción asistida por IA - más información y sugerencias

En la Parte 1, probaste los comandos de Samtools y GATK manualmente en sus respectivos contenedores. A continuación, vamos a envolver esos mismos comandos en un workflow de Nextflow.

Asignación

En esta parte del curso, vamos a desarrollar un workflow que hace lo siguiente:

  1. Generar un archivo de índice para cada archivo BAM de entrada usando Samtools
  2. Ejecutar GATK HaplotypeCaller en cada archivo BAM de entrada para generar llamados de variantes por muestra en formato VCF (Variant Call Format)
BAMSamtools indexBAM indexIntervalsReference+ index & dictGATK HaplotypeCallerVCF + index

Esto automatiza los pasos de la primera sección de Parte 1: Descripción del método, donde ejecutaste estos comandos manualmente en sus contenedores.

Como punto de partida, te proporcionamos un archivo de workflow, genomics.nf, que describe las partes principales del workflow, así como dos archivos de módulo, samtools_index.nf y gatk_haplotypecaller.nf, que describen la estructura de cada proceso.

Archivos esqueleto
genomics.nf
#!/usr/bin/env nextflow

// Module INCLUDE statements

/*
 * Pipeline parameters
 */

// Primary input

workflow {

    main:
    // Create input channel

    // Call processes

    publish:
    // Declare outputs to publish
}

output {
    // Configure publish targets
}
modules/samtools_index.nf
#!/usr/bin/env nextflow

/*
 * Generate BAM index file
 */
process SAMTOOLS_INDEX {

    container

    input:

    output:

    script:
    """

    """
}
modules/gatk_haplotypecaller.nf
#!/usr/bin/env nextflow

/*
 * Call variants with GATK HaplotypeCaller
 */
process GATK_HAPLOTYPECALLER {

    container

    input:

    output:

    script:
    """

    """
}

Estos archivos no son funcionales; su propósito es solo servir como esqueletos para que los completes con las partes interesantes del código.

Plan de la lección

Para hacer el proceso de desarrollo más educativo, lo hemos dividido en cuatro etapas:

  1. Escribir un workflow de una sola etapa que ejecute Samtools index en un archivo BAM. Esto cubre la creación de un módulo, su importación y su llamado en un workflow.
  2. Agregar un segundo proceso para ejecutar GATK HaplotypeCaller en el archivo BAM indexado. Esto introduce el encadenamiento de salidas de procesos a entradas y el manejo de archivos accesorios.
  3. Adaptar el workflow para ejecutarse en un lote de muestras. Esto cubre la ejecución paralela e introduce tuplas para mantener archivos asociados juntos.
  4. Hacer que el workflow acepte un archivo de texto que contenga un lote de archivos de entrada. Esto demuestra un patrón común para proporcionar entradas en lote.

Cada etapa se enfoca en un aspecto específico del desarrollo de workflows.

Consejo

Asegúrate de estar en el directorio de trabajo correcto: cd /workspaces/training/nf4-science/genomics

1. Escribir un workflow de una sola etapa que ejecute Samtools index en un archivo BAM

Esta primera etapa se enfoca en lo básico: cargar un archivo BAM y generar un índice para él.

Recuerda el comando samtools index de la Parte 1:

samtools index '<input_bam>'

El comando toma un archivo BAM como entrada y produce un archivo de índice .bai junto a él. El URI del contenedor era community.wave.seqera.io/library/samtools:1.20--b5dfbd93de237464.

Vamos a tomar esta información y envolverla en Nextflow en tres etapas:

  1. Configurar la entrada
  2. Escribir el proceso de indexación y llamarlo en el workflow
  3. Configurar el manejo de la salida

1.1. Configurar la entrada

Necesitamos declarar un parámetro de entrada, crear un perfil de prueba para proporcionar un valor predeterminado conveniente y crear un canal de entrada.

1.1.1. Agregar una declaración de parámetro de entrada

En el archivo principal del workflow genomics.nf, bajo la sección Pipeline parameters, declara un parámetro CLI llamado input.

genomics.nf
 5
 6
 7
 8
 9
10
11
/*
 * Pipeline parameters
 */
params {
    // Primary input
    input: Path
}
genomics.nf
5
6
7
8
9
/*
 * Pipeline parameters
 */

// Primary input

Eso configura el parámetro CLI, pero no queremos escribir la ruta del archivo cada vez que ejecutemos el workflow durante el desarrollo. Hay múltiples opciones para proporcionar un valor predeterminado; aquí usamos un perfil de prueba.

1.1.2. Crear un perfil de prueba con un valor predeterminado en nextflow.config

Un perfil de prueba proporciona valores predeterminados convenientes para probar un workflow sin especificar entradas en la línea de comandos. Esta es una convención común en el ecosistema de Nextflow (consulta Hello Config para más detalles).

Agrega un bloque profiles a nextflow.config con un perfil test que establezca el parámetro input a uno de los archivos BAM de prueba.

nextflow.config
1
2
3
4
5
6
7
docker.enabled = true

profiles {
    test {
        params.input = "${projectDir}/data/bam/reads_mother.bam"
    }
}
nextflow.config
1
docker.enabled = true

Aquí, estamos usando ${projectDir}, una variable integrada de Nextflow que apunta al directorio donde se encuentra el script del workflow. Esto facilita hacer referencia a archivos de datos y otros recursos sin codificar rutas absolutas.

1.1.3. Configurar el canal de entrada

En el bloque workflow, crea un canal de entrada a partir del valor del parámetro usando el channel factory .fromPath (como se usó en Hello Channels).

genomics.nf
13
14
15
16
17
workflow {

    main:
    // Create input channel (single file via CLI parameter)
    reads_ch = channel.fromPath(params.input)
genomics.nf
13
14
15
16
workflow {

    main:
    // Create input channel

A continuación, necesitaremos crear el proceso para ejecutar la indexación en esta entrada.

1.2. Escribir el proceso de indexación y llamarlo en el workflow

Necesitamos escribir la definición del proceso en el archivo del módulo, importarlo al workflow usando una declaración include y llamarlo con la entrada.

1.2.1. Completar el módulo para el proceso de indexación

Abre modules/samtools_index.nf y examina el esquema de la definición del proceso. Deberías reconocer los elementos estructurales principales; si no, considera leer Hello Nextflow para refrescar la memoria.

Adelante, completa la definición del proceso por ti mismo usando la información proporcionada arriba, luego verifica tu trabajo contra la solución en la pestaña "Después" a continuación.

modules/samtools_index.nf
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
#!/usr/bin/env nextflow

/*
 * Generate BAM index file
 */
process SAMTOOLS_INDEX {

    container

    input:

    output:

    script:
    """

    """
}
modules/samtools_index.nf
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
#!/usr/bin/env nextflow

/*
 * Generate BAM index file
 */
process SAMTOOLS_INDEX {

    container 'community.wave.seqera.io/library/samtools:1.20--b5dfbd93de237464'

    input:
    path input_bam

    output:
    path "${input_bam}.bai"

    script:
    """
    samtools index '$input_bam'
    """
}

Una vez que hayas completado esto, el proceso está completo. Para usarlo en el workflow, necesitarás importar el módulo y agregar un llamado al proceso.

1.2.2. Incluir el módulo

En genomics.nf, agrega una declaración include para hacer que el proceso esté disponible para el workflow:

genomics.nf
3
4
// Module INCLUDE statements
include { SAMTOOLS_INDEX } from './modules/samtools_index.nf'
genomics.nf
3
// Module INCLUDE statements

El proceso ahora está disponible en el ámbito del workflow.

1.2.3. Llamar al proceso de indexación con la entrada

Ahora, agreguemos un llamado a SAMTOOLS_INDEX en el bloque workflow, pasando el canal de entrada como argumento.

genomics.nf
14
15
16
17
18
19
20
21
workflow {

    main:
    // Create input channel (single file via CLI parameter)
    reads_ch = channel.fromPath(params.input)

    // Create index file for input BAM file
    SAMTOOLS_INDEX(reads_ch)
genomics.nf
14
15
16
17
18
19
20
workflow {

    main:
    // Create input channel (single file via CLI parameter)
    reads_ch = channel.fromPath(params.input)

    // Call processes

El workflow ahora carga la entrada y ejecuta el proceso de indexación en ella. A continuación, necesitamos configurar cómo se publica la salida.

1.3. Configurar el manejo de la salida

Necesitamos declarar qué salidas de procesos publicar y especificar dónde deben ir.

1.3.1. Declarar una salida en la sección publish:

La sección publish: dentro del bloque workflow declara qué salidas de procesos deben publicarse. Asigna la salida de SAMTOOLS_INDEX a un destino nombrado llamado bam_index.

genomics.nf
22
23
24
    publish:
    bam_index = SAMTOOLS_INDEX.out
}
genomics.nf
22
23
24
    publish:
    // Declare outputs to publish
}

A continuación, necesitaremos decirle a Nextflow dónde poner la salida publicada.

1.3.2. Configurar el destino de salida en el bloque output {}

El bloque output {} se encuentra fuera del workflow y especifica dónde se publica cada destino nombrado. Agreguemos un destino para bam_index que publique en un subdirectorio bam/.

genomics.nf
26
27
28
29
30
output {
    bam_index {
        path 'bam'
    }
}
genomics.nf
26
27
28
output {
    // Configure publish targets
}

Nota

Por defecto, Nextflow publica archivos de salida como enlaces simbólicos, lo que evita duplicación innecesaria. Aunque los archivos de datos que estamos usando aquí son muy pequeños, en genómica pueden volverse muy grandes. Los enlaces simbólicos se romperán cuando limpies tu directorio work, así que para workflows de producción es posible que quieras sobrescribir el modo de publicación predeterminado a 'copy'.

1.4. Ejecutar el workflow

En este punto, tenemos un workflow de indexación de un paso que debería ser completamente funcional. ¡Probemos que funciona!

Podemos ejecutarlo con -profile test para usar el valor predeterminado configurado en el perfil de prueba y evitar tener que escribir la ruta en la línea de comandos.

nextflow run genomics.nf -profile test
Salida del comando 
N E X T F L O W   ~  version 25.10.2

┃ Launching `genomics.nf` [reverent_sinoussi] DSL2 - revision: 41d43ad7fe

executor >  local (1)
[2a/e69536] SAMTOOLS_INDEX (1) | 1 of 1 ✔

Puedes verificar que el archivo de índice se haya generado correctamente mirando en el directorio de trabajo o en el directorio de resultados.

Contenido del directorio work 
work/2a/e695367b2f60df09cf826b07192dc3
├── reads_mother.bam -> /workspaces/training/nf4-science/genomics/data/bam/reads_mother.bam
└── reads_mother.bam.bai
Contenido del directorio de resultados 
results/
└── bam/
    └── reads_mother.bam.bai -> ...

¡Ahí está!

Conclusión

Sabes cómo crear un módulo que contiene un proceso, importarlo a un workflow, llamarlo con un canal de entrada y publicar los resultados.

¿Qué sigue?

Agregar un segundo paso que tome la salida del proceso de indexación y la use para ejecutar el llamado de variantes.

2. Agregar un segundo proceso para ejecutar GATK HaplotypeCaller en el archivo BAM indexado

Ahora que tenemos un índice para nuestro archivo de entrada, podemos pasar a configurar el paso de llamado de variantes.

Recuerda el comando gatk HaplotypeCaller de la Parte 1:

gatk HaplotypeCaller \
        -R /data/ref/ref.fasta \
        -I /data/bam/reads_mother.bam \
        -O reads_mother.vcf \
        -L /data/ref/intervals.bed

El comando toma un archivo BAM (-I), un genoma de referencia (-R) y un archivo de intervalos (-L), y produce un archivo VCF (-O) junto con su índice. La herramienta también espera que el índice del BAM, el índice de referencia y el diccionario de referencia estén ubicados junto a sus respectivos archivos. El URI del contenedor era community.wave.seqera.io/library/gatk4:4.5.0.0--730ee8817e436867.

Seguimos las mismas tres etapas que antes:

  1. Configurar las entradas
  2. Escribir el proceso de llamado de variantes y llamarlo en el workflow
  3. Configurar el manejo de la salida

2.1. Configurar las entradas

El paso de llamado de variantes requiere varios archivos de entrada adicionales. Necesitamos declarar parámetros para ellos, agregar valores predeterminados al perfil de prueba y crear variables para cargarlos.

2.1.1. Agregar declaraciones de parámetros para entradas accesorias

Dado que nuestro nuevo proceso espera un puñado de archivos adicionales, agrega declaraciones de parámetros para ellos en genomics.nf bajo la sección Pipeline parameters:

genomics.nf
 9
10
11
12
13
14
15
16
17
18
params {
    // Primary input
    input: Path

    // Accessory files
    reference: Path
    reference_index: Path
    reference_dict: Path
    intervals: Path
}
genomics.nf
 9
10
11
12
params {
    // Primary input
    input: Path
}

Como antes, proporcionaremos valores predeterminados a través del perfil de prueba en lugar de en línea.

2.1.2. Agregar valores predeterminados de archivos accesorios al perfil de prueba

Tal como hicimos para input en la sección 1.1.2, agrega valores predeterminados para los archivos accesorios al perfil de prueba en nextflow.config:

nextflow.config
 4
 5
 6
 7
 8
 9
10
test {
    params.input = "${projectDir}/data/bam/reads_mother.bam"
    params.reference = "${projectDir}/data/ref/ref.fasta"
    params.reference_index = "${projectDir}/data/ref/ref.fasta.fai"
    params.reference_dict = "${projectDir}/data/ref/ref.dict"
    params.intervals = "${projectDir}/data/ref/intervals.bed"
}
nextflow.config
4
5
6
test {
    params.input = "${projectDir}/data/bam/reads_mother.bam"
}

A continuación, necesitaremos crear variables que carguen estas rutas de archivos para usarlas en el workflow.

2.1.3. Crear variables para los archivos accesorios

Agrega variables para las rutas de archivos accesorios dentro del bloque workflow:

genomics.nf
21
22
23
24
25
26
27
28
29
30
31
32
33
34
workflow {

    main:
    // Create input channel (single file via CLI parameter)
    reads_ch = channel.fromPath(params.input)

    // Load the file paths for the accessory files (reference and intervals)
    ref_file        = file(params.reference)
    ref_index_file  = file(params.reference_index)
    ref_dict_file   = file(params.reference_dict)
    intervals_file  = file(params.intervals)

    // Create index file for input BAM file
    SAMTOOLS_INDEX(reads_ch)
genomics.nf
21
22
23
24
25
26
27
28
workflow {

    main:
    // Create input channel (single file via CLI parameter)
    reads_ch = channel.fromPath(params.input)

    // Create index file for input BAM file
    SAMTOOLS_INDEX(reads_ch)

La sintaxis file() le dice explícitamente a Nextflow que maneje estas entradas como rutas de archivos. Puedes aprender más sobre esto en la Misión Secundaria Trabajando con archivos.

2.2. Escribir el proceso de llamado de variantes y llamarlo en el workflow

Necesitamos escribir la definición del proceso en el archivo del módulo, importarlo al workflow usando una declaración include y llamarlo con las lecturas de entrada más la salida del paso de indexación y los archivos accesorios.

2.2.1. Completar el módulo para el proceso de llamado de variantes

Abre modules/gatk_haplotypecaller.nf y examina el esquema de la definición del proceso.

Adelante, completa la definición del proceso por ti mismo usando la información proporcionada arriba, luego verifica tu trabajo contra la solución en la pestaña "Después" a continuación.

modules/gatk_haplotypecaller.nf
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
#!/usr/bin/env nextflow

/*
 * Call variants with GATK HaplotypeCaller
 */
process GATK_HAPLOTYPECALLER {

    container

    input:

    output:

    script:
    """

    """
}
modules/gatk_haplotypecaller.nf
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#!/usr/bin/env nextflow

/*
 * Call variants with GATK HaplotypeCaller
 */
process GATK_HAPLOTYPECALLER {

    container "community.wave.seqera.io/library/gatk4:4.5.0.0--730ee8817e436867"

    input:
    path input_bam
    path input_bam_index
    path ref_fasta
    path ref_index
    path ref_dict
    path interval_list

    output:
    path "${input_bam}.vcf"     , emit: vcf
    path "${input_bam}.vcf.idx" , emit: idx

    script:
    """
    gatk HaplotypeCaller \
        -R ${ref_fasta} \
        -I ${input_bam} \
        -O ${input_bam}.vcf \
        -L ${interval_list}
    """
}

Notarás que este proceso tiene más entradas de las que el comando GATK en sí requiere. GATK sabe buscar el archivo de índice del BAM y los archivos accesorios del genoma de referencia basándose en convenciones de nomenclatura, pero Nextflow es agnóstico al dominio y no conoce estas convenciones. Necesitamos listarlos explícitamente para que Nextflow los prepare en el directorio de trabajo en tiempo de ejecución; de lo contrario, GATK arrojará un error sobre archivos faltantes.

De manera similar, listamos explícitamente el archivo de índice del VCF de salida ("${input_bam}.vcf.idx") para que Nextflow lo rastree para pasos posteriores. Usamos la sintaxis emit: para asignar un nombre a cada canal de salida, lo que será útil cuando conectemos las salidas al bloque publish.

Una vez que hayas completado esto, el proceso está completo. Para usarlo en el workflow, necesitarás importar el módulo y agregar un llamado al proceso.

2.2.2. Importar el nuevo módulo

Actualiza genomics.nf para importar el nuevo módulo:

genomics.nf
3
4
5
// Module INCLUDE statements
include { SAMTOOLS_INDEX } from './modules/samtools_index.nf'
include { GATK_HAPLOTYPECALLER } from './modules/gatk_haplotypecaller.nf'
genomics.nf
3
4
// Module INCLUDE statements
include { SAMTOOLS_INDEX } from './modules/samtools_index.nf'

El proceso ahora está disponible en el ámbito del workflow.

2.2.3. Agregar el llamado al proceso

Agrega el llamado al proceso en el cuerpo del workflow, bajo main::

genomics.nf
33
34
35
36
37
38
39
40
41
42
43
44
    // Create index file for input BAM file
    SAMTOOLS_INDEX(reads_ch)

    // Call variants from the indexed BAM file
    GATK_HAPLOTYPECALLER(
        reads_ch,
        SAMTOOLS_INDEX.out,
        ref_file,
        ref_index_file,
        ref_dict_file,
        intervals_file
    )
genomics.nf
33
34
    // Create index file for input BAM file
    SAMTOOLS_INDEX(reads_ch)

Deberías reconocer la sintaxis *.out de la serie de capacitación Hello Nextflow; le estamos diciendo a Nextflow que tome el canal de salida de SAMTOOLS_INDEX y lo conecte al llamado del proceso GATK_HAPLOTYPECALLER.

Nota

Observa que las entradas se proporcionan en el mismo orden exacto en el llamado al proceso que están listadas en el bloque input del proceso. En Nextflow, las entradas son posicionales, lo que significa que debes seguir el mismo orden; y por supuesto debe haber el mismo número de elementos.

2.3. Configurar el manejo de la salida

Necesitamos agregar las nuevas salidas a la declaración publish y configurar dónde van.

2.3.1. Agregar destinos de publicación para las salidas del llamado de variantes

Agrega las salidas VCF e índice a la sección publish::

genomics.nf
45
46
47
48
49
    publish:
    bam_index = SAMTOOLS_INDEX.out
    vcf = GATK_HAPLOTYPECALLER.out.vcf
    vcf_idx = GATK_HAPLOTYPECALLER.out.idx
}
genomics.nf
45
46
47
    publish:
    bam_index = SAMTOOLS_INDEX.out
}

A continuación, necesitaremos decirle a Nextflow dónde poner las nuevas salidas.

2.3.2. Configurar los nuevos destinos de salida

Agrega entradas para los destinos vcf y vcf_idx en el bloque output {}, publicando ambos en un subdirectorio vcf/:

genomics.nf
51
52
53
54
55
56
57
58
59
60
61
output {
    bam_index {
        path 'bam'
    }
    vcf {
        path 'vcf'
    }
    vcf_idx {
        path 'vcf'
    }
}
genomics.nf
49
50
51
52
53
output {
    bam_index {
        path 'bam'
    }
}

El VCF y su índice se publican como destinos separados que ambos van al subdirectorio vcf/.

2.4. Ejecutar el workflow

Ejecuta el workflow expandido, agregando -resume esta vez para que no tengamos que ejecutar el paso de indexación nuevamente.

nextflow run genomics.nf -profile test -resume
Salida del comando 
N E X T F L O W   ~  version 25.10.2

┃ Launching `genomics.nf` [grave_volta] DSL2 - revision: 4790abc96a

executor >  local (1)
[2a/e69536] SAMTOOLS_INDEX (1)       | 1 of 1, cached: 1 ✔
[53/e18e98] GATK_HAPLOTYPECALLER (1) | 1 of 1 ✔

Ahora si miramos la salida de la consola, vemos los dos procesos listados.

El primer proceso se omitió gracias al almacenamiento en caché, como se esperaba, mientras que el segundo proceso se ejecutó ya que es completamente nuevo.

Encontrarás los archivos de salida en el directorio de resultados (como enlaces simbólicos al directorio de trabajo).

Contenido del directorio 
results/
├── bam/
│   └── reads_mother.bam.bai -> ...
└── vcf/
    ├── reads_mother.bam.vcf -> ...
    └── reads_mother.bam.vcf.idx -> ...

Si abres el archivo VCF, deberías ver el mismo contenido que en el archivo que generaste al ejecutar el comando GATK directamente en el contenedor.

Contenido del archivo
reads_mother.bam.vcf
26
27
28
29
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	reads_mother
20_10037292_10066351	3480	.	C	CT	503.03	.	AC=2;AF=1.00;AN=2;DP=23;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=27.95;SOR=1.179	GT:AD:DP:GQ:PL	1/1:0,18:18:54:517,54,0
20_10037292_10066351	3520	.	AT	A	609.03	.	AC=2;AF=1.00;AN=2;DP=18;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=33.83;SOR=0.693	GT:AD:DP:GQ:PL	1/1:0,18:18:54:623,54,0
20_10037292_10066351	3529	.	T	A	155.64	.	AC=1;AF=0.500;AN=2;BaseQRankSum=-0.544;DP=21;ExcessHet=0.0000;FS=1.871;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=7.78;ReadPosRankSum=-1.158;SOR=1.034	GT:AD:DP:GQ:PL	0/1:12,8:20:99:163,0,328

Esta es la salida que nos interesa generar para cada muestra en nuestro estudio.

Conclusión

Sabes cómo hacer un workflow modular de dos pasos que hace trabajo de análisis real y es capaz de lidiar con idiosincrasias de formatos de archivos genómicos como los archivos accesorios.

¿Qué sigue?

Hacer que el workflow maneje múltiples muestras en lote.

3. Adaptar el workflow para ejecutarse en un lote de muestras

Está muy bien tener un workflow que pueda automatizar el procesamiento en una sola muestra, pero ¿qué pasa si tienes 1000 muestras? ¿Necesitas escribir un script bash que recorra todas tus muestras?

¡No, gracias a Dios! Solo haz un pequeño ajuste al código y Nextflow también manejará eso por ti.

3.1. Actualizar la entrada para listar tres muestras

Para ejecutar en múltiples muestras, actualiza el perfil de prueba para proporcionar un arreglo de rutas de archivos en lugar de una sola. Esta es una forma rápida de probar la ejecución de múltiples muestras; en el siguiente paso cambiaremos a un enfoque más escalable usando un archivo de entradas.

Primero, comenta la anotación de tipo en la declaración del parámetro, ya que los arreglos no pueden usar declaraciones tipadas:

genomics.nf
10
11
    // Primary input (array of three samples)
    input //: Path
genomics.nf
10
11
    // Primary input
    input: Path

Luego actualiza el perfil de prueba para listar las tres muestras:

nextflow.config
 4
 5
 6
 7
 8
 9
10
11
12
13
14
test {
    params.input = [
        "${projectDir}/data/bam/reads_mother.bam",
        "${projectDir}/data/bam/reads_father.bam",
        "${projectDir}/data/bam/reads_son.bam"
    ]
    params.reference = "${projectDir}/data/ref/ref.fasta"
    params.reference_index = "${projectDir}/data/ref/ref.fasta.fai"
    params.reference_dict = "${projectDir}/data/ref/ref.dict"
    params.intervals = "${projectDir}/data/ref/intervals.bed"
}
nextflow.config
 4
 5
 6
 7
 8
 9
10
test {
    params.input = "${projectDir}/data/bam/reads_mother.bam"
    params.reference = "${projectDir}/data/ref/ref.fasta"
    params.reference_index = "${projectDir}/data/ref/ref.fasta.fai"
    params.reference_dict = "${projectDir}/data/ref/ref.dict"
    params.intervals = "${projectDir}/data/ref/intervals.bed"
}

El channel factory en el cuerpo del workflow (.fromPath) acepta múltiples rutas de archivos tan bien como una sola, por lo que no se necesitan otros cambios.

3.2. Ejecutar el workflow

Intenta ejecutar el workflow ahora que la configuración está lista para ejecutarse en las tres muestras de prueba.

nextflow run genomics.nf -profile test -resume

Algo curioso: esto podría funcionar, O podría fallar. Por ejemplo, aquí hay una ejecución que tuvo éxito:

Salida del comando 
N E X T F L O W   ~  version 25.10.2

┃ Launching `genomics.nf` [peaceful_yalow] DSL2 - revision: a256d113ad

executor >  local (6)
[4f/7071b0] SAMTOOLS_INDEX (3)       | 3 of 3, cached: 1 ✔
[7a/89bc43] GATK_HAPLOTYPECALLER (2) | 3 of 3, cached: 1 ✔

Si tu ejecución del workflow tuvo éxito, ejecútalo nuevamente hasta que obtengas un error como este:

Salida del comando 
N E X T F L O W   ~  version 25.10.2

┃ Launching `genomics.nf` [loving_pasteur] DSL2 - revision: d2a8e63076

executor >  local (4)
[01/eea165] SAMTOOLS_INDEX (2)       | 3 of 3, cached: 1 ✔
[a5/fa9fd0] GATK_HAPLOTYPECALLER (3) | 1 of 3, cached: 1
ERROR ~ Error executing process > 'GATK_HAPLOTYPECALLER (2)'

Caused by:
  Process `GATK_HAPLOTYPECALLER (2)` terminated with an error exit status (2)

Command executed:

  gatk HaplotypeCaller         -R ref.fasta         -I reads_father.bam         -O reads_father.bam.vcf         -L intervals.bed

Command exit status:
  2

Command error:
  ...
  A USER ERROR has occurred: Traversal by intervals was requested but some input files are not indexed.
  ...

Si miras la salida de error del comando GATK, habrá una línea como esta:

A USER ERROR has occurred: Traversal by intervals was requested but some input files are not indexed.

Bueno, eso es extraño, considerando que indexamos explícitamente los archivos BAM en el primer paso del workflow. ¿Podría haber algo mal con la configuración?

3.3. Solucionar el problema

Inspeccionaremos los directorios de trabajo y usaremos el operador view() para averiguar qué salió mal.

3.3.1. Verificar los directorios de trabajo para los llamados relevantes

Echa un vistazo dentro del directorio de trabajo para el llamado al proceso GATK_HAPLOTYPECALLER fallido listado en la salida de la consola.

Contenido del directorio 
work/a5/fa9fd0994b6beede5fb9ea073596c2
├── intervals.bed -> /workspaces/training/nf4-science/genomics/data/ref/intervals.bed
├── reads_father.bam.bai -> /workspaces/training/nf4-science/genomics/work/01/eea16597bd6e810fb4cf89e60f8c2d/reads_father.bam.bai
├── reads_son.bam -> /workspaces/training/nf4-science/genomics/data/bam/reads_son.bam
├── reads_son.bam.vcf
├── reads_son.bam.vcf.idx
├── ref.dict -> /workspaces/training/nf4-science/genomics/data/ref/ref.dict
├── ref.fasta -> /workspaces/training/nf4-science/genomics/data/ref/ref.fasta
└── ref.fasta.fai -> /workspaces/training/nf4-science/genomics/data/ref/ref.fasta.fai

Presta especial atención a los nombres del archivo BAM y el índice BAM que están listados en este directorio: reads_son.bam y reads_father.bam.bai.

¿Qué demonios? Nextflow ha preparado un archivo de índice en el directorio de trabajo de este llamado al proceso, pero es el incorrecto. ¿Cómo pudo haber pasado esto?

3.3.2. Usar el operador view() para inspeccionar el contenido del canal

Agrega estas dos líneas en el cuerpo del workflow antes del llamado al proceso GATK_HAPLOTYPECALLER para ver el contenido del canal:

genomics.nf
36
37
38
39
40
41
42
43
    SAMTOOLS_INDEX(reads_ch)

    // temporary diagnostics
    reads_ch.view()
    SAMTOOLS_INDEX.out.view()

    // Call variants from the indexed BAM file
    GATK_HAPLOTYPECALLER(
genomics.nf
36
37
38
39
    SAMTOOLS_INDEX(reads_ch)

    // Call variants from the indexed BAM file
    GATK_HAPLOTYPECALLER(

Luego ejecuta el comando del workflow nuevamente.

nextflow run genomics.nf -profile test

Una vez más, esto puede tener éxito o fallar. Aquí está la salida de los dos llamados a .view() para una ejecución fallida:

/workspaces/training/nf4-science/genomics/data/bam/reads_mother.bam
/workspaces/training/nf4-science/genomics/data/bam/reads_father.bam
/workspaces/training/nf4-science/genomics/data/bam/reads_son.bam
/workspaces/training/nf4-science/genomics/work/9c/53492e3518447b75363e1cd951be4b/reads_father.bam.bai
/workspaces/training/nf4-science/genomics/work/cc/37894fffdf6cc84c3b0b47f9b536b7/reads_son.bam.bai
/workspaces/training/nf4-science/genomics/work/4d/dff681a3d137ba7d9866e3d9307bd0/reads_mother.bam.bai

Las primeras tres líneas corresponden al canal de entrada y las segundas, al canal de salida. Puedes ver que los archivos BAM y los archivos de índice para las tres muestras ¡no están listados en el mismo orden!

Nota

Cuando llamas a un proceso de Nextflow en un canal que contiene múltiples elementos, Nextflow intentará paralelizar la ejecución tanto como sea posible, y recolectará las salidas en el orden en que estén disponibles. La consecuencia es que las salidas correspondientes pueden recolectarse en un orden diferente al que se proporcionaron las entradas originales.

Como está escrito actualmente, nuestro script de workflow asume que los archivos de índice saldrán del paso de indexación listados en el mismo orden madre/padre/hijo que se dieron las entradas. Pero eso no está garantizado que sea el caso, por lo que a veces (aunque no siempre) los archivos incorrectos se emparejan en el segundo paso.

Para arreglar esto, necesitamos asegurarnos de que los archivos BAM y sus archivos de índice viajen juntos a través de los canales.

Consejo

Las declaraciones view() en el código del workflow no hacen nada, por lo que no es un problema dejarlas. Sin embargo, desordenarán tu salida de consola, por lo que recomendamos eliminarlas cuando hayas terminado de solucionar el problema.

3.4. Actualizar el workflow para manejar los archivos de índice correctamente

La solución es agrupar cada archivo BAM con su índice en una tupla, luego actualizar el proceso descendente y la configuración del workflow para que coincidan.

3.4.1. Cambiar la salida del módulo SAMTOOLS_INDEX a una tupla

La forma más simple de asegurar que un archivo BAM y su índice permanezcan estrechamente asociados es empaquetarlos juntos en una tupla que salga de la tarea de índice.

Nota

Una tupla es una lista finita y ordenada de elementos que se usa comúnmente para devolver múltiples valores de una función. Las tuplas son particularmente útiles para pasar múltiples entradas o salidas entre procesos mientras se preserva su asociación y orden.

Actualiza la salida en modules/samtools_index.nf para incluir el archivo BAM:

modules/samtools_index.nf
14
15
    output:
    tuple path(input_bam), path("${input_bam}.bai")
modules/samtools_index.nf
14
15
    output:
    path "${input_bam}.bai"

De esta manera, cada archivo de índice estará estrechamente acoplado con su archivo BAM original, y la salida general del paso de indexación será un solo canal que contiene pares de archivos.

3.4.2. Cambiar la entrada del módulo GATK_HAPLOTYPECALLER para aceptar una tupla

Dado que hemos cambiado la 'forma' de la salida del primer proceso, necesitamos actualizar la definición de entrada del segundo proceso para que coincida.

Actualiza modules/gatk_haplotypecaller.nf:

modules/gatk_haplotypecaller.nf
11
12
    input:
    tuple path(input_bam), path(input_bam_index)
modules/gatk_haplotypecaller.nf
11
12
13
    input:
    path input_bam
    path input_bam_index

A continuación, necesitaremos actualizar el workflow para reflejar la nueva estructura de tupla en el llamado al proceso y los destinos de publicación.

3.4.3. Actualizar el llamado a GATK_HAPLOTYPECALLER en el workflow

Ya no necesitamos proporcionar el reads_ch original al proceso GATK_HAPLOTYPECALLER, ya que el archivo BAM ahora está incluido en el canal de salida de SAMTOOLS_INDEX.

Actualiza el llamado en genomics.nf:

genomics.nf
42
43
    GATK_HAPLOTYPECALLER(
        SAMTOOLS_INDEX.out,
genomics.nf
42
43
44
    GATK_HAPLOTYPECALLER(
        reads_ch,
        SAMTOOLS_INDEX.out,

Finalmente, necesitamos actualizar los destinos de publicación para reflejar la nueva estructura de salida.

3.4.4. Actualizar el destino de publicación para la salida del BAM indexado

Dado que la salida de SAMTOOLS_INDEX ahora es una tupla que contiene tanto el archivo BAM como su índice, renombra el destino de publicación de bam_index a indexed_bam para reflejar mejor su contenido:

genomics.nf
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
    publish:
    indexed_bam = SAMTOOLS_INDEX.out
    vcf = GATK_HAPLOTYPECALLER.out.vcf
    vcf_idx = GATK_HAPLOTYPECALLER.out.idx
}

output {
    indexed_bam {
        path 'bam'
    }
    vcf {
        path 'vcf'
    }
    vcf_idx {
        path 'vcf'
    }
}
genomics.nf
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
    publish:
    bam_index = SAMTOOLS_INDEX.out
    vcf = GATK_HAPLOTYPECALLER.out.vcf
    vcf_idx = GATK_HAPLOTYPECALLER.out.idx
}

output {
    bam_index {
        path 'bam'
    }
    vcf {
        path 'vcf'
    }
    vcf_idx {
        path 'vcf'
    }
}

Con estos cambios, el BAM y su índice están garantizados de viajar juntos, por lo que el emparejamiento siempre será correcto.

3.5. Ejecutar el workflow corregido

Ejecuta el workflow nuevamente para asegurarte de que esto funcionará de manera confiable en el futuro.

nextflow run genomics.nf -profile test

Esta vez (y cada vez) todo debería ejecutarse correctamente:

Salida del comando 
N E X T F L O W   ~  version 25.10.2

┃ Launching `genomics.nf` [special_goldstine] DSL2 - revision: 4cbbf6ea3e

executor >  local (6)
[d6/10c2c4] SAMTOOLS_INDEX (1)       | 3 of 3 ✔
[88/1783aa] GATK_HAPLOTYPECALLER (2) | 3 of 3 ✔

El directorio de resultados ahora contiene tanto archivos BAM como BAI para cada muestra (de la tupla), junto con las salidas VCF:

Contenido del directorio de resultados 
results/
├── bam/
│   ├── reads_father.bam -> ...
│   ├── reads_father.bam.bai -> ...
│   ├── reads_mother.bam -> ...
│   ├── reads_mother.bam.bai -> ...
│   ├── reads_son.bam -> ...
│   └── reads_son.bam.bai -> ...
└── vcf/
    ├── reads_father.bam.vcf -> ...
    ├── reads_father.bam.vcf.idx -> ...
    ├── reads_mother.bam.vcf -> ...
    ├── reads_mother.bam.vcf.idx -> ...
    ├── reads_son.bam.vcf -> ...
    └── reads_son.bam.vcf.idx -> ...

Al agrupar archivos asociados en tuplas, nos aseguramos de que los archivos correctos siempre viajen juntos a través del workflow. El workflow ahora procesa cualquier número de muestras de manera confiable, pero listarlas individualmente en la configuración no es muy escalable. En el siguiente paso, cambiaremos a leer entradas desde un archivo.

Conclusión

Sabes cómo hacer que tu workflow se ejecute en múltiples muestras (independientemente).

¿Qué sigue?

Hacer más fácil manejar muestras en lote.

4. Hacer que el workflow acepte un archivo de texto que contenga un lote de archivos de entrada

Una forma muy común de proporcionar múltiples archivos de datos de entrada a un workflow es hacerlo con un archivo de texto que contenga las rutas de los archivos. Puede ser tan simple como un archivo de texto que liste una ruta de archivo por línea y nada más, o el archivo puede contener metadatos adicionales, en cuyo caso a menudo se llama samplesheet.

Aquí vamos a mostrarte cómo hacer el caso simple.

4.1. Examinar el archivo CSV proporcionado que lista las rutas de archivos de entrada

Ya hicimos un archivo CSV que lista las rutas de archivos de entrada, llamado samplesheet.csv, que puedes encontrar en el directorio data/.

samplesheet.csv
sample_id,reads_bam
NA12878,/workspaces/training/nf4-science/genomics/data/bam/reads_mother.bam
NA12877,/workspaces/training/nf4-science/genomics/data/bam/reads_father.bam
NA12882,/workspaces/training/nf4-science/genomics/data/bam/reads_son.bam

Este archivo CSV incluye una línea de encabezado que nombra las columnas.

Advertencia

Las rutas de archivos en el CSV son rutas absolutas que deben coincidir con tu entorno. Si no estás ejecutando esto en el entorno de capacitación que proporcionamos, necesitarás actualizar las rutas para que coincidan con tu sistema.

4.2. Actualizar el parámetro y el perfil de prueba

Cambia el parámetro input para que apunte al archivo samplesheet.csv en lugar de listar muestras individuales.

Restaura la anotación de tipo en el bloque params (ya que es una sola ruta nuevamente):

genomics.nf
10
11
    // Primary input (file of input files, one per line)
    input: Path
genomics.nf
10
11
    // Primary input (array of three samples)
    input

Luego actualiza el perfil de prueba para que apunte al archivo de texto:

nextflow.config
 4
 5
 6
 7
 8
 9
10
test {
    params.input = "${projectDir}/data/samplesheet.csv"
    params.reference = "${projectDir}/data/ref/ref.fasta"
    params.reference_index = "${projectDir}/data/ref/ref.fasta.fai"
    params.reference_dict = "${projectDir}/data/ref/ref.dict"
    params.intervals = "${projectDir}/data/ref/intervals.bed"
}
nextflow.config
 4
 5
 6
 7
 8
 9
10
11
12
13
14
test {
    params.input = [
        "${projectDir}/data/bam/reads_mother.bam",
        "${projectDir}/data/bam/reads_father.bam",
        "${projectDir}/data/bam/reads_son.bam"
    ]
    params.reference = "${projectDir}/data/ref/ref.fasta"
    params.reference_index = "${projectDir}/data/ref/ref.fasta.fai"
    params.reference_dict = "${projectDir}/data/ref/ref.dict"
    params.intervals = "${projectDir}/data/ref/intervals.bed"
}

La lista de archivos ya no vive en el código en absoluto, lo cual es un gran paso en la dirección correcta.

4.3. Actualizar el channel factory para analizar entrada CSV

Actualmente, nuestro channel factory de entrada trata cualquier archivo que le demos como las entradas de datos que queremos alimentar al proceso de indexación. Dado que ahora le estamos dando un archivo que lista rutas de archivos de entrada, necesitamos cambiar su comportamiento para analizar el archivo y tratar las rutas de archivos que contiene como las entradas de datos.

Podemos hacer esto usando el mismo patrón que usamos en la Parte 2 de Hello Nextflow: aplicando el operador splitCsv() para analizar el archivo, luego una operación map para extraer la ruta del archivo de cada fila.

genomics.nf
24
25
26
27
    // Create input channel from a CSV file listing input file paths
    reads_ch = channel.fromPath(params.input)
            .splitCsv(header: true)
            .map { row -> file(row.reads_bam) }
genomics.nf
24
25
    // Create input channel (single file via CLI parameter)
    reads_ch = channel.fromPath(params.input)

Una cosa que es nueva comparada con lo que encontraste en el curso Hello Nextflow es que este CSV tiene una línea de encabezado, por lo que agregamos header: true al llamado splitCsv(). Eso nos permite referenciar columnas por nombre en la operación map: row.reads_bam extrae la ruta del archivo de la columna reads_bam de cada fila.

Consejo

Si no estás seguro de entender qué están haciendo los operadores aquí, esta es otra gran oportunidad para usar el operador .view() para ver cómo se ven los contenidos del canal antes y después de aplicarlos.

4.4. Ejecutar el workflow

Ejecuta el workflow una vez más.

nextflow run genomics.nf -profile test
Salida del comando 
N E X T F L O W   ~  version 25.10.2

┃ Launching `genomics.nf` [sick_albattani] DSL2 - revision: 46d84642f6

executor >  local (6)
[18/23b4bb] SAMTOOLS_INDEX (1)       | 3 of 3 ✔
[12/f727bb] GATK_HAPLOTYPECALLER (3) | 3 of 3 ✔

Esto debería producir el mismo resultado que antes. Nuestro workflow simple de llamado de variantes ahora tiene todas las características básicas que queríamos.

Conclusión

Sabes cómo hacer un workflow modular de múltiples pasos para indexar un archivo BAM y aplicar llamado de variantes por muestra usando GATK.

Más generalmente, has aprendido cómo usar componentes y lógica esenciales de Nextflow para construir un pipeline de genómica simple que hace trabajo real, teniendo en cuenta las idiosincrasias de los formatos de archivos genómicos y los requisitos de las herramientas.

¿Qué sigue?

¡Toma un descanso! Eso fue mucho.

Cuando te sientas renovado, dirígete a la Parte 3, donde aprenderás cómo transformar este workflow simple de llamado de variantes por muestra para aplicar llamado de variantes conjunto a los datos.