भाग 2: कोहोर्ट पर संयुक्त कॉलिंग¶

AI-सहायता प्राप्त अनुवाद - अधिक जानें और सुधार सुझाएं

इस कोर्स के पहले भाग में, आपने एक वेरिएंट कॉलिंग पाइपलाइन बनाई जो पूरी तरह से रेखीय थी और प्रत्येक नमूने के डेटा को दूसरों से स्वतंत्र रूप से प्रोसेस करती थी। हालाँकि, वास्तविक जीनोमिक्स उपयोग में, आपको आमतौर पर एक साथ कई नमूनों के वेरिएंट कॉल देखने की आवश्यकता होगी।

इस दूसरे भाग में, हम आपको दिखाते हैं कि भाग 1 की पाइपलाइन पर आधारित, GATK के साथ संयुक्त वेरिएंट कॉलिंग को लागू करने के लिए channels और channel operators का उपयोग कैसे करें।

विधि अवलोकन¶

GATK वेरिएंट कॉलिंग विधि जिसका हमने इस कोर्स के पहले भाग में उपयोग किया था, केवल प्रति नमूना वेरिएंट कॉल उत्पन्न करती थी। यदि आप केवल प्रत्येक नमूने से वेरिएंट को अलगाव में देखना चाहते हैं तो यह ठीक है, लेकिन यह सीमित जानकारी देता है। कई नमूनों में वेरिएंट कॉल कैसे भिन्न होते हैं यह देखना अक्सर अधिक दिलचस्प होता है, और ऐसा करने के लिए, GATK एक वैकल्पिक विधि प्रदान करता है जिसे संयुक्त वेरिएंट कॉलिंग कहा जाता है, जिसे हम यहाँ प्रदर्शित करते हैं।

संयुक्त वेरिएंट कॉलिंग में प्रत्येक नमूने के लिए एक विशेष प्रकार का वेरिएंट आउटपुट उत्पन्न करना शामिल है जिसे GVCF (Genomic VCF के लिए) कहा जाता है, फिर सभी नमूनों से GVCF डेटा को जोड़ना और अंत में, 'संयुक्त जीनोटाइपिंग' सांख्यिकीय विश्लेषण चलाना।

संयुक्त विश्लेषण

नमूने के GVCF में विशेष यह है कि इसमें जीनोम के लक्षित क्षेत्र में सभी स्थितियों के बारे में अनुक्रम डेटा सांख्यिकी का सारांश देने वाले रिकॉर्ड होते हैं, न कि केवल उन स्थितियों के जहाँ प्रोग्राम को भिन्नता के साक्ष्य मिले। यह संयुक्त जीनोटाइपिंग गणना के लिए महत्वपूर्ण है (अधिक पढ़ें)।

GVCF को GATK HaplotypeCaller द्वारा उत्पन्न किया जाता है, वही टूल जिसका हमने भाग 1 में उपयोग किया था, एक अतिरिक्त पैरामीटर (-ERC GVCF) के साथ। GVCFs को जोड़ना GATK GenomicsDBImport के साथ किया जाता है, जो प्रति-नमूना कॉल को एक डेटा स्टोर (डेटाबेस के समान) में जोड़ता है, फिर वास्तविक 'संयुक्त जीनोटाइपिंग' विश्लेषण GATK GenotypeGVCFs के साथ किया जाता है।

वर्कफ़्लो¶

तो संक्षेप में, कोर्स के इस भाग में, हम एक वर्कफ़्लो विकसित करने जा रहे हैं जो निम्नलिखित करता है:

Samtools का उपयोग करके प्रत्येक BAM इनपुट फ़ाइल के लिए एक इंडेक्स फ़ाइल उत्पन्न करें
प्रति-नमूना जीनोमिक वेरिएंट कॉल का GVCF उत्पन्न करने के लिए प्रत्येक BAM इनपुट फ़ाइल पर GATK HaplotypeCaller चलाएँ
सभी GVCFs एकत्र करें और उन्हें GenomicsDB डेटा स्टोर में जोड़ें
कोहोर्ट-स्तरीय VCF उत्पन्न करने के लिए संयुक्त GVCF डेटा स्टोर पर संयुक्त जीनोटाइपिंग चलाएँ

हम इसे भाग 1 के समान डेटासेट पर लागू करेंगे।

0. वार्मअप: Samtools और GATK को सीधे चलाएँ¶

पहले की तरह, हम उन्हें वर्कफ़्लो में लपेटने का प्रयास करने से पहले कमांड को मैन्युअल रूप से आज़माना चाहते हैं।

Note

सुनिश्चित करें कि आप सही कार्य डायरेक्टरी में हैं: cd /workspaces/training/nf4-science/genomics

0.1. Samtools के साथ BAM इनपुट फ़ाइल को इंडेक्स करें¶

यह पहला चरण भाग 1 के समान है, इसलिए यह बहुत परिचित लगना चाहिए, लेकिन इस बार हमें तीनों नमूनों के लिए ऐसा करना होगा।

Note

हमने तकनीकी रूप से अपनी पाइपलाइन के माध्यम से तीनों नमूनों के लिए पहले से ही इंडेक्स फ़ाइलें उत्पन्न की हैं, इसलिए हम उन्हें results डायरेक्टरी से निकाल सकते हैं। हालाँकि, इसे मैन्युअल रूप से फिर से करना साफ़-सुथरा है, और इसमें केवल एक मिनट लगेगा।

0.1.1. Samtools कंटेनर को इंटरैक्टिव रूप से स्पिन करें¶

docker run -it -v ./data:/data community.wave.seqera.io/library/samtools:1.20--b5dfbd93de237464

0.1.2. तीनों नमूनों के लिए इंडेक्सिंग कमांड चलाएँ¶

samtools index /data/bam/reads_mother.bam
samtools index /data/bam/reads_father.bam
samtools index /data/bam/reads_son.bam

पहले की तरह, यह संबंधित BAM फ़ाइलों के समान डायरेक्टरी में इंडेक्स फ़ाइलें उत्पन्न करनी चाहिए।

डायरेक्टरी सामग्री

data/bam/
├── reads_father.bam
├── reads_father.bam.bai
├── reads_mother.bam
├── reads_mother.bam.bai
├── reads_son.bam
└── reads_son.bam.bai

अब जब हमारे पास तीनों नमूनों के लिए इंडेक्स फ़ाइलें हैं, तो हम उनमें से प्रत्येक के लिए GVCFs उत्पन्न करने के लिए आगे बढ़ सकते हैं।

0.1.3. Samtools कंटेनर से बाहर निकलें¶

exit

0.2. GVCF मोड में GATK HaplotypeCaller के साथ वेरिएंट कॉल करें¶

यह दूसरा चरण उस चीज़ के समान है जो हमने भाग 1: Hello Genomics में किया था, लेकिन हम अब GATK को 'GVCF मोड' में चलाने जा रहे हैं।

0.2.1. GATK कंटेनर को इंटरैक्टिव रूप से स्पिन करें¶

docker run -it -v ./data:/data community.wave.seqera.io/library/gatk4:4.5.0.0--730ee8817e436867

0.2.2. GVCF विकल्प के साथ वेरिएंट कॉलिंग कमांड चलाएँ¶

एक जीनोमिक VCF (GVCF) उत्पन्न करने के लिए, हम बेस कमांड में -ERC GVCF विकल्प जोड़ते हैं, जो HaplotypeCaller के GVCF मोड को चालू करता है।

हम आउटपुट फ़ाइल के लिए फ़ाइल एक्सटेंशन को .vcf से .g.vcf में भी बदलते हैं। यह तकनीकी रूप से आवश्यकता नहीं है, लेकिन यह एक दृढ़ता से अनुशंसित परंपरा है।

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

यह कंटेनर में वर्तमान कार्य डायरेक्टरी में GVCF आउटपुट फ़ाइल reads_mother.g.vcf बनाता है।

यदि आप सामग्री देखने के लिए इसे cat करते हैं, तो आप देखेंगे कि यह भाग 1 में हमारे द्वारा उत्पन्न समतुल्य VCF की तुलना में बहुत लंबा है। आप फ़ाइल की शुरुआत तक भी स्क्रॉल नहीं कर सकते, और अधिकांश पंक्तियाँ भाग 1 में VCF में हमने जो देखा था उससे काफी अलग दिखती हैं।

Output
20_10037292_10066351    14714   .       T       <NON_REF>       .       .       END=14718       GT:DP:GQ:MIN_DP:PL       0/0:37:99:37:0,99,1192
20_10037292_10066351    14719   .       T       <NON_REF>       .       .       END=14719       GT:DP:GQ:MIN_DP:PL       0/0:36:82:36:0,82,1087
20_10037292_10066351    14720   .       T       <NON_REF>       .       .       END=14737       GT:DP:GQ:MIN_DP:PL       0/0:42:99:37:0,100,1160

ये गैर-वेरिएंट क्षेत्रों का प्रतिनिधित्व करते हैं जहाँ वेरिएंट कॉलर को भिन्नता का कोई साक्ष्य नहीं मिला, इसलिए इसने भिन्नता की अनुपस्थिति में अपने विश्वास के स्तर का वर्णन करने वाली कुछ सांख्यिकी कैप्चर कीं। यह दो बहुत अलग केस आंकड़ों के बीच अंतर करना संभव बनाता है: (1) अच्छी गुणवत्ता का डेटा है जो दर्शाता है कि नमूना homozygous-reference है, और (2) किसी भी तरह से निर्धारण करने के लिए पर्याप्त अच्छा डेटा उपलब्ध नहीं है।

GVCF में, आमतौर पर इस तरह की बहुत सारी गैर-वेरिएंट लाइनें होती हैं, जिनके बीच कम संख्या में वेरिएंट रिकॉर्ड बिखरे होते हैं। वास्तविक वेरिएंट कॉल खोजने के लिए फ़ाइल की केवल पहली 176 लाइनों को लोड करने के लिए GVCF पर head -176 चलाने का प्रयास करें।

Output
20_10037292_10066351    3479    .       T       <NON_REF>       .       .       END=3479        GT:DP:GQ:MIN_DP:PL       0/0:34:36:34:0,36,906
20_10037292_10066351    3480    .       C       CT,<NON_REF>    503.03  .       DP=23;ExcessHet=0.0000;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQandDP=82800,23    GT:AD:DP:GQ:PL:SB       1/1:0,18,0:18:54:517,54,0,517,54,517:0,0,7,11
20_10037292_10066351    3481    .       T       <NON_REF>       .       .       END=3481        GT:DP:GQ:MIN_DP:PL       0/0:21:51:21:0,51,765

दूसरी लाइन फ़ाइल में पहला वेरिएंट रिकॉर्ड दिखाती है, जो भाग 1 में हमने देखी VCF फ़ाइल में पहले वेरिएंट से मेल खाती है।

मूल VCF की तरह, आउटपुट GVCF फ़ाइल भी एक इंडेक्स फ़ाइल के साथ आती है, जिसे reads_mother.g.vcf.idx कहा जाता है।

0.2.3. अन्य दो नमूनों पर प्रक्रिया दोहराएँ¶

संयुक्त जीनोटाइपिंग चरण का परीक्षण करने के लिए, हमें तीनों नमूनों के लिए GVCFs की आवश्यकता है, तो आइए अभी मैन्युअल रूप से उन्हें उत्पन्न करें।

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

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

एक बार यह पूरा हो जाने पर, आपके वर्तमान डायरेक्टरी में .g.vcf में समाप्त होने वाली तीन फ़ाइलें (प्रति नमूना एक) और .g.vcf.idx में समाप्त होने वाली उनकी संबंधित इंडेक्स फ़ाइलें होनी चाहिए।

0.3. संयुक्त जीनोटाइपिंग चलाएँ¶

अब जब हमारे पास सभी GVCFs हैं, तो हम अंततः नमूनों के कोहोर्ट के लिए वेरिएंट कॉल उत्पन्न करने के लिए संयुक्त जीनोटाइपिंग दृष्टिकोण को आज़मा सकते हैं। एक अनुस्मारक के रूप में, यह एक दो-चरण विधि है जिसमें सभी GVCFs से डेटा को डेटा स्टोर में जोड़ना, फिर संयुक्त-कॉल किए गए वेरिएंट के अंतिम VCF को उत्पन्न करने के लिए संयुक्त जीनोटाइपिंग विश्लेषण उचित रूप से चलाना शामिल है।

0.3.1. सभी प्रति-नमूना GVCFs को जोड़ें¶

यह पहला चरण एक अन्य GATK टूल का उपयोग करता है, जिसे GenomicsDBImport कहा जाता है, सभी GVCFs से डेटा को GenomicsDB डेटा स्टोर में जोड़ने के लिए।

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

इस चरण का आउटपुट प्रभावी रूप से एक डायरेक्टरी है जिसमें कई अलग-अलग फ़ाइलों के रूप में संयुक्त वेरिएंट डेटा रखने वाली आगे नेस्टेड डायरेक्टरियों का एक सेट होता है। आप इसके चारों ओर देख सकते हैं लेकिन आप जल्दी देखेंगे कि इस डेटा स्टोर प्रारूप को मनुष्यों द्वारा सीधे पढ़ने के लिए नहीं बनाया गया है।

Note

GATK में ऐसे टूल शामिल हैं जो आवश्यकतानुसार डेटा स्टोर से वेरिएंट कॉल डेटा का निरीक्षण और निष्कर्षण करना संभव बनाते हैं।

0.3.2. संयुक्त जीनोटाइपिंग विश्लेषण उचित रूप से चलाएँ¶

यह दूसरा चरण एक और GATK टूल का उपयोग करता है, जिसे GenotypeGVCFs कहा जाता है, कोहोर्ट में सभी नमूनों में उपलब्ध डेटा के आलोक में वेरिएंट सांख्यिकी और व्यक्तिगत जीनोटाइप की पुनर्गणना करने के लिए।

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

यह कंटेनर में वर्तमान कार्य डायरेक्टरी में VCF आउटपुट फ़ाइल family_trio.vcf बनाता है। यह एक और उचित रूप से छोटी फ़ाइल है इसलिए आप इसकी सामग्री देखने के लिए इस फ़ाइल को cat कर सकते हैं, और पहली कुछ वेरिएंट लाइनों को खोजने के लिए ऊपर स्क्रॉल कर सकते हैं।

family_trio.vcf
#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

यह भाग 1 में हमारे द्वारा उत्पन्न मूल VCF की तरह अधिक दिखता है, सिवाय इसके कि इस बार हमारे पास तीनों नमूनों के लिए जीनोटाइप-स्तरीय जानकारी है। फ़ाइल में अंतिम तीन कॉलम नमूनों के लिए जीनोटाइप ब्लॉक हैं, जो वर्णानुक्रम में सूचीबद्ध हैं।

यदि हम अपने परीक्षण परिवार ट्रायो के लिए पहले वेरिएंट के लिए कॉल किए गए जीनोटाइप देखते हैं, तो हम देखते हैं कि पिता heterozygous-variant (0/1) हैं, और माँ और बेटा दोनों homozygous-variant (1/1) हैं।

यह अंततः वह जानकारी है जिसे हम डेटासेट से निकालने की कोशिश कर रहे हैं! तो आइए इन सभी को Nextflow वर्कफ़्लो में लपेटें ताकि हम इसे बड़े पैमाने पर कर सकें।

0.3.3. GATK कंटेनर से बाहर निकलें¶

exit

निष्कर्ष¶

आप जानते हैं कि टर्मिनल में संयुक्त वेरिएंट कॉलिंग से जुड़े व्यक्तिगत कमांड को कैसे चलाना है ताकि यह सत्यापित किया जा सके कि वे आपकी इच्छित जानकारी उत्पन्न करेंगे।

आगे क्या है?¶

इन कमांड को वास्तविक पाइपलाइन में लपेटें।

1. GVCF उत्पन्न करने के लिए प्रति-नमूना वेरिएंट कॉलिंग चरण को संशोधित करें¶

अच्छी खबर यह है कि हमें सब कुछ फिर से शुरू करने की आवश्यकता नहीं है, क्योंकि हमने पहले से ही भाग 1 में एक वर्कफ़्लो लिखा है जो इस काम का कुछ हिस्सा करता है। हालाँकि, वह पाइपलाइन VCF फ़ाइलें उत्पन्न करती है, जबकि अब हम संयुक्त जीनोटाइपिंग करने के लिए GVCF फ़ाइलें चाहते हैं। इसलिए हमें GVCF वेरिएंट कॉलिंग मोड को चालू करके और आउटपुट फ़ाइल एक्सटेंशन को अपडेट करके शुरू करने की आवश्यकता है।

Note

सुविधा के लिए, हम GATK वर्कफ़्लो की एक नई प्रतिलिपि के साथ काम करने जा रहे हैं जैसा कि यह भाग 1 के अंत में खड़ा है, लेकिन एक अलग नाम के तहत: genomics-2.nf।

1.1. HaplotypeCaller को GVCF उत्सर्जित करने के लिए कहें और आउटपुट एक्सटेंशन अपडेट करें¶

आइए कोड एडिटर में genomics-2.nf फ़ाइल खोलें। यह बहुत परिचित लगनी चाहिए, लेकिन यदि आप खुद को संतुष्ट करना चाहते हैं कि यह अपेक्षित रूप से चलती है तो इसे चलाने के लिए स्वतंत्र महसूस करें।

हम दो परिवर्तन करके शुरू करने जा रहे हैं:

GATK HaplotypeCaller कमांड में -ERC GVCF पैरामीटर जोड़ें;
GATK परंपरा के अनुसार, संबंधित .g.vcf एक्सटेंशन का उपयोग करने के लिए आउटपुट फ़ाइल पथ अपडेट करें।

जब आप -ERC GVCF जोड़ें तो सुनिश्चित करें कि आप पिछली लाइन के अंत में बैकस्लैश (\) जोड़ें।

बाद मेंपहले

genomics-2.nf
    """
    gatk HaplotypeCaller \
        -R ${ref_fasta} \
        -I ${input_bam} \
        -O ${input_bam}.g.vcf \
        -L ${interval_list} \
        -ERC GVCF
    """

genomics-2.nf
    """
    gatk HaplotypeCaller \
        -R ${ref_fasta} \
        -I ${input_bam} \
        -O ${input_bam}.vcf \
        -L ${interval_list}
    """

और VCFs के बजाय GVCFs उत्पन्न करने के लिए HaplotypeCaller को स्विच करने के लिए बस इतना ही है, है ना?

1.2. सत्यापित करने के लिए पाइपलाइन चलाएँ कि आप GVCFs उत्पन्न कर सकते हैं¶

Nextflow निष्पादन कमांड पहले की तरह ही है, वर्कफ़्लो फ़ाइल नाम को छोड़कर। सुनिश्चित करें कि आप इसे उचित रूप से अपडेट करें।

nextflow run genomics-2.nf

कमांड आउटपुट

N E X T F L O W   ~  version 25.10.2

┃ Launching `genomics-2.nf` [crazy_venter] DSL2 - revision: a2d6f6f09f

executor >  local (6)
[f1/8d8486] SAMTOOLS_INDEX (1)       | 3 of 3 ✔
[72/3249ca] GATK_HAPLOTYPECALLER (3) | 0 of 3
ERROR ~ Error executing process > 'GATK_HAPLOTYPECALLER (2)'

Caused by:
  Missing output file(s) `reads_son.bam.vcf` expected by process `GATK_HAPLOTYPECALLER (2)`

Command executed:

  gatk HaplotypeCaller         -R ref.fasta         -I reads_son.bam         -O reads_son.bam.g.vcf         -L intervals.bed         -ERC GVCF

और आउटपुट है... सब लाल! ओह नहीं।

जो कमांड निष्पादित की गई थी वह सही है, इसलिए हम सही थे कि GATK टूल के व्यवहार को बदलने के लिए यह पर्याप्त था। लेकिन लापता आउटपुट फ़ाइल के बारे में उस लाइन को देखें। कुछ भी नोटिस करें?

यह सही है, हम Nextflow को बताना भूल गए कि अपेक्षित आउटपुट फ़ाइल नाम बदल गया है। उफ़्फ़।

1.3. process outputs ब्लॉक में भी आउटपुट फ़ाइल एक्सटेंशन अपडेट करें¶

क्योंकि टूल कमांड में ही फ़ाइल एक्सटेंशन बदलना पर्याप्त नहीं है, आपको Nextflow को यह भी बताना होगा कि अपेक्षित आउटपुट फ़ाइल नाम बदल गया है।

बाद मेंपहले

genomics-2.nf
    output:
    path "${input_bam}.g.vcf"     , emit: vcf
    path "${input_bam}.g.vcf.idx" , emit: idx

genomics-2.nf
    output:
    path "${input_bam}.vcf"     , emit: vcf
    path "${input_bam}.vcf.idx" , emit: idx

1.4. नए GVCF आउटपुट के लिए प्रकाशन लक्ष्य अपडेट करें¶

चूँकि हम अब VCFs के बजाय GVCFs उत्पन्न कर रहे हैं, हमें अधिक वर्णनात्मक नामों का उपयोग करने के लिए वर्कफ़्लो के publish: अनुभाग को अपडेट करना चाहिए। हम स्पष्टता के लिए GVCF फ़ाइलों को उनकी अपनी उपनिर्देशिका में भी व्यवस्थित करेंगे।

बाद मेंपहले

genomics-2.nf
    publish:
    indexed_bam = SAMTOOLS_INDEX.out
    gvcf = GATK_HAPLOTYPECALLER.out.vcf
    gvcf_idx = GATK_HAPLOTYPECALLER.out.idx

genomics-2.nf
    publish:
    indexed_bam = SAMTOOLS_INDEX.out
    vcf = GATK_HAPLOTYPECALLER.out.vcf
    vcf_idx = GATK_HAPLOTYPECALLER.out.idx

1.5. नई डायरेक्टरी संरचना के लिए आउटपुट ब्लॉक अपडेट करें¶

हमें GVCF फ़ाइलों को gvcf उपनिर्देशिका में रखने के लिए output ब्लॉक को भी अपडेट करने की आवश्यकता है।

बाद मेंपहले

genomics-2.nf
output {
    indexed_bam {
        path 'indexed_bam'
    }
    gvcf {
        path 'gvcf'
    }
    gvcf_idx {
        path 'gvcf'
    }
}

genomics-2.nf
output {
    indexed_bam {
        path '.'
    }
    vcf {
        path '.'
    }
    vcf_idx {
        path '.'
    }
}

1.6. पाइपलाइन फिर से चलाएँ¶

आइए इस बार इसे -resume के साथ चलाएँ।

nextflow run genomics-2.nf -resume

कमांड आउटपुट

N E X T F L O W   ~  version 25.10.2

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

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

इस बार यह काम करता है।

Nextflow आउटपुट स्वयं (सामान्य VCF मोड में सफल रन की तुलना में) कोई अलग नहीं दिखता है, लेकिन अब हम .g.vcf फ़ाइलें और उनकी संबंधित इंडेक्स फ़ाइलें, तीनों नमूनों के लिए, उपनिर्देशिकाओं में व्यवस्थित पा सकते हैं।

डायरेक्टरी सामग्री (symlinks छोटी की गई)

results_genomics/
├── 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

यदि आप GVCF फ़ाइलों में से एक को खोलते हैं और इसके माध्यम से स्क्रॉल करते हैं, तो आप सत्यापित कर सकते हैं कि GATK HaplotypeCaller ने अनुरोध के अनुसार GVCF फ़ाइलें उत्पन्न कीं।

निष्कर्ष¶

ठीक है, यह Nextflow सीखने के मामले में न्यूनतम था... लेकिन process output ब्लॉक के महत्व को दोहराने का यह एक अच्छा अवसर था!

आगे क्या है?¶

सभी नमूनों में GVCF डेटा एकत्र करना और जोड़ना सीखें।

2. सभी नमूनों में GVCF डेटा एकत्र करें और जोड़ें¶

अब हमें सभी प्रति-नमूना GVCFs से डेटा को एक ऐसे रूप में जोड़ने की आवश्यकता है जो हम जो संयुक्त जीनोटाइपिंग विश्लेषण करना चाहते हैं उसका समर्थन करता है।

2.1. उस process को परिभाषित करें जो GVCFs को जोड़ेगा¶

वार्मअप अनुभाग में हमने पहले जो किया था उसकी याद दिलाने के लिए, GVCFs को जोड़ना GATK टूल GenomicsDBImport के लिए एक काम है, जो तथाकथित GenomicsDB प्रारूप में एक डेटा स्टोर उत्पन्न करेगा।

आइए वार्मअप अनुभाग में पहले उपयोग की गई कमांड के आधार पर यह परिभाषित करने के लिए एक नई प्रक्रिया लिखें कि यह कैसे काम करने वाला है।

genomics-2.nf
/*
 * GVCFs को GenomicsDB datastore में जोड़ें
 */
process GATK_GENOMICSDB {

    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

    output:
    path "${cohort_name}_gdb"

    script:
    """
    gatk GenomicsDBImport \
        -V ${all_gvcfs} \
        -L ${interval_list} \
        --genomicsdb-workspace-path ${cohort_name}_gdb
    """
}

आपको क्या लगता है, उचित लग रहा है?

आइए इसे वायर करें और देखें कि क्या होता है।

2.2. डिफ़ॉल्ट मान के साथ `cohort_name` पैरामीटर जोड़ें¶

हमें कोहोर्ट के लिए एक मनमाना नाम प्रदान करने की आवश्यकता है। प्रशिक्षण श्रृंखला में बाद में आप इस तरह की चीज़ के लिए नमूना मेटाडेटा का उपयोग करना सीखेंगे, लेकिन अभी के लिए हम सुविधा के लिए params का उपयोग करके एक CLI पैरामीटर घोषित करते हैं और इसे एक डिफ़ॉल्ट मान देते हैं।

genomics-2.nf
    // अंतिम आउटपुट फ़ाइल के लिए आधार नाम
    cohort_name: String = "family_trio"

2.3. नमूनों में GATK_HAPLOTYPECALLER के आउटपुट एकत्र करें¶

यदि हम GATK_HAPLOTYPECALLER process के आउटपुट channel को जैसा है वैसा ही प्लग करते हैं, तो Nextflow प्रत्येक नमूना GVCF पर अलग से process को कॉल करेगा। हालाँकि, हम सभी तीन GVCFs (और उनकी इंडेक्स फ़ाइलों) को इस तरह से बंडल करना चाहते हैं कि Nextflow उन सभी को एक साथ एक process कॉल को सौंपे।

अच्छी खबर: हम collect() channel operator का उपयोग करके ऐसा कर सकते हैं। आइए GATK_HAPLOTYPECALLER की कॉल के ठीक बाद, workflow बॉडी में निम्नलिखित लाइनें जोड़ें:

genomics-2.nf
// नमूनों में वेरिएंट कॉलिंग आउटपुट एकत्र करें
all_gvcfs_ch = GATK_HAPLOTYPECALLER.out.vcf.collect()
all_idxs_ch = GATK_HAPLOTYPECALLER.out.idx.collect()

क्या यह थोड़ा जटिल लगता है? आइए इसे तोड़ें और इसे सादी भाषा में अनुवादित करें।

हम GATK_HAPLOTYPECALLER process से आउटपुट channel ले रहे हैं, जिसे .out गुण का उपयोग करके संदर्भित किया गया है।
channel से निकलने वाला प्रत्येक 'तत्व' फ़ाइलों की एक जोड़ी है: GVCF और उसकी इंडेक्स फ़ाइल, उस क्रम में क्योंकि वह क्रम है जिसमें वे process output ब्लॉक में सूचीबद्ध हैं। सुविधाजनक रूप से, क्योंकि पिछले सत्र में हमने इस process के आउटपुट का नाम दिया था (emit: का उपयोग करके), हम एक हाथ में .out गुण के बाद .vcf जोड़कर GVCFs को चुन सकते हैं और दूसरी ओर .idx जोड़कर इंडेक्स फ़ाइलों को। यदि हमने उन आउटपुट का नाम नहीं दिया होता, तो हमें उन्हें क्रमशः .out[0] और .out[1] द्वारा संदर्भित करना होता।
हम सभी GVCF फ़ाइलों को all_gvcfs_ch नामक एक नए channel में एक एकल तत्व में बंडल करने के लिए collect() channel operator जोड़ते हैं, और all_idxs_ch नामक नया channel बनाने के लिए इंडेक्स फ़ाइलों के साथ भी ऐसा ही करते हैं।

Tip

यदि आपको यह कल्पना करना कठिन हो रहा है कि यहाँ वास्तव में क्या हो रहा है, तो याद रखें कि आप channel operators को लागू करने से पहले और बाद में channels की सामग्री का निरीक्षण करने के लिए view() operator का उपयोग कर सकते हैं।

परिणामी all_gvcfs_ch और all_idxs_ch channels वे हैं जिन्हें हम अभी लिखे गए GATK_GENOMICSDB process में प्लग करने जा रहे हैं।

Note

यदि आप सोच रहे थे, तो हम GVCFs और उनकी इंडेक्स फ़ाइलों को अलग से एकत्र करते हैं क्योंकि GATK GenomicsDBImport कमांड केवल GVCF फ़ाइल पथ देखना चाहता है। सौभाग्य से, चूँकि Nextflow निष्पादन के लिए सभी फ़ाइलों को एक साथ मंचित करेगा, इसलिए हमें भाग 1 में BAMs और उनके इंडेक्स के लिए जैसी फ़ाइलों के क्रम के बारे में चिंता करने की आवश्यकता नहीं है।

2.4. GATK_GENOMICSDB चलाने के लिए वर्कफ़्लो ब्लॉक में एक कॉल जोड़ें¶

हमारे पास एक process है, और हमारे पास इनपुट channels हैं। हमें बस process कॉल जोड़ने की आवश्यकता है।

genomics-2.nf
    // GVCFs को GenomicsDB डेटा स्टोर में जोड़ें
    GATK_GENOMICSDB(
        all_gvcfs_ch,
        all_idxs_ch,
        intervals_file,
        params.cohort_name
    )

ठीक है, सब कुछ वायर किया गया है।

2.5. वर्कफ़्लो चलाएँ¶

आइए देखें कि क्या यह काम करता है।

nextflow run genomics-2.nf -resume

कमांड आउटपुट

N E X T F L O W   ~  version 25.10.2

┃ Launching `genomics-2.nf` [disturbed_bell] DSL2 - revision: 57942246cc

executor >  local (1)
[f1/8d8486] SAMTOOLS_INDEX (1)       | 3 of 3, cached: 3 ✔
[e4/4ed55e] GATK_HAPLOTYPECALLER (2) | 3 of 3, cached: 3 ✔
[51/d350ea] GATK_GENOMICSDB          | 0 of 1
ERROR ~ Error executing process > 'GATK_GENOMICSDB'

Caused by:
  Process `GATK_GENOMICSDB` terminated with an error exit status (1)

Command executed:

  gatk GenomicsDBImport         -V reads_son.bam.g.vcf reads_father.bam.g.vcf reads_mother.bam.g.vcf         -L intervals.bed         --genomicsdb-workspace-path family_trio_gdb

यह काफी तेज़ी से चलता है, क्योंकि हम -resume के साथ चल रहे हैं, लेकिन यह विफल हो जाता है!

आह। सकारात्मक पक्ष पर, हम देखते हैं कि Nextflow ने GATK_GENOMICSDB process को उठाया है, और विशेष रूप से इसे केवल एक बार कॉल किया है। यह सुझाव देता है कि collect() दृष्टिकोण काम किया, एक सीमा तक। लेकिन, और यह एक बड़ा है, process कॉल विफल रहा।

जब हम ऊपर कंसोल आउटपुट में खोदते हैं, तो हम देख सकते हैं कि निष्पादित कमांड सही नहीं है।

क्या आप त्रुटि देख सकते हैं? इस बिट को देखें: -V reads_father.bam.g.vcf reads_son.bam.g.vcf reads_mother.bam.g.vcf

हमने gatk GenomicsDBImport को एकल -V तर्क के लिए कई GVCF फ़ाइलें दीं, लेकिन टूल प्रत्येक GVCF फ़ाइल के लिए एक अलग -V तर्क की अपेक्षा करता है।

एक अनुस्मारक के रूप में, यह वह कमांड थी जो हमने कंटेनर में चलाई थी:

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

तो इसका मतलब है कि हमें किसी तरह GVCF फ़ाइलों के अपने बंडल को उचित रूप से स्वरूपित कमांड स्ट्रिंग में बदलने की आवश्यकता है।

2.6. प्रत्येक इनपुट GVCF के लिए एक अलग `-V` तर्क के साथ एक कमांड लाइन बनाएँ¶

यह वह जगह है जहाँ Groovy पर आधारित Nextflow होना काम आता है, क्योंकि यह हमें आवश्यक कमांड स्ट्रिंग बनाने के लिए कुछ काफी सीधे स्ट्रिंग हेरफेर का उपयोग करने की अनुमति देने वाला है।

विशेष रूप से, इस सिंटैक्स का उपयोग करते हुए: all_gvcfs.collect { gvcf -> "-V ${gvcf}" }.join(' ')

एक बार फिर, आइए इसे इसके घटकों में तोड़ें।

सबसे पहले, हम all_gvcfs input channel की सामग्री लेते हैं और उस पर .collect() लागू करते हैं (पहले की तरह)।
यह हमें बंडल में प्रत्येक व्यक्तिगत GVCF फ़ाइल पथ को closure में पास करने की अनुमति देता है, { gvcf -> "-V ${gvcf}" }, जहाँ gvcf उस GVCF फ़ाइल पथ को संदर्भित करता है। closure एक मिनी-फ़ंक्शन है जिसका उपयोग हम फ़ाइल पथ में -V जोड़ने के लिए करते हैं, "-V ${gvcf}" के रूप में।
फिर हम सभी तीन स्ट्रिंग्स को विभाजक के रूप में एकल स्थान के साथ संयोजित करने के लिए .join(' ') का उपयोग करते हैं।

एक ठोस उदाहरण के साथ, यह इस तरह दिखता है:

हमारे पास तीन फ़ाइलें हैं:

[A.ext, B.ext, C.ext]

closure प्रत्येक को स्ट्रिंग्स बनाने के लिए संशोधित करता है:

"-V A.ext", "-V B.ext", "-V C.ext"

.join(' ') ऑपरेशन अंतिम स्ट्रिंग उत्पन्न करता है:

"-V A.ext -V B.ext -V C.ext"

एक बार जब हमारे पास वह स्ट्रिंग हो, तो हम इसे एक स्थानीय वेरिएबल, gvcfs_line को असाइन कर सकते हैं, जिसे def कीवर्ड के साथ परिभाषित किया गया है:

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

ठीक है, तो हमारे पास हमारी स्ट्रिंग हेरफेर चीज़ है। हम इसे कहाँ रखते हैं?

हम चाहते हैं कि यह process परिभाषा में कहीं जाए, क्योंकि हम चाहते हैं कि यह बाद में करें जब हमने GVCF फ़ाइल पथों को process में channeled किया हो। ऐसा इसलिए है क्योंकि Nextflow को उन्हें फ़ाइल पथों के रूप में देखना चाहिए ताकि निष्पादन के लिए फ़ाइलों को स्वयं सही ढंग से मंचित किया जा सके।

लेकिन process में कहाँ

भाग 2: कोहोर्ट पर संयुक्त कॉलिंग¶

विधि अवलोकन¶

वर्कफ़्लो¶

0. वार्मअप: Samtools और GATK को सीधे चलाएँ¶

0.1. Samtools के साथ BAM इनपुट फ़ाइल को इंडेक्स करें¶

0.1.1. Samtools कंटेनर को इंटरैक्टिव रूप से स्पिन करें¶

0.1.2. तीनों नमूनों के लिए इंडेक्सिंग कमांड चलाएँ¶

0.1.3. Samtools कंटेनर से बाहर निकलें¶

0.2. GVCF मोड में GATK HaplotypeCaller के साथ वेरिएंट कॉल करें¶

0.2.1. GATK कंटेनर को इंटरैक्टिव रूप से स्पिन करें¶

0.2.2. GVCF विकल्प के साथ वेरिएंट कॉलिंग कमांड चलाएँ¶

0.2.3. अन्य दो नमूनों पर प्रक्रिया दोहराएँ¶

0.3. संयुक्त जीनोटाइपिंग चलाएँ¶

0.3.1. सभी प्रति-नमूना GVCFs को जोड़ें¶

0.3.2. संयुक्त जीनोटाइपिंग विश्लेषण उचित रूप से चलाएँ¶

0.3.3. GATK कंटेनर से बाहर निकलें¶

निष्कर्ष¶

आगे क्या है?¶

1. GVCF उत्पन्न करने के लिए प्रति-नमूना वेरिएंट कॉलिंग चरण को संशोधित करें¶

1.1. HaplotypeCaller को GVCF उत्सर्जित करने के लिए कहें और आउटपुट एक्सटेंशन अपडेट करें¶

1.2. सत्यापित करने के लिए पाइपलाइन चलाएँ कि आप GVCFs उत्पन्न कर सकते हैं¶

1.3. process outputs ब्लॉक में भी आउटपुट फ़ाइल एक्सटेंशन अपडेट करें¶

1.4. नए GVCF आउटपुट के लिए प्रकाशन लक्ष्य अपडेट करें¶

1.5. नई डायरेक्टरी संरचना के लिए आउटपुट ब्लॉक अपडेट करें¶

1.6. पाइपलाइन फिर से चलाएँ¶

निष्कर्ष¶

आगे क्या है?¶

2. सभी नमूनों में GVCF डेटा एकत्र करें और जोड़ें¶

2.1. उस process को परिभाषित करें जो GVCFs को जोड़ेगा¶

2.2. डिफ़ॉल्ट मान के साथ cohort_name पैरामीटर जोड़ें¶

2.3. नमूनों में GATK_HAPLOTYPECALLER के आउटपुट एकत्र करें¶

2.4. GATK_GENOMICSDB चलाने के लिए वर्कफ़्लो ब्लॉक में एक कॉल जोड़ें¶

2.5. वर्कफ़्लो चलाएँ¶

2.6. प्रत्येक इनपुट GVCF के लिए एक अलग -V तर्क के साथ एक कमांड लाइन बनाएँ¶

2.2. डिफ़ॉल्ट मान के साथ `cohort_name` पैरामीटर जोड़ें¶

2.6. प्रत्येक इनपुट GVCF के लिए एक अलग `-V` तर्क के साथ एक कमांड लाइन बनाएँ¶