भाग 3: कोहोर्ट पर जॉइंट कॉलिंग¶
भाग 2 में, तुमने एक per-sample वेरिएंट कॉलिंग पाइपलाइन बनाई थी जो हर सैंपल का डेटा स्वतंत्र रूप से प्रोसेस करती है। अब हम इसे एक्सटेंड करके जॉइंट वेरिएंट कॉलिंग को इम्प्लीमेंट करने जा रहे हैं, जैसा कि भाग 1 में बताया गया था।
असाइनमेंट¶
कोर्स के इस भाग में, हम वर्कफ़्लो को निम्नलिखित कार्यों के लिए एक्सटेंड करने जा रहे हैं:
- हर BAM इनपुट फ़ाइल के लिए Samtools का उपयोग करके एक इंडेक्स फ़ाइल जेनरेट करो
- हर BAM इनपुट फ़ाइल पर GATK HaplotypeCaller को रन करो ताकि per-sample जीनोमिक वेरिएंट कॉल्स का एक GVCF जेनरेट हो
- सभी GVCFs को इकट्ठा करो और उन्हें एक GenomicsDB डेटा स्टोर में कंबाइन करो
- कंबाइन किए गए GVCF डेटा स्टोर पर जॉइंट जेनोटाइपिंग रन करके कोहोर्ट-लेवल VCF बनाओ
यह भाग सीधे भाग 2 द्वारा बनाए गए वर्कफ़्लो पर आधारित है।
इस सेक्शन से कैसे शुरू करें
कोर्स का यह सेक्शन मानता है कि तुमने भाग 2: Per-sample वेरिएंट कॉलिंग पूरा कर लिया है और तुम्हारे पास एक काम करती हुई genomics.nf पाइपलाइन है।
अगर तुमने भाग 2 पूरा नहीं किया है या इस भाग के लिए नए सिरे से शुरू करना चाहते हो, तो तुम भाग 2 के सॉल्यूशन को अपने शुरुआती बिंदु के रूप में उपयोग कर सकते हो।
nf4-science/genomics/ डायरेक्टरी के अंदर से ये कमांड्स रन करो:
cp solutions/part2/genomics-2.nf genomics.nf
cp solutions/part2/nextflow.config .
cp solutions/part2/modules/* modules/
इससे तुम्हें एक पूर्ण per-sample वेरिएंट कॉलिंग वर्कफ़्लो मिल जाएगा। तुम यह टेस्ट कर सकते हो कि यह सफलतापूर्वक रन हो रहा है, निम्नलिखित कमांड रन करके:
पाठ योजना¶
हमने इसे दो चरणों में विभाजित किया है:
- Per-sample वेरिएंट कॉलिंग स्टेप को संशोधित करके एक GVCF बनाओ। इसमें प्रोसेस कमांड्स और आउटपुट्स को अपडेट करना शामिल है।
- एक जॉइंट जेनोटाइपिंग स्टेप जोड़ो जो per-sample GVCFs को कंबाइन और जेनोटाइप करे।
इसमें
collect()ऑपरेटर, कमांड-लाइन कंस्ट्रक्शन के लिए Groovy closures, और मल्टी-कमांड प्रोसेसेस का परिचय होगा।
Note
सुनिश्चित करो कि तुम सही वर्किंग डायरेक्टरी में हो:
cd /workspaces/training/nf4-science/genomics
1. Per-sample वेरिएंट कॉलिंग स्टेप को GVCF बनाने के लिए संशोधित करो¶
भाग 2 की पाइपलाइन VCF फ़ाइलें बनाती है, लेकिन जॉइंट कॉलिंग के लिए GVCF फ़ाइलें चाहिए। हमें GVCF वेरिएंट कॉलिंग मोड को स्विच ऑन करना होगा और आउटपुट फ़ाइल एक्सटेंशन को अपडेट करना होगा।
भाग 1 से GVCF वेरिएंट कॉलिंग कमांड याद करो:
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
भाग 2 में हमने जिस बेस HaplotypeCaller कमांड को रैप किया था, उसकी तुलना में अंतर -ERC GVCF पैरामीटर और .g.vcf आउटपुट एक्सटेंशन हैं।
1.1. HaplotypeCaller को GVCF emit करने के लिए बताओ और आउटपुट एक्सटेंशन को अपडेट करो¶
modules/gatk_haplotypecaller.nf मॉड्यूल फ़ाइल खोलो और दो बदलाव करो:
- GATK HaplotypeCaller कमांड में
-ERC GVCFपैरामीटर जोड़ो; - आउटपुट फ़ाइल पाथ को अपडेट करके संबंधित
.g.vcfएक्सटेंशन का उपयोग करो, GATK कन्वेंशन के अनुसार।
सुनिश्चित करो कि जब तुम -ERC GVCF जोड़ो तो पिछली लाइन के अंत में बैकस्लैश (\) जोड़ो।
हमें आउटपुट ब्लॉक को भी नए फ़ाइल एक्सटेंशन से मैच करने के लिए अपडेट करना होगा।
चूंकि हमने कमांड आउटपुट को .vcf से .g.vcf में बदल दिया है, प्रोसेस output: ब्लॉक को भी वही बदलाव दर्शाना चाहिए।
1.2. प्रोसेस आउटपुट्स ब्लॉक में आउटपुट फ़ाइल एक्सटेंशन को अपडेट करो¶
हमें नए GVCF आउटपुट्स को दर्शाने के लिए वर्कफ़्लो के publish और output कॉन्फ़िगरेशन को भी अपडेट करना होगा।
1.3. नए GVCF आउटपुट्स के लिए publish टार्गेट्स को अपडेट करो¶
चूंकि अब हम VCFs के बजाय GVCFs बना रहे हैं, हमें वर्कफ़्लो के publish: सेक्शन को अधिक वर्णनात्मक नामों का उपयोग करने के लिए अपडेट करना चाहिए।
हम स्पष्टता के लिए GVCF फ़ाइलों को अपनी सबडायरेक्टरी में व्यवस्थित भी करेंगे।
अब आउटपुट ब्लॉक को मैच करने के लिए अपडेट करो।
1.4. नई डायरेक्टरी स्ट्रक्चर के लिए आउटपुट ब्लॉक को अपडेट करो¶
हमें GVCF फ़ाइलों को एक gvcf सबडायरेक्टरी में रखने के लिए output ब्लॉक को भी अपडेट करना होगा।
मॉड्यूल, publish टार्गेट्स, और आउटपुट ब्लॉक सभी अपडेट होने के साथ, हम बदलावों को टेस्ट कर सकते हैं।
1.5. पाइपलाइन को रन करो¶
बदलावों को वेरिफाई करने के लिए वर्कफ़्लो को रन करो।
कमांड आउटपुट
Nextflow आउटपुट पहले जैसा दिखता है, लेकिन .g.vcf फ़ाइलें और उनकी इंडेक्स फ़ाइलें अब सबडायरेक्टरीज़ में व्यवस्थित हैं।
डायरेक्टरी सामग्री (symlinks छोटे किए गए)
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
अगर तुम किसी एक GVCF फ़ाइल को खोलो और उसे स्क्रॉल करो, तो तुम वेरिफाई कर सकते हो कि GATK HaplotypeCaller ने अनुरोध के अनुसार GVCF फ़ाइलें बनाई हैं।
सारांश¶
जब तुम किसी टूल कमांड की आउटपुट फ़ाइलनेम बदलते हो, तो प्रोसेस output: ब्लॉक और publish/output कॉन्फ़िगरेशन को मैच करने के लिए अपडेट करना होगा।
आगे क्या है?¶
चैनल की सामग्री को इकट्ठा करना और उन्हें अगली प्रोसेस में एक सिंगल इनपुट के रूप में पास करना सीखो।
2. एक जॉइंट जेनोटाइपिंग स्टेप जोड़ो¶
अब हमें per-sample GVCFs को इकट्ठा करना होगा, उन्हें एक GenomicsDB डेटा स्टोर में कंबाइन करना होगा, और कोहोर्ट-लेवल VCF बनाने के लिए जॉइंट जेनोटाइपिंग रन करना होगा।
जैसा कि भाग 1 में बताया गया है, यह दो-टूल ऑपरेशन है: GenomicsDBImport GVCFs को कंबाइन करता है, फिर GenotypeGVCFs अंतिम वेरिएंट कॉल्स बनाता है।
हम दोनों टूल्स को एक सिंगल प्रोसेस में रैप करेंगे जिसे GATK_JOINTGENOTYPING कहा जाएगा।
भाग 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
पहली कमांड per-sample GVCFs और एक intervals फ़ाइल लेती है, और एक GenomicsDB डेटा स्टोर बनाती है।
दूसरी उस डेटा स्टोर, एक रेफरेंस जीनोम लेती है, और अंतिम कोहोर्ट-लेवल VCF बनाती है।
कंटेनर URI HaplotypeCaller के समान है: community.wave.seqera.io/library/gatk4:4.5.0.0--730ee8817e436867।
2.1. इनपुट्स सेट अप करो¶
जॉइंट जेनोटाइपिंग प्रोसेस को दो तरह के इनपुट्स चाहिए जो हमारे पास अभी नहीं हैं: एक arbitrary कोहोर्ट नाम, और सभी सैंपल्स से इकट्ठे किए गए GVCF आउटपुट्स एक साथ बंडल किए हुए।
2.1.1. एक cohort_name पैरामीटर जोड़ो¶
हमें कोहोर्ट के लिए एक arbitrary नाम प्रदान करना होगा।
ट्रेनिंग सीरीज़ में बाद में तुम सीखोगे कि इस तरह की चीज़ों के लिए सैंपल मेटाडेटा का उपयोग कैसे करें, लेकिन अभी हम सिर्फ params का उपयोग करके एक CLI पैरामीटर डिक्लेयर करते हैं और सुविधा के लिए इसे एक डिफ़ॉल्ट वैल्यू देते हैं।
2.1.2. सैंपल्स में HaplotypeCaller आउटपुट्स को इकट्ठा करो¶
अगर हम GATK_HAPLOTYPECALLER से आउटपुट चैनल को सीधे नई प्रोसेस में प्लग करें, तो Nextflow प्रोसेस को हर सैंपल GVCF पर अलग से कॉल करेगा।
हम सभी तीन GVCFs (और उनकी इंडेक्स फ़ाइलों) को बंडल करना चाहते हैं ताकि Nextflow उन सभी को एक साथ एक सिंगल प्रोसेस कॉल में दे।
हम collect() चैनल ऑपरेटर का उपयोग करके ऐसा कर सकते हैं।
GATK_HAPLOTYPECALLER की कॉल के ठीक बाद workflow बॉडी में निम्नलिखित लाइनें जोड़ो:
इसे तोड़ते हुए:
- हम
.outप्रॉपर्टी का उपयोग करकेGATK_HAPLOTYPECALLERसे आउटपुट चैनल लेते हैं। - क्योंकि हमने सेक्शन 1 में
emit:का उपयोग करके आउटपुट्स को नाम दिया था, हम.vcfसे GVCFs और.idxसे इंडेक्स फ़ाइलों को सिलेक्ट कर सकते हैं। नामित आउटपुट्स के बिना, हमें.out[0]और.out[1]का उपयोग करना पड़ता। collect()ऑपरेटर सभी फ़ाइलों को एक सिंगल एलिमेंट में बंडल करता है, इसलिएall_gvcfs_chमें सभी तीन GVCFs एक साथ होते हैं, औरall_idxs_chमें सभी तीन इंडेक्स फ़ाइलें एक साथ होती हैं।
हम GVCFs और उनकी इंडेक्स फ़ाइलों को अलग से इकट्ठा कर सकते हैं (tuples में उन्हें एक साथ रखने के विपरीत) क्योंकि Nextflow execution के लिए सभी इनपुट फ़ाइलों को एक साथ stage करेगा, इसलिए इंडेक्स फ़ाइलें GVCFs के साथ मौजूद होंगी।
Tip
चैनल ऑपरेटर्स लागू करने से पहले और बाद में चैनलों की सामग्री को इंस्पेक्ट करने के लिए तुम view() ऑपरेटर का उपयोग कर सकते हो।
2.2. जॉइंट जेनोटाइपिंग प्रोसेस लिखो और इसे वर्कफ़्लो में कॉल करो¶
भाग 2 में उपयोग किए गए समान पैटर्न का पालन करते हुए, हम मॉड्यूल फ़ाइल में प्रोसेस डेफिनिशन लिखेंगे, इसे वर्कफ़्लो में इंपोर्ट करेंगे, और इसे जिन इनपुट्स को हमने अभी तैयार किया है, उन पर कॉल करेंगे।
2.2.1. हर GVCF को -V आर्ग्युमेंट देने के लिए एक स्ट्रिंग कंस्ट्रक्ट करो¶
प्रोसेस डेफिनिशन को भरना शुरू करने से पहले, एक चीज़ काम करनी है।
GenomicsDBImport कमांड हर GVCF फ़ाइल के लिए एक अलग -V आर्ग्युमेंट की उम्मीद करता है, इस तरह:
gatk GenomicsDBImport \
-V reads_mother.bam.g.vcf \
-V reads_father.bam.g.vcf \
-V reads_son.bam.g.vcf \
...
अगर हम -V ${all_gvcfs_ch} लिखें, तो Nextflow बस फ़ाइलनेम्स को concatenate कर देगा और कमांड का वह हिस्सा इस तरह दिखेगा:
लेकिन हमें स्ट्रिंग को इस तरह दिखना चाहिए:
महत्वपूर्ण रूप से, हमें इस स्ट्रिंग को collected चैनल में जो भी फ़ाइलें हैं, उनसे dynamically कंस्ट्रक्ट करनी होगी। Nextflow (Groovy के माध्यम से) ऐसा करने का एक संक्षिप्त तरीका प्रदान करता है:
इसे तोड़ते हुए:
all_gvcfs.collect { gvcf -> "-V ${gvcf}" }हर फ़ाइल पाथ पर iterate करता है और उसमें-Vprepend करता है,["-V A.g.vcf", "-V B.g.vcf", "-V C.g.vcf"]बनाता है।.join(' ')उन्हें स्पेस के साथ concatenate करता है:"-V A.g.vcf -V B.g.vcf -V C.g.vcf"।- परिणाम एक लोकल वेरिएबल
gvcfs_lineको असाइन किया जाता है (defके साथ डिफाइन किया गया), जिसे हम कमांड टेम्पलेट में interpolate कर सकते हैं।
यह लाइन प्रोसेस के script: ब्लॉक के अंदर जाती है, कमांड टेम्पलेट से पहले।
तुम कमांड टेम्पलेट के opening """ से पहले script: के बीच arbitrary Groovy कोड रख सकते हो।
फिर तुम प्रोसेस के script: ब्लॉक में उस पूरी स्ट्रिंग को gvcfs_line के रूप में refer कर पाओगे।
2.2.2. जॉइंट जेनोटाइपिंग प्रोसेस के लिए मॉड्यूल को भरो¶
अब हम पूर्ण प्रोसेस लिखने की ओर बढ़ सकते हैं।
modules/gatk_jointgenotyping.nf खोलो और प्रोसेस डेफिनिशन की outline को एग्ज़ामिन करो।
ऊपर दी गई जानकारी का उपयोग करके प्रोसेस डेफिनिशन को भरो, फिर नीचे "बाद में" टैब में सॉल्यूशन के खिलाफ अपना काम चेक करो।
यहाँ कई चीज़ें उल्लेखनीय हैं।
पहले की तरह, कई इनपुट्स लिस्ट किए गए हैं भले ही कमांड्स उन्हें सीधे रेफरेंस नहीं करते: all_idxs, ref_index, और ref_dict।
उन्हें लिस्ट करना सुनिश्चित करता है कि Nextflow इन फ़ाइलों को वर्किंग डायरेक्टरी में उन फ़ाइलों के साथ stage करे जो कमांड्स में दिखाई देती हैं, जिन्हें GATK naming conventions के आधार पर खोजने की उम्मीद करता है।
gvcfs_line वेरिएबल GenomicsDBImport के लिए -V आर्ग्युमेंट्स कंस्ट्रक्ट करने के लिए ऊपर वर्णित Groovy closure का उपयोग करता है।
यह प्रोसेस serial में दो कमांड्स रन करता है, ठीक जैसे तुम टर्मिनल में करते। GenomicsDBImport per-sample GVCFs को एक डेटा स्टोर में कंबाइन करता है, फिर GenotypeGVCFs उस डेटा स्टोर को पढ़ता है और अंतिम कोहोर्ट-लेवल VCF बनाता है।
GenomicsDB डेटा स्टोर (${cohort_name}_gdb) एक intermediate artifact है जिसका उपयोग केवल प्रोसेस के अंदर किया जाता है; यह आउटपुट ब्लॉक में दिखाई नहीं देता।
एक बार जब तुम यह पूरा कर लो, तो प्रोसेस उपयोग के लिए तैयार है। वर्कफ़्लो में इसका उपयोग करने के लिए, तुम्हें मॉड्यूल को इंपोर्ट करना होगा और एक प्रोसेस कॉल जोड़नी होगी।
2.2.3. मॉड्यूल को इंपोर्ट करो¶
मौजूदा इंपोर्ट स्टेटमेंट्स के नीचे genomics.nf में इंपोर्ट स्टेटमेंट जोड़ो:
प्रोसेस अब वर्कफ़्लो scope में उपलब्ध है।
2.2.4. प्रोसेस कॉल जोड़ो¶
collect() लाइनों के बाद वर्कफ़्लो बॉडी में GATK_JOINTGENOTYPING की कॉल जोड़ो:
प्रोसेस अब पूरी तरह से wired up है। अगला, हम configure करेंगे कि आउटपुट्स कैसे publish होते हैं।
2.3. आउटपुट हैंडलिंग को कॉन्फ़िगर करो¶
हमें जॉइंट VCF आउटपुट्स को publish करना होगा। जॉइंट जेनोटाइपिंग परिणामों के लिए publish टार्गेट्स और आउटपुट ब्लॉक एंट्रीज़ जोड़ो।
2.3.1. जॉइंट VCF के लिए publish टार्गेट्स जोड़ो¶
वर्कफ़्लो के publish: सेक्शन में जॉइंट VCF और उसका इंडेक्स जोड़ो:
अब मैच करने के लिए आउटपुट ब्लॉक को अपडेट करो।
2.3.2. जॉइंट VCF के लिए आउटपुट ब्लॉक एंट्रीज़ जोड़ो¶
जॉइंट VCF फ़ाइलों के लिए एंट्रीज़ जोड़ो। हम उन्हें results डायरेक्टरी की रूट पर रखेंगे क्योंकि यह अंतिम आउटपुट है।
प्रोसेस, publish टार्गेट्स, और आउटपुट ब्लॉक सभी के साथ, हम पूर्ण वर्कफ़्लो को टेस्ट कर सकते हैं।
2.4. वर्कफ़्लो को रन करो¶
सब कुछ काम करता है वेरिफाई करने के लिए वर्कफ़्लो को रन करो।
कमांड आउटपुट
पहले दो स्टेप्स पिछले रन से cached हैं, और नया GATK_JOINTGENOTYPING स्टेप सभी तीन सैंपल्स से collected इनपुट्स पर एक बार रन होता है।
अंतिम आउटपुट फ़ाइल, family_trio.joint.vcf (और उसका इंडेक्स), results डायरेक्टरी में है।
डायरेक्टरी सामग्री (symlinks छोटे किए गए)
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
अगर तुम जॉइंट VCF फ़ाइल खोलो, तो तुम वेरिफाई कर सकते हो कि वर्कफ़्लो ने अपेक्षित वेरिएंट कॉल्स बनाए।
अब तुम्हारे पास एक automated, पूर्ण रूप से reproducible जॉइंट वेरिएंट कॉलिंग वर्कफ़्लो है!
Note
ध्यान रखो कि हमने तुम्हें जो डेटा फ़ाइलें दी हैं वे chromosome 20 के केवल एक छोटे से हिस्से को कवर करती हैं। वेरिएंट कॉलसेट का वास्तविक साइज़ लाखों वेरिएंट्स में गिना जाएगा। इसीलिए हम ट्रेनिंग उद्देश्यों के लिए केवल छोटे डेटा सबसेट्स का उपयोग करते हैं!
सारांश¶
तुम जानते हो कि चैनल से आउटपुट्स कैसे इकट्ठा करें और उन्हें किसी अन्य प्रोसेस के लिए एक सिंगल इनपुट के रूप में बंडल करें। तुम Groovy closures का उपयोग करके कमांड लाइन कैसे कंस्ट्रक्ट करें, और एक सिंगल प्रोसेस में मल्टिपल कमांड्स कैसे रन करें, यह भी जानते हो।
आगे क्या है?¶
अपनी सफलता का जश्न मनाओ और एक अच्छा ब्रेक लो।
इस कोर्स के अगले भाग में, तुम सीखोगे कि nf-core से एक production-ready वेरिएंट कॉलिंग पाइपलाइन कैसे रन करें और इसकी तुलना उस पाइपलाइन से कैसे करें जो तुमने मैन्युअली बनाई है।