Vai al contenuto

Parte 3: Joint calling su una coorte

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

In precedenza, avete costruito una pipeline di variant calling per-campione che processava i dati di ciascun campione in modo indipendente. Ora la estenderemo per implementare il joint variant calling, come descritto nella Parte 1 (caso d'uso multi-campione).

Come iniziare da questa sezione

Questa sezione del corso presuppone che abbiate completato la Parte 1: Panoramica del metodo, la Parte 2: Variant calling per-campione e che abbiate una pipeline genomics.nf funzionante.

Se non avete completato la Parte 2 o volete ripartire da zero per questa parte, potete usare la soluzione della Parte 2 come punto di partenza. Eseguite questi comandi dall'interno della directory nf4-science/genomics/:

cp solutions/part2/genomics-2.nf genomics.nf
cp solutions/part2/nextflow.config .
cp solutions/part2/modules/* modules/

Questo vi fornisce un flusso di lavoro completo di variant calling per-campione. Potete verificare che funzioni correttamente eseguendo il seguente comando:

nextflow run genomics.nf -profile test

Compito

In questa parte del corso, estenderemo il flusso di lavoro per fare quanto segue:

  1. Generare un file indice per ciascun file BAM di input usando Samtools
  2. Eseguire GATK HaplotypeCaller su ciascun file BAM di input per generare un GVCF delle chiamate di varianti genomiche per-campione
  3. Raccogliere tutti i GVCF e combinarli in un data store GenomicsDB
  4. Eseguire il joint genotyping sul data store GVCF combinato per produrre un VCF a livello di coorte
BAMSamtools indexBAM indexIntervalsReference+ index & dictGATK HaplotypeCallerGVCF + indexx multiple samplesGVCF + indexGVCF + indexGVCF + indexGenomicsDBvariant storeGATK GenomicsDBImportGATK GenotypeGVCFsJoint-calledVCFGVCF mode

Questo implementa il metodo descritto nella Parte 1: Panoramica del metodo (seconda sezione che copre il caso d'uso multi-campione) e si basa direttamente sul flusso di lavoro prodotto dalla Parte 2: Variant calling per-campione.

Piano della lezione

Abbiamo suddiviso questo in due fasi:

  1. Modificare il passaggio di variant calling per-campione per produrre un GVCF. Questo copre l'aggiornamento dei comandi e degli output del processo.
  2. Aggiungere un passaggio di joint genotyping che combina e genotipizza i GVCF per-campione. Questo introduce l'operatore collect(), le closure Groovy per la costruzione della riga di comando e i processi multi-comando.

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

Suggerimento

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

1. Modificare il passaggio di variant calling per-campione per produrre un GVCF

La pipeline della Parte 2 produce file VCF, ma il joint calling richiede file GVCF. Dobbiamo attivare la modalità di variant calling GVCF e aggiornare l'estensione del file di output.

Ricordate il comando di variant calling GVCF dalla Parte 1:

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

Rispetto al comando base HaplotypeCaller che abbiamo incapsulato nella Parte 2, le differenze sono il parametro -ERC GVCF e l'estensione di output .g.vcf.

1.1. Dire a HaplotypeCaller di emettere un GVCF e aggiornare l'estensione di output

Aprite il file modulo modules/gatk_haplotypecaller.nf per apportare due modifiche:

  • Aggiungere il parametro -ERC GVCF al comando GATK HaplotypeCaller;
  • Aggiornare il percorso del file di output per usare l'estensione corrispondente .g.vcf, secondo la convenzione GATK.

Assicuratevi di aggiungere una barra rovesciata (\) alla fine della riga precedente quando aggiungete -ERC GVCF.

modules/gatk_haplotypecaller.nf
22
23
24
25
26
27
28
29
    """
    gatk HaplotypeCaller \
        -R ${ref_fasta} \
        -I ${input_bam} \
        -O ${input_bam}.g.vcf \
        -L ${interval_list} \
        -ERC GVCF
    """
modules/gatk_haplotypecaller.nf
22
23
24
25
26
27
28
    """
    gatk HaplotypeCaller \
        -R ${ref_fasta} \
        -I ${input_bam} \
        -O ${input_bam}.vcf \
        -L ${interval_list}
    """

Dobbiamo anche aggiornare il blocco output per corrispondere alla nuova estensione del file. Poiché abbiamo cambiato l'output del comando da .vcf a .g.vcf, il blocco output: del processo deve riflettere lo stesso cambiamento.

1.2. Aggiornare l'estensione del file di output nel blocco outputs del processo

modules/gatk_haplotypecaller.nf
17
18
19
    output:
    path "${input_bam}.g.vcf"     , emit: vcf
    path "${input_bam}.g.vcf.idx" , emit: idx
modules/gatk_haplotypecaller.nf
17
18
19
    output:
    path "${input_bam}.vcf"     , emit: vcf
    path "${input_bam}.vcf.idx" , emit: idx

Dobbiamo anche aggiornare la configurazione di publish e output del flusso di lavoro per riflettere i nuovi output GVCF.

1.3. Aggiornare i target di publish per i nuovi output GVCF

Poiché ora stiamo producendo GVCF invece di VCF, dovremmo aggiornare la sezione publish: del flusso di lavoro per usare nomi più descrittivi. Organizzeremo anche i file GVCF nella loro sottodirectory per maggiore chiarezza.

genomics.nf
47
48
49
50
    publish:
    indexed_bam = SAMTOOLS_INDEX.out
    gvcf = GATK_HAPLOTYPECALLER.out.vcf
    gvcf_idx = GATK_HAPLOTYPECALLER.out.idx
genomics.nf
47
48
49
50
    publish:
    indexed_bam = SAMTOOLS_INDEX.out
    vcf = GATK_HAPLOTYPECALLER.out.vcf
    vcf_idx = GATK_HAPLOTYPECALLER.out.idx

Ora aggiorniamo il blocco output per corrispondere.

1.4. Aggiornare il blocco output per la nuova struttura di directory

Dobbiamo anche aggiornare il blocco output per mettere i file GVCF in una sottodirectory gvcf.

genomics.nf
53
54
55
56
57
58
59
60
61
62
63
output {
    indexed_bam {
        path 'indexed_bam'
    }
    gvcf {
        path 'gvcf'
    }
    gvcf_idx {
        path 'gvcf'
    }
}
genomics.nf
53
54
55
56
57
58
59
60
61
62
63
output {
    indexed_bam {
        path 'bam'
    }
    vcf {
        path 'vcf'
    }
    vcf_idx {
        path 'vcf'
    }
}

Con il modulo, i target di publish e il blocco output tutti aggiornati, possiamo testare le modifiche.

1.5. Eseguire la pipeline

Eseguite il flusso di lavoro per verificare che le modifiche funzionino.

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

┃ Launching `genomics.nf` [nostalgic_franklin] DSL2 - revision: f2c0a93c6a

executor >  local (6)
[cc/fbc705] SAMTOOLS_INDEX (3)       | 3 of 3 ✔
[27/0d7eb9] GATK_HAPLOTYPECALLER (2) | 3 of 3 ✔

L'output di Nextflow appare uguale a prima, ma i file .g.vcf e i loro file indice sono ora organizzati in sottodirectory.

Contenuto della directory (symlink abbreviati) 
results/
├── gvcf/
│   ├── reads_father.bam.g.vcf -> */27/0d7eb9*/reads_father.bam.g.vcf
│   ├── reads_father.bam.g.vcf.idx -> */27/0d7eb9*/reads_father.bam.g.vcf.idx
│   ├── reads_mother.bam.g.vcf -> */e4/4ed55e*/reads_mother.bam.g.vcf
│   ├── reads_mother.bam.g.vcf.idx -> */e4/4ed55e*/reads_mother.bam.g.vcf.idx
│   ├── reads_son.bam.g.vcf -> */08/e95962*/reads_son.bam.g.vcf
│   └── reads_son.bam.g.vcf.idx -> */08/e95962*/reads_son.bam.g.vcf.idx
└── indexed_bam/
    ├── reads_father.bam -> */9a/c7a873*/reads_father.bam
    ├── reads_father.bam.bai -> */9a/c7a873*/reads_father.bam.bai
    ├── reads_mother.bam -> */f1/8d8486*/reads_mother.bam
    ├── reads_mother.bam.bai -> */f1/8d8486*/reads_mother.bam.bai
    ├── reads_son.bam -> */cc/fbc705*/reads_son.bam
    └── reads_son.bam.bai -> */cc/fbc705*/reads_son.bam.bai

Se aprite uno dei file GVCF e scorrete attraverso di esso, potete verificare che GATK HaplotypeCaller abbia prodotto file GVCF come richiesto.

Takeaway

Quando modificate il nome del file di output di un comando di uno strumento, il blocco output: del processo e la configurazione publish/output devono essere aggiornati per corrispondere.

Cosa c'è dopo?

Imparate a raccogliere i contenuti di un canale e passarli al processo successivo come singolo input.

2. Aggiungere un passaggio di joint genotyping

Ora dobbiamo raccogliere i GVCF per-campione, combinarli in un data store GenomicsDB ed eseguire il joint genotyping per produrre un VCF a livello di coorte. Come trattato nella Parte 1, questa è un'operazione a due strumenti: GenomicsDBImport combina i GVCF, poi GenotypeGVCFs produce le chiamate di varianti finali. Incapsuleremo entrambi gli strumenti in un singolo processo chiamato GATK_JOINTGENOTYPING.

Ricordate i due comandi dalla Parte 1:

gatk GenomicsDBImport \
    -V reads_mother.g.vcf \
    -V reads_father.g.vcf \
    -V reads_son.g.vcf \
    -L /data/ref/intervals.bed \
    --genomicsdb-workspace-path family_trio_gdb

gatk GenotypeGVCFs \
    -R /data/ref/ref.fasta \
    -V gendb://family_trio_gdb \
    -O family_trio.vcf

Il primo comando prende i GVCF per-campione e un file di intervalli, e produce un data store GenomicsDB. Il secondo prende quel data store, un genoma di riferimento, e produce il VCF finale a livello di coorte. L'URI del container è lo stesso di HaplotypeCaller: community.wave.seqera.io/library/gatk4:4.5.0.0--730ee8817e436867.

2.1. Configurare gli input

Il processo di joint genotyping necessita di due tipi di input che non abbiamo ancora: un nome arbitrario di coorte e gli output GVCF raccolti da tutti i campioni raggruppati insieme.

2.1.1. Aggiungere un parametro cohort_name

Dobbiamo fornire un nome arbitrario per la coorte. Più avanti nella serie di formazione imparerete come usare i metadati dei campioni per questo tipo di cose, ma per ora dichiariamo semplicemente un parametro CLI usando params.

genomics.nf
18
19
20
21
22
    intervals: Path

    // Nome base per il file di output finale
    cohort_name: String
}
genomics.nf
18
19
    intervals: Path
}

Useremo questo parametro per nominare il file di output finale.

2.1.2. Aggiungere un valore predefinito per cohort_name nel profilo test

Aggiungiamo anche un valore predefinito per il parametro cohort_name nel profilo test:

nextflow.config
 4
 5
 6
 7
 8
 9
10
11
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"
    params.cohort_name = "family_trio"
}
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"
}

Successivamente, dovremo raccogliere gli output per-campione in modo che possano essere processati insieme.

2.1.3. Raccogliere gli output di HaplotypeCaller attraverso i campioni

Se dovessimo collegare il canale di output da GATK_HAPLOTYPECALLER direttamente al nuovo processo, Nextflow chiamerebbe il processo su ciascun GVCF di campione separatamente. Vogliamo raggruppare tutti e tre i GVCF (e i loro file indice) in modo che Nextflow li passi tutti insieme a una singola chiamata di processo.

Possiamo farlo usando l'operatore di canale collect(). Aggiungete le seguenti righe al corpo del workflow, subito dopo la chiamata a GATK_HAPLOTYPECALLER:

genomics.nf
48
49
50
51
52
53
        intervals_file
    )

    // Raccoglie gli output di variant calling attraverso i campioni
    all_gvcfs_ch = GATK_HAPLOTYPECALLER.out.vcf.collect()
    all_idxs_ch = GATK_HAPLOTYPECALLER.out.idx.collect()
genomics.nf
48
49
        intervals_file
    )

Scomponendo questo:

  1. Prendiamo il canale di output da GATK_HAPLOTYPECALLER usando la proprietà .out.
  2. Poiché abbiamo nominato gli output usando emit: nella sezione 1, possiamo selezionare i GVCF con .vcf e i file indice con .idx. Senza output nominati, dovremmo usare .out[0] e .out[1].
  3. L'operatore collect() raggruppa tutti i file in un singolo elemento, quindi all_gvcfs_ch contiene tutti e tre i GVCF insieme, e all_idxs_ch contiene tutti e tre i file indice insieme.

Possiamo raccogliere i GVCF e i loro file indice separatamente (invece di tenerli insieme in tuple) perché Nextflow metterà in stage tutti i file di input insieme per l'esecuzione, quindi i file indice saranno presenti accanto ai GVCF.

Suggerimento

Potete usare l'operatore view() per ispezionare i contenuti dei canali prima e dopo l'applicazione degli operatori di canale.

2.2. Scrivere il processo di joint genotyping e chiamarlo nel flusso di lavoro

Seguendo lo stesso schema che abbiamo usato nella Parte 2, scriveremo la definizione del processo in un file modulo, lo importeremo nel flusso di lavoro e lo chiameremo sugli input che abbiamo appena preparato.

2.2.1. Costruire una stringa per dare a ciascun GVCF un argomento -V

Prima di iniziare a riempire la definizione del processo, c'è una cosa da risolvere. Il comando GenomicsDBImport si aspetta un argomento -V separato per ciascun file GVCF, così:

gatk GenomicsDBImport \
    -V reads_mother.bam.g.vcf \
    -V reads_father.bam.g.vcf \
    -V reads_son.bam.g.vcf \
    ...

Se dovessimo scrivere -V ${all_gvcfs_ch}, Nextflow concatenerebbe semplicemente i nomi dei file e quella parte del comando apparirebbe così:

-V reads_mother.bam.g.vcf reads_father.bam.g.vcf reads_son.bam.g.vcf

Ma abbiamo bisogno che la stringa appaia così:

-V reads_mother.bam.g.vcf -V reads_father.bam.g.vcf -V reads_son.bam.g.vcf

Importante, dobbiamo costruire questa stringa dinamicamente da qualunque file sia nel canale raccolto. Nextflow (tramite Groovy) fornisce un modo conciso per farlo:

def gvcfs_line = all_gvcfs.collect { gvcf -> "-V ${gvcf}" }.join(' ')

Scomponendo questo:

  1. all_gvcfs.collect { gvcf -> "-V ${gvcf}" } itera su ciascun percorso di file e antepone -V ad esso, producendo ["-V A.g.vcf", "-V B.g.vcf", "-V C.g.vcf"].
  2. .join(' ') li concatena con spazi: "-V A.g.vcf -V B.g.vcf -V C.g.vcf".
  3. Il risultato è assegnato a una variabile locale gvcfs_line (definita con def), che possiamo interpolare nel template del comando.

Questa riga va all'interno del blocco script: del processo, prima del template del comando. Potete inserire codice Groovy arbitrario tra script: e l'apertura """ del template del comando.

Poi sarete in grado di riferirvi a quella intera stringa come gvcfs_line nel blocco script: del processo.

2.2.2. Riempire il modulo per il processo di joint genotyping

Successivamente, possiamo affrontare la scrittura del processo completo.

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

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

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

/*
 * Combina i GVCF in un datastore GenomicsDB ed esegue il joint genotyping per produrre chiamate a livello di coorte
 */
process GATK_JOINTGENOTYPING {

    container

    input:

    output:

    script:
    """

    """
}
modules/gatk_jointgenotyping.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
31
32
33
34
35
36
37
#!/usr/bin/env nextflow

/*
 * Combina i GVCF in un datastore GenomicsDB ed esegue il joint genotyping per produrre chiamate a livello di coorte
 */
process GATK_JOINTGENOTYPING {

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

    input:
    path all_gvcfs
    path all_idxs
    path interval_list
    val cohort_name
    path ref_fasta
    path ref_index
    path ref_dict

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

    script:
    def gvcfs_line = all_gvcfs.collect { gvcf -> "-V ${gvcf}" }.join(' ')
    """
    gatk GenomicsDBImport \
        ${gvcfs_line} \
        -L ${interval_list} \
        --genomicsdb-workspace-path ${cohort_name}_gdb

    gatk GenotypeGVCFs \
        -R ${ref_fasta} \
        -V gendb://${cohort_name}_gdb \
        -L ${interval_list} \
        -O ${cohort_name}.joint.vcf
    """
}

Ci sono diverse cose che vale la pena evidenziare qui.

Come in precedenza, diversi input sono elencati anche se i comandi non vi fanno riferimento direttamente: all_idxs, ref_index e ref_dict. Elencarli assicura che Nextflow metta in stage questi file nella directory di lavoro accanto ai file che appaiono nei comandi, che GATK si aspetta di trovare in base alle convenzioni di denominazione.

La variabile gvcfs_line usa la closure Groovy descritta sopra per costruire gli argomenti -V per GenomicsDBImport.

Questo processo esegue due comandi in serie, proprio come fareste nel terminale. GenomicsDBImport combina i GVCF per-campione in un data store, poi GenotypeGVCFs legge quel data store e produce il VCF finale a livello di coorte. Il data store GenomicsDB (${cohort_name}_gdb) è un artefatto intermedio usato solo all'interno del processo; non appare nel blocco output.

Una volta completato questo, il processo è pronto per l'uso. Per usarlo nel flusso di lavoro, dovrete importare il modulo e aggiungere una chiamata di processo.

2.2.3. Importare il modulo

Aggiungete l'istruzione import a genomics.nf, sotto le istruzioni import esistenti:

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

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

2.2.4. Aggiungere la chiamata di processo

Aggiungete la chiamata a GATK_JOINTGENOTYPING nel corpo del flusso di lavoro, dopo le righe collect():

genomics.nf
53
54
55
56
57
58
59
60
61
62
63
64
    all_idxs_ch = GATK_HAPLOTYPECALLER.out.idx.collect()

    // Combina i GVCF in un data store GenomicsDB e applica il joint genotyping
    GATK_JOINTGENOTYPING(
        all_gvcfs_ch,
        all_idxs_ch,
        intervals_file,
        params.cohort_name,
        ref_file,
        ref_index_file,
        ref_dict_file
    )
genomics.nf
53
    all_idxs_ch = GATK_HAPLOTYPECALLER.out.idx.collect()

Il processo è ora completamente collegato. Successivamente, configuriamo come vengono pubblicati gli output.

2.3. Configurare la gestione degli output

Dobbiamo pubblicare gli output del VCF joint. Aggiungete target di publish e voci del blocco output per i risultati del joint genotyping.

2.3.1. Aggiungere target di publish per il VCF joint

Aggiungete il VCF joint e il suo indice alla sezione publish: del flusso di lavoro:

genomics.nf
66
67
68
69
70
71
    publish:
    indexed_bam = SAMTOOLS_INDEX.out
    gvcf = GATK_HAPLOTYPECALLER.out.vcf
    gvcf_idx = GATK_HAPLOTYPECALLER.out.idx
    joint_vcf = GATK_JOINTGENOTYPING.out.vcf
    joint_vcf_idx = GATK_JOINTGENOTYPING.out.idx
genomics.nf
66
67
68
69
    publish:
    indexed_bam = SAMTOOLS_INDEX.out
    gvcf = GATK_HAPLOTYPECALLER.out.vcf
    gvcf_idx = GATK_HAPLOTYPECALLER.out.idx

Ora aggiornate il blocco output per corrispondere.

2.3.2. Aggiungere voci del blocco output per il VCF joint

Aggiungete voci per i file VCF joint. Li metteremo alla radice della directory results poiché questo è l'output finale.

genomics.nf
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
output {
    indexed_bam {
        path 'indexed_bam'
    }
    gvcf {
        path 'gvcf'
    }
    gvcf_idx {
        path 'gvcf'
    }
    joint_vcf {
        path '.'
    }
    joint_vcf_idx {
        path '.'
    }
}
genomics.nf
74
75
76
77
78
79
80
81
82
83
84
output {
    indexed_bam {
        path 'indexed_bam'
    }
    gvcf {
        path 'gvcf'
    }
    gvcf_idx {
        path 'gvcf'
    }
}

Con il processo, i target di publish e il blocco output tutti a posto, possiamo testare il flusso di lavoro completo.

2.4. Eseguire il flusso di lavoro

Eseguite il flusso di lavoro per verificare che tutto funzioni.

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` [crazy_marconi] DSL2 - revision: 5da9afc841

executor >  local (1)
[9a/c7a873] SAMTOOLS_INDEX (2)       | 3 of 3, cached: 3 ✔
[e4/4ed55e] GATK_HAPLOTYPECALLER (2) | 3 of 3, cached: 3 ✔
[a6/7cc8ed] GATK_JOINTGENOTYPING     | 1 of 1 ✔

I primi due passaggi sono in cache dall'esecuzione precedente, e il nuovo passaggio GATK_JOINTGENOTYPING viene eseguito una volta sugli input raccolti da tutti e tre i campioni. Il file di output finale, family_trio.joint.vcf (e il suo indice), sono nella directory results.

Contenuto della directory (symlink abbreviati) 
results/
├── family_trio.joint.vcf -> */a6/7cc8ed*/family_trio.joint.vcf
├── family_trio.joint.vcf.idx -> */a6/7cc8ed*/family_trio.joint.vcf.idx
├── gvcf/
│   ├── reads_father.bam.g.vcf -> */27/0d7eb9*/reads_father.bam.g.vcf
│   ├── reads_father.bam.g.vcf.idx -> */27/0d7eb9*/reads_father.bam.g.vcf.idx
│   ├── reads_mother.bam.g.vcf -> */e4/4ed55e*/reads_mother.bam.g.vcf
│   ├── reads_mother.bam.g.vcf.idx -> */e4/4ed55e*/reads_mother.bam.g.vcf.idx
│   ├── reads_son.bam.g.vcf -> */08/e95962*/reads_son.bam.g.vcf
│   └── reads_son.bam.g.vcf.idx -> */08/e95962*/reads_son.bam.g.vcf.idx
└── indexed_bam/
    ├── reads_father.bam -> */9a/c7a873*/reads_father.bam
    ├── reads_father.bam.bai -> */9a/c7a873*/reads_father.bam.bai
    ├── reads_mother.bam -> */f1/8d8486*/reads_mother.bam
    ├── reads_mother.bam.bai -> */f1/8d8486*/reads_mother.bam.bai
    ├── reads_son.bam -> */cc/fbc705*/reads_son.bam
    └── reads_son.bam.bai -> */cc/fbc705*/reads_son.bam.bai

Se aprite il file VCF joint, potete verificare che il flusso di lavoro abbia prodotto le chiamate di varianti attese.

family_trio.joint.vcf
40
41
42
43
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	reads_father	reads_mother	reads_son
20_10037292_10066351	3480	.	C	CT	1625.89	.	AC=5;AF=0.833;AN=6;BaseQRankSum=0.220;DP=85;ExcessHet=0.0000;FS=2.476;MLEAC=5;MLEAF=0.833;MQ=60.00;MQRankSum=0.00;QD=21.68;ReadPosRankSum=-1.147e+00;SOR=0.487	GT:AD:DP:GQ:PL	0/1:15,16:31:99:367,0,375	1/1:0,18:18:54:517,54,0	1/1:0,26:26:78:756,78,0
20_10037292_10066351	3520	.	AT	A	1678.89	.	AC=5;AF=0.833;AN=6;BaseQRankSum=1.03;DP=80;ExcessHet=0.0000;FS=2.290;MLEAC=5;MLEAF=0.833;MQ=60.00;MQRankSum=0.00;QD=22.39;ReadPosRankSum=0.701;SOR=0.730	GT:AD:DP:GQ:PL	0/1:18,13:31:99:296,0,424	1/1:0,18:18:54:623,54,0	1/1:0,26:26:78:774,78,0
20_10037292_10066351	3529	.	T	A	154.29	.	AC=1;AF=0.167;AN=6;BaseQRankSum=-5.440e-01;DP=104;ExcessHet=0.0000;FS=1.871;MLEAC=1;MLEAF=0.167;MQ=60.00;MQRankSum=0.00;QD=7.71;ReadPosRankSum=-1.158e+00;SOR=1.034	GT:AD:DP:GQ:PL	0/0:44,0:44:99:0,112,1347	0/1:12,8:20:99:163,0,328	0/0:39,0:39:99:0,105,1194

Ora avete un flusso di lavoro di joint variant calling automatizzato e completamente riproducibile!

Nota

Tenete presente che i file di dati che vi abbiamo fornito coprono solo una piccola porzione del cromosoma 20. La dimensione reale di un callset di varianti sarebbe contata in milioni di varianti. Ecco perché usiamo solo piccoli sottoinsiemi di dati per scopi di formazione!

Takeaway

Sapete come raccogliere output da un canale e raggrupparli come singolo input per un altro processo. Sapete anche come costruire una riga di comando usando le closure Groovy, e come eseguire comandi multipli in un singolo processo.

Cosa c'è dopo?

Datevi una grande pacca sulla spalla! Avete completato il corso Nextflow for Genomics.

Passate al riepilogo finale del corso per rivedere ciò che avete imparato e scoprire cosa viene dopo.