Vai al contenuto

Parte 2: Variant calling per campione

Traduzione assistita da IA - scopri di più e suggerisci miglioramenti

Nella Parte 1, avete testato manualmente i comandi di Samtools e GATK nei rispettivi container. Successivamente, andremo a integrare quegli stessi comandi in un flusso di lavoro Nextflow.

Obiettivo

In questa parte del corso, svilupperemo un flusso di lavoro che esegue le seguenti operazioni:

  1. Generare un file di indice per ogni file BAM di input utilizzando Samtools
  2. Eseguire GATK HaplotypeCaller su ogni file BAM di input per generare chiamate di varianti per campione in formato VCF (Variant Call Format)
BAMSamtools indexBAM indexIntervalsReference+ index & dictGATK HaplotypeCallerVCF + index

Questo automatizza i passaggi della prima sezione della Parte 1: Panoramica del metodo, dove avete eseguito questi comandi manualmente nei container.

Come punto di partenza, vi forniamo un file di flusso di lavoro, genomics.nf, che delinea le parti principali del flusso di lavoro, insieme a due file modulo, samtools_index.nf e gatk_haplotypecaller.nf, che delineano la struttura di ogni processo.

File di struttura
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:
    """

    """
}

Questi file non sono funzionali; il loro scopo è solo quello di servire da struttura da riempire con le parti interessanti del codice.

Piano della lezione

Per rendere il processo di sviluppo più educativo, abbiamo suddiviso tutto in quattro fasi:

  1. Scrivere un flusso di lavoro a singolo stadio che esegue Samtools index su un file BAM. Questo copre la creazione di un modulo, la sua importazione e la sua chiamata in un flusso di lavoro.
  2. Aggiungere un secondo processo per eseguire GATK HaplotypeCaller sul file BAM indicizzato. Questo introduce il concatenamento degli output dei processi agli input e la gestione dei file accessori.
  3. Adattare il flusso di lavoro per eseguirlo su un batch di campioni. Questo copre l'esecuzione parallela e introduce le tuple per mantenere insieme i file associati.
  4. Far accettare al flusso di lavoro un file di testo contenente un batch di file di input. Questo dimostra un pattern comune per fornire input in blocco.

Ogni fase si concentra su un aspetto specifico dello sviluppo del flusso di lavoro.

Suggerimento

Assicuratevi di essere nella directory di lavoro corretta: cd /workspaces/training/nf4-science/genomics

1. Scrivere un flusso di lavoro a singolo stadio che esegue Samtools index su un file BAM

Questa prima fase si concentra sulle basi: caricare un file BAM e generare un indice per esso.

Ricordate il comando samtools index dalla Parte 1:

samtools index '<input_bam>'

Il comando prende un file BAM come input e produce un file di indice .bai accanto ad esso. L'URI del container era community.wave.seqera.io/library/samtools:1.20--b5dfbd93de237464.

Prenderemo queste informazioni e le integreremo in Nextflow in tre fasi:

  1. Configurare l'input
  2. Scrivere il processo di indicizzazione e chiamarlo nel flusso di lavoro
  3. Configurare la gestione dell'output

1.1. Configurare l'input

Dobbiamo dichiarare un parametro di input, creare un profilo di test per fornire un valore predefinito conveniente e creare un canale di input.

1.1.1. Aggiungere una dichiarazione di parametro di input

Nel file principale del flusso di lavoro genomics.nf, sotto la sezione Pipeline parameters, dichiarate un parametro CLI chiamato 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

Questo configura il parametro CLI, ma non vogliamo digitare il percorso del file ogni volta che eseguiamo il flusso di lavoro durante lo sviluppo. Ci sono diverse opzioni per fornire un valore predefinito; qui utilizziamo un profilo di test.

1.1.2. Creare un profilo di test con un valore predefinito in nextflow.config

Un profilo di test fornisce valori predefiniti convenienti per provare un flusso di lavoro senza specificare input dalla riga di comando. Questa è una convenzione comune nell'ecosistema Nextflow (vedi Hello Config per maggiori dettagli).

Aggiungete un blocco profiles a nextflow.config con un profilo test che imposta il parametro input su uno dei file BAM di test.

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

Qui stiamo usando ${projectDir}, una variabile integrata di Nextflow che punta alla directory dove si trova lo script del flusso di lavoro. Questo rende facile fare riferimento a file di dati e altre risorse senza codificare percorsi assoluti.

1.1.3. Configurare il canale di input

Nel blocco workflow, create un canale di input dal valore del parametro utilizzando la fabbrica di canali .fromPath (come usato in 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

Successivamente, dovremo creare il processo per eseguire l'indicizzazione su questo input.

1.2. Scrivere il processo di indicizzazione e chiamarlo nel flusso di lavoro

Dobbiamo scrivere la definizione del processo nel file modulo, importarlo nel flusso di lavoro usando un'istruzione include e chiamarlo sull'input.

1.2.1. Completare il modulo per il processo di indicizzazione

Aprite modules/samtools_index.nf ed esaminate la struttura della definizione del processo. Dovreste riconoscere gli elementi strutturali principali; in caso contrario, considerate di leggere Hello Nextflow per un ripasso.

Procedete e completate la definizione del processo da soli utilizzando le informazioni fornite sopra, poi verificate il vostro lavoro confrontandolo con la soluzione nella scheda "Dopo" qui sotto.

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 volta completato questo, il processo è completo. Per usarlo nel flusso di lavoro, dovrete importare il modulo e aggiungere una chiamata al processo.

1.2.2. Includere il modulo

In genomics.nf, aggiungete un'istruzione include per rendere il processo disponibile al flusso di lavoro:

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

Il processo è ora disponibile nello scope del flusso di lavoro.

1.2.3. Chiamare il processo di indicizzazione sull'input

Ora aggiungiamo una chiamata a SAMTOOLS_INDEX nel blocco workflow, passando il canale di input come argomento.

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

Il flusso di lavoro ora carica l'input ed esegue il processo di indicizzazione su di esso. Successivamente, dobbiamo configurare come viene pubblicato l'output.

1.3. Configurare la gestione dell'output

Dobbiamo dichiarare quali output del processo pubblicare e specificare dove devono andare.

1.3.1. Dichiarare un output nella sezione publish:

La sezione publish: all'interno del blocco workflow dichiara quali output del processo devono essere pubblicati. Assegnate l'output di SAMTOOLS_INDEX a un target denominato bam_index.

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

Successivamente, dovremo dire a Nextflow dove mettere l'output pubblicato.

1.3.2. Configurare il target di output nel blocco output {}

Il blocco output {} si trova fuori dal flusso di lavoro e specifica dove viene pubblicato ogni target denominato. Aggiungiamo un target per bam_index che pubblica in una sottodirectory bam/.

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

Nota

Per impostazione predefinita, Nextflow pubblica i file di output come link simbolici, il che evita duplicazioni non necessarie. Anche se i file di dati che stiamo usando qui sono molto piccoli, in genomica possono diventare molto grandi. I link simbolici si romperanno quando pulite la vostra directory work, quindi per flussi di lavoro in produzione potreste voler sovrascrivere la modalità di pubblicazione predefinita con 'copy'.

1.4. Eseguire il flusso di lavoro

A questo punto, abbiamo un flusso di lavoro di indicizzazione a un passaggio che dovrebbe essere completamente funzionale. Verifichiamo che funzioni!

Possiamo eseguirlo con -profile test per usare il valore predefinito configurato nel profilo di test ed evitare di dover scrivere il percorso sulla riga di comando.

nextflow run genomics.nf -profile test
Output 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 ✔

Potete verificare che il file di indice sia stato generato correttamente guardando nella directory di lavoro o nella directory dei risultati.

Contenuti della directory di lavoro 
work/2a/e695367b2f60df09cf826b07192dc3
├── reads_mother.bam -> /workspaces/training/nf4-science/genomics/data/bam/reads_mother.bam
└── reads_mother.bam.bai
Contenuti della directory dei risultati 
results/
└── bam/
    └── reads_mother.bam.bai -> ...

Eccolo!

Takeaway

Sapete come creare un modulo contenente un processo, importarlo in un flusso di lavoro, chiamarlo con un canale di input e pubblicare i risultati.

Cosa c'è dopo?

Aggiungere un secondo passaggio che prende l'output del processo di indicizzazione e lo usa per eseguire il variant calling.

2. Aggiungere un secondo processo per eseguire GATK HaplotypeCaller sul file BAM indicizzato

Ora che abbiamo un indice per il nostro file di input, possiamo passare alla configurazione del passaggio di variant calling.

Ricordate il comando gatk HaplotypeCaller dalla Parte 1:

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

Il comando prende un file BAM (-I), un genoma di riferimento (-R) e un file di intervalli (-L), e produce un file VCF (-O) insieme al suo indice. Lo strumento si aspetta anche che l'indice BAM, l'indice di riferimento e il dizionario di riferimento siano co-localizzati con i rispettivi file. L'URI del container era community.wave.seqera.io/library/gatk4:4.5.0.0--730ee8817e436867.

Seguiamo le stesse tre fasi di prima:

  1. Configurare gli input
  2. Scrivere il processo di variant calling e chiamarlo nel flusso di lavoro
  3. Configurare la gestione dell'output

2.1. Configurare gli input

Il passaggio di variant calling richiede diversi file di input aggiuntivi. Dobbiamo dichiarare parametri per essi, aggiungere valori predefiniti al profilo di test e creare variabili per caricarli.

2.1.1. Aggiungere dichiarazioni di parametri per gli input accessori

Poiché il nostro nuovo processo si aspetta una manciata di file aggiuntivi da fornire, aggiungete dichiarazioni di parametri per essi in genomics.nf sotto la sezione 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
}

Come prima, forniremo valori predefiniti attraverso il profilo di test piuttosto che inline.

2.1.2. Aggiungere valori predefiniti per i file accessori al profilo di test

Proprio come abbiamo fatto per input nella sezione 1.1.2, aggiungete valori predefiniti per i file accessori al profilo di test in 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"
}

Successivamente, dovremo creare variabili che caricano questi percorsi di file per l'uso nel flusso di lavoro.

2.1.3. Creare variabili per i file accessori

Aggiungete variabili per i percorsi dei file accessori all'interno del blocco 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 sintassi file() dice esplicitamente a Nextflow di gestire questi input come percorsi di file. Potete saperne di più su questo nella Side Quest Working with files.

2.2. Scrivere il processo di variant calling e chiamarlo nel flusso di lavoro

Dobbiamo scrivere la definizione del processo nel file modulo, importarlo nel flusso di lavoro usando un'istruzione include e chiamarlo sull'input reads più l'output del passaggio di indicizzazione e i file accessori.

2.2.1. Completare il modulo per il processo di variant calling

Aprite modules/gatk_haplotypecaller.nf ed esaminate la struttura della definizione del processo.

Procedete e completate la definizione del processo da soli utilizzando le informazioni fornite sopra, poi verificate il vostro lavoro confrontandolo con la soluzione nella scheda "Dopo" qui sotto.

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}
    """
}

Noterete che questo processo ha più input di quelli che il comando GATK stesso richiede. GATK sa cercare il file di indice BAM e i file accessori del genoma di riferimento in base alle convenzioni di denominazione, ma Nextflow è domain-agnostic e non conosce queste convenzioni. Dobbiamo elencarli esplicitamente in modo che Nextflow li prepari nella directory di lavoro al momento dell'esecuzione; altrimenti GATK genererà un errore sui file mancanti.

Allo stesso modo, elenchiamo esplicitamente il file di indice dell'output VCF ("${input_bam}.vcf.idx") in modo che Nextflow ne tenga traccia per i passaggi successivi. Usiamo la sintassi emit: per assegnare un nome a ogni canale di output, che diventerà utile quando collegheremo gli output nel blocco publish.

Una volta completato questo, il processo è completo. Per usarlo nel flusso di lavoro, dovrete importare il modulo e aggiungere una chiamata al processo.

2.2.2. Importare il nuovo modulo

Aggiornate genomics.nf per importare il nuovo modulo:

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'

Il processo è ora disponibile nello scope del flusso di lavoro.

2.2.3. Aggiungere la chiamata al processo

Aggiungete la chiamata al processo nel corpo del workflow, sotto 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)

Dovreste riconoscere la sintassi *.out dalla serie di formazione Hello Nextflow; stiamo dicendo a Nextflow di prendere il canale in output da SAMTOOLS_INDEX e collegarlo alla chiamata del processo GATK_HAPLOTYPECALLER.

Nota

Notate che gli input sono forniti nello stesso identico ordine nella chiamata al processo come sono elencati nel blocco input del processo. In Nextflow, gli input sono posizionali, il che significa che dovete seguire lo stesso ordine; e ovviamente ci deve essere lo stesso numero di elementi.

2.3. Configurare la gestione dell'output

Dobbiamo aggiungere i nuovi output alla dichiarazione publish e configurare dove vanno.

2.3.1. Aggiungere target di pubblicazione per gli output del variant calling

Aggiungete gli output VCF e indice alla sezione 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
}

Successivamente, dovremo dire a Nextflow dove mettere i nuovi output.

2.3.2. Configurare i nuovi target di output

Aggiungete voci per i target vcf e vcf_idx nel blocco output {}, pubblicando entrambi in una sottodirectory 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'
    }
}

Il VCF e il suo indice sono pubblicati come target separati che vanno entrambi nella sottodirectory vcf/.

2.4. Eseguire il flusso di lavoro

Eseguite il flusso di lavoro espanso, aggiungendo -resume questa volta in modo da non dover eseguire nuovamente il passaggio di indicizzazione.

nextflow run genomics.nf -profile test -resume
Output 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 ✔

Ora se guardiamo l'output della console, vediamo i due processi elencati.

Il primo processo è stato saltato grazie alla cache, come previsto, mentre il secondo processo è stato eseguito poiché è completamente nuovo.

Troverete i file di output nella directory dei risultati (come link simbolici alla directory di lavoro).

Contenuto della directory 
results/
├── bam/
│   └── reads_mother.bam.bai -> ...
└── vcf/
    ├── reads_mother.bam.vcf -> ...
    └── reads_mother.bam.vcf.idx -> ...

Se aprite il file VCF, dovreste vedere lo stesso contenuto del file che avete generato eseguendo il comando GATK direttamente nel container.

Contenuto del file
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

Questo è l'output che ci interessa generare per ogni campione nel nostro studio.

Takeaway

Sapete come creare un flusso di lavoro modulare a due passaggi che esegue un lavoro di analisi reale ed è in grado di gestire le idiosincrasie dei formati di file genomici come i file accessori.

Cosa c'è dopo?

Far gestire al flusso di lavoro più campioni in blocco.

3. Adattare il flusso di lavoro per eseguirlo su un batch di campioni

Va benissimo avere un flusso di lavoro che può automatizzare l'elaborazione su un singolo campione, ma cosa succede se avete 1000 campioni? Dovete scrivere uno script bash che cicla attraverso tutti i vostri campioni?

No, per fortuna! Basta fare una piccola modifica al codice e Nextflow gestirà anche questo per voi.

3.1. Aggiornare l'input per elencare tre campioni

Per eseguire su più campioni, aggiornate il profilo di test per fornire un array di percorsi di file invece di uno singolo. Questo è un modo rapido per testare l'esecuzione multi-campione; nel prossimo passaggio passeremo a un approccio più scalabile usando un file di input.

Prima, commentate l'annotazione di tipo nella dichiarazione del parametro, poiché gli array non possono usare dichiarazioni tipizzate:

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

Poi aggiornate il profilo di test per elencare tutti e tre i campioni:

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"
}

La fabbrica di canali nel corpo del workflow (.fromPath) accetta più percorsi di file così come uno singolo, quindi non sono necessarie altre modifiche.

3.2. Eseguire il flusso di lavoro

Provate a eseguire il flusso di lavoro ora che l'impianto è configurato per eseguire su tutti e tre i campioni di test.

nextflow run genomics.nf -profile test -resume

Cosa curiosa: questo potrebbe funzionare, OPPURE potrebbe fallire. Ad esempio, ecco un'esecuzione che ha avuto successo:

Output 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 ✔

Se l'esecuzione del vostro flusso di lavoro ha avuto successo, eseguitelo di nuovo finché non ottenete un errore come questo:

Output 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.
  ...

Se guardate l'output di errore del comando GATK, ci sarà una riga come questa:

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

Beh, è strano, considerando che abbiamo esplicitamente indicizzato i file BAM nel primo passaggio del flusso di lavoro. Potrebbe esserci qualcosa di sbagliato nell'impianto?

3.3. Risolvere il problema

Ispezioniamo le directory di lavoro e usiamo l'operatore view() per capire cosa è andato storto.

3.3.1. Controllare le directory di lavoro per le chiamate rilevanti

Date un'occhiata all'interno della directory di lavoro per la chiamata al processo GATK_HAPLOTYPECALLER fallita elencata nell'output della console.

Contenuto della directory 
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

Prestate particolare attenzione ai nomi del file BAM e dell'indice BAM che sono elencati in questa directory: reads_son.bam e reads_father.bam.bai.

Che diavolo? Nextflow ha preparato un file di indice nella directory di lavoro di questa chiamata al processo, ma è quello sbagliato. Come può essere successo?

3.3.2. Usare l'operatore view() per ispezionare il contenuto dei canali

Aggiungete queste due righe nel corpo del workflow prima della chiamata al processo GATK_HAPLOTYPECALLER per visualizzare il contenuto del canale:

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(

Poi eseguite di nuovo il comando del flusso di lavoro.

nextflow run genomics.nf -profile test

Ancora una volta, questo può avere successo o fallire. Ecco come appare l'output delle due chiamate .view() per un'esecuzione fallita:

/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

Le prime tre righe corrispondono al canale di input e le seconde al canale di output. Potete vedere che i file BAM e i file di indice per i tre campioni non sono elencati nello stesso ordine!

Nota

Quando chiamate un processo Nextflow su un canale contenente più elementi, Nextflow cercherà di parallelizzare l'esecuzione il più possibile e raccoglierà gli output in qualsiasi ordine diventino disponibili. La conseguenza è che gli output corrispondenti possono essere raccolti in un ordine diverso da quello in cui sono stati forniti gli input originali.

Come attualmente scritto, il nostro script del flusso di lavoro presume che i file di indice usciranno dal passaggio di indicizzazione elencati nello stesso ordine madre/padre/figlio in cui sono stati forniti gli input. Ma questo non è garantito, motivo per cui a volte (anche se non sempre) i file sbagliati vengono accoppiati nel secondo passaggio.

Per risolvere questo, dobbiamo assicurarci che i file BAM e i loro file di indice viaggino insieme attraverso i canali.

Suggerimento

Le istruzioni view() nel codice del flusso di lavoro non fanno nulla, quindi non è un problema lasciarle. Tuttavia ingombreranno l'output della console, quindi raccomandiamo di rimuoverle quando avete finito di risolvere il problema.

3.4. Aggiornare il flusso di lavoro per gestire correttamente i file di indice

La soluzione è raggruppare ogni file BAM con il suo indice in una tupla, poi aggiornare il processo downstream e l'impianto del flusso di lavoro per corrispondere.

3.4.1. Cambiare l'output del modulo SAMTOOLS_INDEX in una tupla

Il modo più semplice per garantire che un file BAM e il suo indice rimangano strettamente associati è impacchettarli insieme in una tupla in uscita dall'attività di indice.

Nota

Una tupla è una lista finita e ordinata di elementi che è comunemente usata per restituire più valori da una funzione. Le tuple sono particolarmente utili per passare più input o output tra processi preservando la loro associazione e ordine.

Aggiornate l'output in modules/samtools_index.nf per includere il file 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"

In questo modo, ogni file di indice sarà strettamente accoppiato con il suo file BAM originale, e l'output complessivo del passaggio di indicizzazione sarà un singolo canale contenente coppie di file.

3.4.2. Cambiare l'input del modulo GATK_HAPLOTYPECALLER per accettare una tupla

Poiché abbiamo cambiato la 'forma' dell'output del primo processo, dobbiamo aggiornare la definizione dell'input del secondo processo per corrispondere.

Aggiornate 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

Successivamente, dovremo aggiornare il flusso di lavoro per riflettere la nuova struttura della tupla nella chiamata al processo e nei target di pubblicazione.

3.4.3. Aggiornare la chiamata a GATK_HAPLOTYPECALLER nel flusso di lavoro

Non abbiamo più bisogno di fornire il reads_ch originale al processo GATK_HAPLOTYPECALLER, poiché il file BAM è ora raggruppato nel canale di output da SAMTOOLS_INDEX.

Aggiornate la chiamata in genomics.nf:

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

Infine, dobbiamo aggiornare i target di pubblicazione per riflettere la nuova struttura di output.

3.4.4. Aggiornare il target di pubblicazione per l'output BAM indicizzato

Poiché l'output di SAMTOOLS_INDEX è ora una tupla contenente sia il file BAM che il suo indice, rinominate il target di pubblicazione da bam_index a indexed_bam per riflettere meglio il suo contenuto:

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 queste modifiche, il BAM e il suo indice sono garantiti viaggiare insieme, quindi l'accoppiamento sarà sempre corretto.

3.5. Eseguire il flusso di lavoro corretto

Eseguite di nuovo il flusso di lavoro per assicurarvi che funzioni in modo affidabile andando avanti.

nextflow run genomics.nf -profile test

Questa volta (e ogni volta) tutto dovrebbe funzionare correttamente:

Output 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 ✔

La directory dei risultati ora contiene sia file BAM che BAI per ogni campione (dalla tupla), insieme agli output VCF:

Contenuti della directory dei risultati 
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 -> ...

Raggruppando i file associati in tuple, abbiamo garantito che i file corretti viaggino sempre insieme attraverso il flusso di lavoro. Il flusso di lavoro ora elabora qualsiasi numero di campioni in modo affidabile, ma elencarli individualmente nella configurazione non è molto scalabile. Nel prossimo passaggio, passeremo alla lettura degli input da un file.

Takeaway

Sapete come far eseguire il vostro flusso di lavoro su più campioni (indipendentemente).

Cosa c'è dopo?

Rendere più facile gestire i campioni in blocco.

4. Far accettare al flusso di lavoro un file CSV contenente un batch di file di input

Un modo molto comune per fornire più file di dati di input a un flusso di lavoro è farlo con un file di testo contenente i percorsi dei file. Può essere semplice come un file di testo che elenca un percorso di file per riga e nient'altro, oppure il file può contenere metadati aggiuntivi, nel qual caso è spesso chiamato samplesheet.

Qui vi mostreremo come fare il caso semplice.

4.1. Esaminare il file CSV fornito che elenca i percorsi dei file di input

Abbiamo già creato un file CSV che elenca i percorsi dei file di input, chiamato samplesheet.csv, che potete trovare nella directory 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

Questo file CSV include una riga di intestazione che nomina le colonne.

Avviso

I percorsi dei file nel CSV sono percorsi assoluti che devono corrispondere al vostro ambiente. Se non state eseguendo questo nell'ambiente di formazione che forniamo, dovrete aggiornare i percorsi per corrispondere al vostro sistema.

4.2. Aggiornare il parametro e il profilo di test

Cambiate il parametro input per puntare al file samplesheet.csv invece di elencare i singoli campioni.

Ripristinate l'annotazione di tipo nel blocco params (poiché è di nuovo un singolo percorso):

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

Poi aggiornate il profilo di test per puntare al file di testo:

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"
}

L'elenco dei file non vive più nel codice, il che è un grande passo nella giusta direzione.

4.3. Aggiornare la fabbrica di canali per analizzare l'input CSV

Attualmente, la nostra fabbrica di canali di input tratta qualsiasi file che le diamo come gli input di dati che vogliamo fornire al processo di indicizzazione. Poiché ora le stiamo dando un file che elenca i percorsi dei file di input, dobbiamo cambiare il suo comportamento per analizzare il file e trattare i percorsi dei file che contiene come gli input di dati.

Possiamo farlo usando lo stesso pattern che abbiamo usato nella Parte 2 di Hello Nextflow: applicando l'operatore splitCsv() per analizzare il file, poi un'operazione map per estrarre il percorso del file da ogni riga.

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 nuova rispetto a quanto avete incontrato nel corso Hello Nextflow è che questo CSV ha una riga di intestazione, quindi aggiungiamo header: true alla chiamata splitCsv(). Questo ci permette di fare riferimento alle colonne per nome nell'operazione map: row.reads_bam estrae il percorso del file dalla colonna reads_bam di ogni riga.

Suggerimento

Se non siete sicuri di capire cosa stanno facendo gli operatori qui, questa è un'altra grande opportunità per usare l'operatore .view() per vedere come appaiono i contenuti del canale prima e dopo averli applicati.

4.4. Eseguire il flusso di lavoro

Eseguite il flusso di lavoro un'ultima volta.

nextflow run genomics.nf -profile test
Output 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 ✔

Questo dovrebbe produrre lo stesso risultato di prima. Il nostro semplice flusso di lavoro di variant calling ora ha tutte le funzionalità di base che volevamo.

Takeaway

Sapete come creare un flusso di lavoro modulare multi-passaggio per indicizzare un file BAM e applicare il variant calling per campione usando GATK.

Più in generale, avete imparato come usare componenti e logica essenziali di Nextflow per costruire una semplice pipeline genomica che fa un lavoro reale, tenendo conto delle idiosincrasie dei formati di file genomici e dei requisiti degli strumenti.

Cosa c'è dopo?

Prendetevi una pausa! È stato molto.

Quando vi sentite riposati, passate alla Parte 3, dove imparerete come trasformare questo semplice flusso di lavoro di variant calling per campione per applicare il joint variant calling ai dati.