Przejdź do treści

Część 3: Wspólne wywoływanie wariantów dla kohorty

Tłumaczenie wspomagane przez AI - dowiedz się więcej i zasugeruj ulepszenia

Wcześniej zbudowałeś pipeline wywoływania wariantów dla pojedynczych próbek, który przetwarzał dane każdej próbki niezależnie. Teraz rozszerzymy go, aby zaimplementować wspólne wywoływanie wariantów, jak opisano w Części 1 (przypadek użycia dla wielu próbek).

Jak zacząć od tej sekcji

Ta sekcja kursu zakłada, że ukończyłeś Część 1: Przegląd metody, Część 2: Wywoływanie wariantów dla pojedynczych próbek i masz działający pipeline genomics.nf.

Jeśli nie ukończyłeś Części 2 lub chcesz zacząć od nowa w tej części, możesz użyć rozwiązania z Części 2 jako punktu wyjścia. Uruchom te polecenia z katalogu nf4-science/genomics/:

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

To da Ci kompletny workflow wywoływania wariantów dla pojedynczych próbek. Możesz sprawdzić, czy działa poprawnie, uruchamiając następujące polecenie:

nextflow run genomics.nf -profile test

Zadanie

W tej części kursu rozszerzymy workflow'a, aby wykonywał następujące operacje:

  1. Wygenerować plik indeksu dla każdego pliku BAM wejściowego przy użyciu Samtools
  2. Uruchomić GATK HaplotypeCaller na każdym pliku BAM wejściowym, aby wygenerować GVCF z wywołaniami wariantów genomowych dla pojedynczej próbki
  3. Zebrać wszystkie pliki GVCF i połączyć je w magazyn danych GenomicsDB
  4. Uruchomić wspólne genotypowanie na połączonym magazynie danych GVCF, aby wygenerować plik VCF na poziomie kohorty
BAMSamtools indexBAM indexIntervalsReference+ index & dictGATK HaplotypeCallerGVCF + indexx multiple samplesGVCF + indexGVCF + indexGVCF + indexGenomicsDBvariant storeGATK GenomicsDBImportGATK GenotypeGVCFsJoint-calledVCFGVCF mode

Implementuje to metodę opisaną w Części 1: Przegląd metody (druga sekcja dotycząca przypadku użycia dla wielu próbek) i buduje bezpośrednio na workflow'ie stworzonym w Części 2: Wywoływanie wariantów dla pojedynczych próbek.

Plan lekcji

Podzieliliśmy to na dwa etapy:

  1. Zmodyfikować krok wywoływania wariantów dla pojedynczych próbek, aby generował GVCF. Obejmuje to aktualizację poleceń i wyjść procesu.
  2. Dodać krok wspólnego genotypowania, który łączy i genotypuje pliki GVCF dla pojedynczych próbek. Wprowadza to operator collect(), domknięcia Groovy do konstruowania linii poleceń oraz procesy z wieloma poleceniami.

Automatyzuje to kroki z drugiej sekcji Części 1: Przegląd metody, gdzie uruchamiałeś te polecenia ręcznie w ich kontenerach.

Wskazówka

Upewnij się, że jesteś we właściwym katalogu roboczym: cd /workspaces/training/nf4-science/genomics

1. Zmodyfikuj krok wywoływania wariantów dla pojedynczych próbek, aby generował GVCF

Pipeline z Części 2 generuje pliki VCF, ale wspólne wywoływanie wymaga plików GVCF. Musimy włączyć tryb wywoływania wariantów GVCF i zaktualizować rozszerzenie pliku wyjściowego.

Przypomnij sobie polecenie wywoływania wariantów GVCF z Części 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

W porównaniu z podstawowym poleceniem HaplotypeCaller, które opakowano w Części 2, różnice to parametr -ERC GVCF i rozszerzenie wyjściowe .g.vcf.

1.1. Powiedz HaplotypeCaller, aby emitował GVCF i zaktualizuj rozszerzenie wyjściowe

Otwórz plik modułu modules/gatk_haplotypecaller.nf, aby wprowadzić dwie zmiany:

  • Dodaj parametr -ERC GVCF do polecenia GATK HaplotypeCaller;
  • Zaktualizuj ścieżkę pliku wyjściowego, aby używała odpowiadającego rozszerzenia .g.vcf, zgodnie z konwencją GATK.

Upewnij się, że dodajesz ukośnik odwrotny (\) na końcu poprzedniej linii, gdy dodajesz -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}
    """

Musimy również zaktualizować blok wyjściowy, aby pasował do nowego rozszerzenia pliku. Ponieważ zmieniliśmy wyjście polecenia z .vcf na .g.vcf, blok output: procesu musi odzwierciedlać tę samą zmianę.

1.2. Zaktualizuj rozszerzenie pliku wyjściowego w bloku wyjść procesu

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

Musimy również zaktualizować konfigurację publikowania i wyjścia workflow'a, aby odzwierciedlała nowe wyjścia GVCF.

1.3. Zaktualizuj cele publikowania dla nowych wyjść GVCF

Ponieważ teraz generujemy pliki GVCF zamiast VCF, powinniśmy zaktualizować sekcję publish: workflow'a, aby używała bardziej opisowych nazw. Zorganizujemy również pliki GVCF w ich własnym podkatalogu dla przejrzystości.

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

Teraz zaktualizuj blok wyjściowy, aby pasował.

1.4. Zaktualizuj blok wyjściowy dla nowej struktury katalogów

Musimy również zaktualizować blok output, aby umieścić pliki GVCF w podkatalogu 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'
    }
}

Po zaktualizowaniu modułu, celów publikowania i bloku wyjściowego możemy przetestować zmiany.

1.5. Uruchom pipeline'a

Uruchom workflow'a, aby zweryfikować, że zmiany działają.

nextflow run genomics.nf -profile test
Wyjście polecenia 
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 ✔

Wyjście Nextflow'a wygląda tak samo jak wcześniej, ale pliki .g.vcf i ich pliki indeksów są teraz zorganizowane w podkatalogach.

Zawartość katalogu (skrócone dowiązania symboliczne) 
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

Jeśli otworzysz jeden z plików GVCF i przewiniesz go, możesz zweryfikować, że GATK HaplotypeCaller wygenerował pliki GVCF zgodnie z żądaniem.

Podsumowanie

Gdy zmieniasz nazwę pliku wyjściowego polecenia narzędzia, blok output: procesu oraz konfiguracja publikowania/wyjścia muszą zostać zaktualizowane, aby pasowały.

Co dalej?

Naucz się zbierać zawartość kanału i przekazywać ją do następnego procesu jako pojedyncze wejście.

2. Dodaj krok wspólnego genotypowania

Teraz musimy zebrać pliki GVCF dla pojedynczych próbek, połączyć je w magazyn danych GenomicsDB i uruchomić wspólne genotypowanie, aby wygenerować plik VCF na poziomie kohorty. Jak omówiono w Części 1, jest to operacja dwóch narzędzi: GenomicsDBImport łączy pliki GVCF, a następnie GenotypeGVCFs generuje końcowe wywołania wariantów. Opakujemy oba narzędzia w jeden proces o nazwie GATK_JOINTGENOTYPING.

Przypomnij sobie dwa polecenia z Części 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

Pierwsze polecenie przyjmuje pliki GVCF dla pojedynczych próbek i plik interwałów, a następnie generuje magazyn danych GenomicsDB. Drugie przyjmuje ten magazyn danych, genom referencyjny i generuje końcowy plik VCF na poziomie kohorty. URI kontenera jest taki sam jak dla HaplotypeCaller: community.wave.seqera.io/library/gatk4:4.5.0.0--730ee8817e436867.

2.1. Skonfiguruj wejścia

Proces wspólnego genotypowania potrzebuje dwóch rodzajów wejść, których jeszcze nie mamy: dowolnej nazwy kohorty oraz zebranych wyjść GVCF ze wszystkich próbek połączonych razem.

2.1.1. Dodaj parametr cohort_name

Musimy podać dowolną nazwę dla kohorty. Później w serii szkoleń nauczysz się, jak używać metadanych próbek do tego rodzaju rzeczy, ale na razie po prostu deklarujemy parametr CLI używając params.

genomics.nf
18
19
20
21
22
    intervals: Path

    // Nazwa bazowa dla końcowego pliku wyjściowego
    cohort_name: String
}
genomics.nf
18
19
    intervals: Path
}

Użyjemy tego parametru do nazwania końcowego pliku wyjściowego.

2.1.2. Dodaj wartość domyślną dla cohort_name w profilu testowym

Dodajemy również wartość domyślną dla parametru cohort_name w profilu testowym:

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

Następnie będziemy musieli zebrać wyjścia dla pojedynczych próbek, aby mogły być przetwarzane razem.

2.1.3. Zbierz wyjścia HaplotypeCaller z wszystkich próbek

Gdybyśmy podłączyli kanał wyjściowy z GATK_HAPLOTYPECALLER bezpośrednio do nowego procesu, Nextflow wywołałby proces na każdym pliku GVCF próbki osobno. Chcemy połączyć wszystkie trzy pliki GVCF (i ich pliki indeksów), aby Nextflow przekazał je wszystkie razem do pojedynczego wywołania procesu.

Możemy to zrobić używając operatora kanału collect(). Dodaj następujące linie do ciała workflow, zaraz po wywołaniu GATK_HAPLOTYPECALLER:

genomics.nf
48
49
50
51
52
53
        intervals_file
    )

    // Zbierz wyjścia wywoływania wariantów z wszystkich próbek
    all_gvcfs_ch = GATK_HAPLOTYPECALLER.out.vcf.collect()
    all_idxs_ch = GATK_HAPLOTYPECALLER.out.idx.collect()
genomics.nf
48
49
        intervals_file
    )

Rozbijając to na części:

  1. Pobieramy kanał wyjściowy z GATK_HAPLOTYPECALLER używając właściwości .out.
  2. Ponieważ nazwaliśmy wyjścia używając emit: w sekcji 1, możemy wybrać pliki GVCF za pomocą .vcf, a pliki indeksów za pomocą .idx. Bez nazwanych wyjść musielibyśmy użyć .out[0] i .out[1].
  3. Operator collect() łączy wszystkie pliki w jeden element, więc all_gvcfs_ch zawiera wszystkie trzy pliki GVCF razem, a all_idxs_ch zawiera wszystkie trzy pliki indeksów razem.

Możemy zbierać pliki GVCF i ich pliki indeksów osobno (w przeciwieństwie do trzymania ich razem w krotkach), ponieważ Nextflow przygotuje wszystkie pliki wejściowe razem do wykonania, więc pliki indeksów będą obecne obok plików GVCF.

Wskazówka

Możesz użyć operatora view(), aby sprawdzić zawartość kanałów przed i po zastosowaniu operatorów kanału.

2.2. Napisz proces wspólnego genotypowania i wywołaj go w workflow'ie

Postępując zgodnie z tym samym wzorcem, którego użyliśmy w Części 2, napiszemy definicję procesu w pliku modułu, zaimportujemy ją do workflow'a i wywołamy na wejściach, które właśnie przygotowaliśmy.

2.2.1. Skonstruuj ciąg znaków, aby nadać każdemu GVCF argument -V

Zanim zaczniemy wypełniać definicję procesu, jest jedna rzecz do rozwiązania. Polecenie GenomicsDBImport oczekuje osobnego argumentu -V dla każdego pliku GVCF, w ten sposób:

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

Gdybyśmy napisali -V ${all_gvcfs_ch}, Nextflow po prostu połączyłby nazwy plików i ta część polecenia wyglądałaby tak:

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

Ale potrzebujemy, aby ciąg wyglądał tak:

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

Co ważne, musimy skonstruować ten ciąg dynamicznie z dowolnych plików znajdujących się w zebranym kanale. Nextflow (poprzez Groovy) zapewnia zwięzły sposób, aby to zrobić:

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

Rozbijając to na części:

  1. all_gvcfs.collect { gvcf -> "-V ${gvcf}" } iteruje po każdej ścieżce pliku i dodaje przed nią -V, generując ["-V A.g.vcf", "-V B.g.vcf", "-V C.g.vcf"].
  2. .join(' ') łączy je spacjami: "-V A.g.vcf -V B.g.vcf -V C.g.vcf".
  3. Wynik jest przypisywany do zmiennej lokalnej gvcfs_line (zdefiniowanej za pomocą def), którą możemy interpolować w szablonie polecenia.

Ta linia trafia do bloku script: procesu, przed szablonem polecenia. Możesz umieścić dowolny kod Groovy między script: a otwierającym """ szablonu polecenia.

Następnie będziesz mógł odwoływać się do całego tego ciągu jako gvcfs_line w bloku script: procesu.

2.2.2. Wypełnij moduł dla procesu wspólnego genotypowania

Następnie możemy zająć się napisaniem pełnego procesu.

Otwórz modules/gatk_jointgenotyping.nf i przeanalizuj zarys definicji procesu.

Wypełnij definicję procesu używając informacji podanych powyżej, a następnie sprawdź swoją pracę z rozwiązaniem w zakładce "Po" poniżej.

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

/*
 * Połącz pliki GVCF w magazyn danych GenomicsDB i uruchom wspólne genotypowanie, aby wygenerować wywołania na poziomie kohorty
 */
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

/*
 * Połącz pliki GVCF w magazyn danych GenomicsDB i uruchom wspólne genotypowanie, aby wygenerować wywołania na poziomie kohorty
 */
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
    """
}

Jest kilka rzeczy wartych uwagi.

Jak wcześniej, kilka wejść jest wymienionych, mimo że polecenia nie odwołują się do nich bezpośrednio: all_idxs, ref_index i ref_dict. Wymienienie ich zapewnia, że Nextflow przygotuje te pliki w katalogu roboczym obok plików, które pojawiają się w poleceniach, czego GATK oczekuje znaleźć na podstawie konwencji nazewnictwa.

Zmienna gvcfs_line używa domknięcia Groovy opisanego powyżej do skonstruowania argumentów -V dla GenomicsDBImport.

Ten proces uruchamia dwa polecenia szeregowo, tak jak zrobiłbyś to w terminalu. GenomicsDBImport łączy pliki GVCF dla pojedynczych próbek w magazyn danych, a następnie GenotypeGVCFs odczytuje ten magazyn danych i generuje końcowy plik VCF na poziomie kohorty. Magazyn danych GenomicsDB (${cohort_name}_gdb) jest artefaktem pośrednim używanym tylko w ramach procesu; nie pojawia się w bloku wyjściowym.

Po ukończeniu tego proces jest gotowy do użycia. Aby użyć go w workflow'ie, musisz zaimportować moduł i dodać wywołanie procesu.

2.2.3. Zaimportuj moduł

Dodaj instrukcję importu do genomics.nf, poniżej istniejących instrukcji importu:

genomics.nf
3
4
5
6
// Instrukcje INCLUDE modułów
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
// Instrukcje INCLUDE modułów
include { SAMTOOLS_INDEX } from './modules/samtools_index.nf'
include { GATK_HAPLOTYPECALLER } from './modules/gatk_haplotypecaller.nf'

Proces jest teraz dostępny w zakresie workflow'a.

2.2.4. Dodaj wywołanie procesu

Dodaj wywołanie GATK_JOINTGENOTYPING w ciele workflow'a, po liniach collect():

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

    // Połącz pliki GVCF w magazyn danych GenomicsDB i zastosuj wspólne genotypowanie
    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()

Proces jest teraz w pełni podłączony. Następnie skonfigurujemy sposób publikowania wyjść.

2.3. Skonfiguruj obsługę wyjść

Musimy opublikować wyjścia wspólnego VCF. Dodaj cele publikowania i wpisy bloku wyjściowego dla wyników wspólnego genotypowania.

2.3.1. Dodaj cele publikowania dla wspólnego VCF

Dodaj wspólny VCF i jego indeks do sekcji publish: workflow'a:

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

Teraz zaktualizuj blok wyjściowy, aby pasował.

2.3.2. Dodaj wpisy bloku wyjściowego dla wspólnego VCF

Dodaj wpisy dla plików wspólnego VCF. Umieścimy je w katalogu głównym katalogu wyników, ponieważ jest to końcowe wyjście.

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

Mając proces, cele publikowania i blok wyjściowy na miejscu, możemy przetestować kompletny workflow.

2.4. Uruchom workflow'a

Uruchom workflow'a, aby zweryfikować, że wszystko działa.

nextflow run genomics.nf -profile test -resume
Wyjście polecenia 
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 ✔

Pierwsze dwa kroki są zbuforowane z poprzedniego uruchomienia, a nowy krok GATK_JOINTGENOTYPING uruchamia się raz na zebranych wejściach ze wszystkich trzech próbek. Końcowy plik wyjściowy, family_trio.joint.vcf (i jego indeks), znajduje się w katalogu wyników.

Zawartość katalogu (skrócone dowiązania symboliczne) 
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

Jeśli otworzysz wspólny plik VCF, możesz zweryfikować, że workflow wygenerował oczekiwane wywołania wariantów.

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

Masz teraz zautomatyzowany, w pełni odtwarzalny workflow wspólnego wywoływania wariantów!

Uwaga

Pamiętaj, że pliki danych, które Ci przekazaliśmy, obejmują tylko niewielką część chromosomu 20. Rzeczywisty rozmiar zestawu wywołań wariantów byłby liczony w milionach wariantów. Dlatego używamy tylko małych podzbiorów danych do celów szkoleniowych!

Podsumowanie

Wiesz, jak zbierać wyjścia z kanału i łączyć je jako pojedyncze wejście do innego procesu. Wiesz również, jak konstruować linię poleceń używając domknięć Groovy oraz jak uruchamiać wiele poleceń w jednym procesie.

Co dalej?

Pogratuluj sobie! Ukończyłeś kurs Nextflow dla Genomiki.

Przejdź do końcowego podsumowania kursu, aby przejrzeć to, czego się nauczyłeś i dowiedzieć się, co będzie dalej.