भाग 1: विधि का अवलोकन¶
AI-सहायता प्राप्त अनुवाद - अधिक जानें और सुधार सुझाएं
वेरिएंट कॉलिंग एक जीनोमिक विश्लेषण विधि है जिसका उद्देश्य संदर्भ जीनोम की तुलना में जीनोम अनुक्रम में भिन्नता की पहचान करना है। यहां हम होल-जीनोम सीक्वेंसिंग डेटा में शॉर्ट जर्मलाइन वेरिएंट, यानी SNPs और indels, को कॉल करने के लिए डिज़ाइन किए गए टूल्स और विधियों का उपयोग करने जा रहे हैं।

एक पूर्ण वेरिएंट कॉलिंग पाइपलाइन में आमतौर पर कई चरण शामिल होते हैं, जिनमें रेफरेंस से मैपिंग (कभी-कभी जीनोम संरेखण के रूप में संदर्भित) और वेरिएंट फ़िल्टरिंग और प्राथमिकता शामिल हैं। सरलता के लिए, इस कोर्स में हम सिर्फ वेरिएंट कॉलिंग भाग पर ध्यान केंद्रित करने जा रहे हैं।
विधियां¶
हम तुम्हें होल-जीनोम सीक्वेंसिंग नमूनों पर वेरिएंट कॉलिंग लागू करने के दो तरीके दिखाने जा रहे हैं ताकि जर्मलाइन SNPs और indels की पहचान की जा सके। पहले हम एक सरल प्रति-नमूना दृष्टिकोण से शुरुआत करेंगे जो प्रत्येक नमूने से स्वतंत्र रूप से वेरिएंट को कॉल करता है। फिर हम तुम्हें एक अधिक परिष्कृत संयुक्त कॉलिंग दृष्टिकोण दिखाएंगे जो कई नमूनों का एक साथ विश्लेषण करता है, और अधिक सटीक और जानकारीपूर्ण परिणाम उत्पन्न करता है।
इससे पहले कि हम किसी भी दृष्टिकोण के लिए कोई वर्कफ़्लो कोड लिखना शुरू करें, हम कुछ परीक्षण डेटा पर मैनुअल रूप से कमांड्स को आज़माने जा रहे हैं।
डेटासेट¶
हम निम्नलिखित डेटा और संबंधित संसाधन प्रदान करते हैं:
- एक रेफरेंस जीनोम जो मानव क्रोमोसोम 20 (hg19/b37 से) के एक छोटे से क्षेत्र और इसकी सहायक फ़ाइलों (इंडेक्स और अनुक्रम शब्दकोश) से बना है।
- तीन होल जीनोम सीक्वेंसिंग नमूने जो एक परिवार ट्रायो (मां, पिता और बेटा) से संबंधित हैं, जिन्हें क्रोमोसोम 20 पर डेटा के एक छोटे से हिस्से तक सीमित कर दिया गया है ताकि फ़ाइल का आकार छोटा रहे। यह Illumina शॉर्ट-रीड सीक्वेंसिंग डेटा है जिसे पहले से ही रेफरेंस जीनोम से मैप किया जा चुका है, BAM फॉर्मेट में प्रदान किया गया है (Binary Alignment Map, SAM का संकुचित संस्करण, Sequence Alignment Map)।
- जीनोमिक अंतराल की एक सूची, यानी जीनोम पर निर्देशांक जहां हमारे नमूनों में वेरिएंट कॉल करने के लिए उपयुक्त डेटा है, BED फॉर्मेट में प्रदान किया गया है।
सॉफ्टवेयर¶
शामिल दो मुख्य टूल हैं Samtools, अनुक्रम संरेखण फ़ाइलों में हेरफेर करने के लिए व्यापक रूप से उपयोग किया जाने वाला टूलकिट, और GATK (Genome Analysis Toolkit), Broad Institute में विकसित वेरिएंट खोज के लिए टूल्स का एक सेट।
ये टूल्स GitHub Codespaces वातावरण में इंस्टॉल नहीं हैं, इसलिए हम उन्हें Seqera Containers सेवा के माध्यम से प्राप्त कंटेनरों के माध्यम से उपयोग करेंगे (Hello Containers देखें)।
सुझाव
सुनिश्चित करो कि तुम nf4-science/genomics डायरेक्टरी में हो ताकि जब तुम pwd टाइप करो तो पथ का अंतिम भाग genomics दिखे।
1. प्रति-नमूना वेरिएंट कॉलिंग¶
प्रति-नमूना वेरिएंट कॉलिंग प्रत्येक नमूने को स्वतंत्र रूप से प्रोसेस करती है: वेरिएंट कॉलर एक समय में एक नमूने के लिए सीक्वेंसिंग डेटा की जांच करता है और उन पोजीशन्स की पहचान करता है जहां नमूना रेफरेंस से अलग है।
इस सेक्शन में हम प्रति-नमूना वेरिएंट कॉलिंग दृष्टिकोण के दो कमांड्स का परीक्षण करते हैं: Samtools के साथ BAM फ़ाइल को इंडेक्सिंग करना और GATK HaplotypeCaller के साथ वेरिएंट कॉल करना। ये वे कमांड्स हैं जिन्हें हम इस कोर्स के भाग 2 में Nextflow वर्कफ़्लो में wrap करेंगे।
- Samtools का उपयोग करके BAM इनपुट फ़ाइल के लिए एक इंडेक्स फ़ाइल जेनरेट करें
- इंडेक्स की गई BAM फ़ाइल पर GATK HaplotypeCaller चलाएं ताकि VCF (Variant Call Format) में प्रति-नमूना वेरिएंट कॉल जेनरेट हो सकें
हम सिर्फ एक नमूने पर दो कमांड्स का परीक्षण करके शुरुआत करते हैं।
1.1. Samtools के साथ BAM इनपुट फ़ाइल को इंडेक्स करें¶
इंडेक्स फ़ाइलें बायोइन्फॉर्मेटिक्स फ़ाइल फॉर्मेट की एक सामान्य विशेषता हैं; इनमें मुख्य फ़ाइल की संरचना के बारे में जानकारी होती है जो GATK जैसे टूल्स को पूरी फ़ाइल पढ़े बिना डेटा के एक सबसेट तक पहुंचने की अनुमति देती है। यह महत्वपूर्ण है क्योंकि ये फ़ाइलें कितनी बड़ी हो सकती हैं।
BAM फ़ाइलें अक्सर इंडेक्स के बिना प्रदान की जाती हैं, इसलिए कई विश्लेषण वर्कफ़्लो में पहला चरण samtools index का उपयोग करके एक जेनरेट करना है।
हम एक Samtools कंटेनर को पुल करने जा रहे हैं, इसे इंटरएक्टिव रूप से स्पिन अप करेंगे और BAM फ़ाइलों में से एक पर samtools index कमांड चलाएंगे।
1.1.1. Samtools कंटेनर को पुल करें¶
Samtools कंटेनर इमेज डाउनलोड करने के लिए docker pull कमांड चलाएं:
कमांड आउटपुट
1.20--b5dfbd93de237464: Pulling from library/samtools
6360b3717211: Pull complete
2ec3f7ad9b3c: Pull complete
7716ca300600: Pull complete
4f4fb700ef54: Pull complete
8c61d418774c: Pull complete
03dae77ff45c: Pull complete
aab7f787139d: Pull complete
4f4fb700ef54: Pull complete
837d55536720: Pull complete
897362c12ca7: Pull complete
3893cbe24e91: Pull complete
d1b61e94977b: Pull complete
c72ff66fb90f: Pull complete
0e0388f29b6d: Pull complete
Digest: sha256:bbfc45b4f228975bde86cba95e303dd94ecf2fdacea5bfb2e2f34b0d7b141e41
Status: Downloaded newer image for community.wave.seqera.io/library/samtools:1.20--b5dfbd93de237464
community.wave.seqera.io/library/samtools:1.20--b5dfbd93de237464
अगर तुमने यह इमेज पहले डाउनलोड नहीं की है, तो इसे पूरा होने में एक मिनट लग सकता है। एक बार यह हो जाने के बाद, तुम्हारे पास कंटेनर इमेज की एक लोकल कॉपी है।
1.1.2. Samtools कंटेनर को इंटरएक्टिव रूप से स्पिन अप करें¶
कंटेनर को इंटरएक्टिव रूप से चलाने के लिए, -it फ्लैग के साथ docker run का उपयोग करो।
-v ./data:/data विकल्प लोकल data डायरेक्टरी को कंटेनर में माउंट करता है ताकि टूल्स इनपुट फ़ाइलों तक पहुंच सकें।
तुम्हारा प्रॉम्प्ट (base) root@a1b2c3d4e5f6:/tmp# जैसे कुछ में बदल जाता है, जो दर्शाता है कि तुम अब कंटेनर के अंदर हो।
सत्यापित करो कि तुम /data/bam के तहत सीक्वेंस डेटा फ़ाइलें देख सकते हो:
इसके साथ, तुम अपनी पहली कमांड आज़माने के लिए तैयार हो।
1.1.3. इंडेक्सिंग कमांड चलाएं¶
Samtools डॉक्यूमेंटेशन हमें BAM फ़ाइल को इंडेक्स करने के लिए चलाने वाली कमांड लाइन देता है।
हमें केवल इनपुट फ़ाइल प्रदान करने की आवश्यकता है; टूल स्वचालित रूप से इनपुट फ़ाइलनाम में .bai जोड़कर आउटपुट के लिए एक नाम जेनरेट करेगा।
एक डेटा फ़ाइल पर samtools index कमांड चलाएं:
कमांड टर्मिनल में कोई आउटपुट उत्पन्न नहीं करती है, लेकिन तुम्हें अब मूल BAM इनपुट फ़ाइल के समान डायरेक्टरी में reads_mother.bam.bai नामक एक फ़ाइल दिखनी चाहिए।
डायरेक्टरी सामग्री
यह पहले चरण के परीक्षण को पूरा करता है।
1.1.4. Samtools कंटेनर से बाहर निकलें¶
कंटेनर से बाहर निकलने के लिए, exit टाइप करो।
तुम्हारा प्रॉम्प्ट अब वापस वैसा होना चाहिए जैसा कंटेनर शुरू करने से पहले था।
1.2. GATK HaplotypeCaller के साथ वेरिएंट कॉल करें¶
हम BAM फ़ाइल पर gatk HaplotypeCaller कमांड चलाना चाहते हैं जिसे हमने अभी इंडेक्स किया है।
1.2.1. GATK कंटेनर को पुल करें¶
पहले, आओ GATK कंटेनर इमेज डाउनलोड करने के लिए docker pull कमांड चलाएं:
कमांड आउटपुट
कुछ लेयर्स Already exists दिखाती हैं क्योंकि वे Samtools कंटेनर इमेज के साथ शेयर की गई हैं जिसे हमने पहले पुल किया था।
4.5.0.0--730ee8817e436867: Pulling from library/gatk4
6360b3717211: Already exists
2ec3f7ad9b3c: Already exists
7716ca300600: Already exists
4f4fb700ef54: Already exists
8c61d418774c: Already exists
03dae77ff45c: Already exists
aab7f787139d: Already exists
4f4fb700ef54: Already exists
837d55536720: Already exists
897362c12ca7: Already exists
3893cbe24e91: Already exists
d1b61e94977b: Already exists
e5c558f54708: Pull complete
087cce32d294: Pull complete
Digest: sha256:e33413b9100f834fcc62fd5bc9edc1e881e820aafa606e09301eac2303d8724b
Status: Downloaded newer image for community.wave.seqera.io/library/gatk4:4.5.0.0--730ee8817e436867
community.wave.seqera.io/library/gatk4:4.5.0.0--730ee8817e436867
यह पहले पुल की तुलना में तेज़ होना चाहिए क्योंकि दोनों कंटेनर इमेजेस अपनी अधिकांश लेयर्स शेयर करती हैं।
1.2.2. GATK कंटेनर को इंटरएक्टिव रूप से स्पिन अप करें¶
GATK कंटेनर को डेटा डायरेक्टरी माउंट के साथ इंटरएक्टिव रूप से स्पिन अप करो, जैसे हमने Samtools के लिए किया था।
तुम्हारा प्रॉम्प्ट बदल जाता है जो दर्शाता है कि तुम अब GATK कंटेनर के अंदर हो।
1.2.3. वेरिएंट कॉलिंग कमांड चलाएं¶
GATK डॉक्यूमेंटेशन हमें BAM फ़ाइल पर वेरिएंट कॉलिंग करने के लिए चलाने वाली कमांड लाइन देता है।
हमें BAM इनपुट फ़ाइल (-I) के साथ-साथ रेफरेंस जीनोम (-R), आउटपुट फ़ाइल के लिए एक नाम (-O) और विश्लेषण करने के लिए जीनोमिक अंतराल की एक सूची (-L) प्रदान करने की आवश्यकता है।
हालांकि, हमें इंडेक्स फ़ाइल का पथ निर्दिष्ट करने की आवश्यकता नहीं है; टूल स्वचालित रूप से स्थापित नामकरण और सह-स्थान सम्मेलन के आधार पर इसे उसी डायरेक्टरी में देखेगा।
यही रेफरेंस जीनोम की सहायक फ़ाइलों (इंडेक्स और अनुक्रम शब्दकोश फ़ाइलें, *.fai और *.dict) पर भी लागू होता है।
gatk HaplotypeCaller \
-R /data/ref/ref.fasta \
-I /data/bam/reads_mother.bam \
-O /data/vcf/reads_mother.vcf \
-L /data/ref/intervals.bed
कमांड आउटपुट
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar HaplotypeCaller -R /data/ref/ref.fasta -I /data/bam/reads_mother.bam -O reads_mother.vcf -L /data/ref/intervals.bed
00:27:50.687 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
00:27:50.854 INFO HaplotypeCaller - ------------------------------------------------------------
00:27:50.858 INFO HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.5.0.0
00:27:50.858 INFO HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
00:27:50.858 INFO HaplotypeCaller - Executing as root@a1fe8ff42d07 on Linux v6.10.14-linuxkit amd64
00:27:50.858 INFO HaplotypeCaller - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
00:27:50.859 INFO HaplotypeCaller - Start Date/Time: February 8, 2026 at 12:27:50 AM GMT
00:27:50.859 INFO HaplotypeCaller - ------------------------------------------------------------
00:27:50.859 INFO HaplotypeCaller - ------------------------------------------------------------
00:27:50.861 INFO HaplotypeCaller - HTSJDK Version: 4.1.0
00:27:50.861 INFO HaplotypeCaller - Picard Version: 3.1.1
00:27:50.861 INFO HaplotypeCaller - Built for Spark Version: 3.5.0
00:27:50.862 INFO HaplotypeCaller - HTSJDK Defaults.COMPRESSION_LEVEL : 2
00:27:50.862 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
00:27:50.862 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
00:27:50.863 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
00:27:50.864 INFO HaplotypeCaller - Deflater: IntelDeflater
00:27:50.864 INFO HaplotypeCaller - Inflater: IntelInflater
00:27:50.864 INFO HaplotypeCaller - GCS max retries/reopens: 20
00:27:50.864 INFO HaplotypeCaller - Requester pays: disabled
00:27:50.865 INFO HaplotypeCaller - Initializing engine
00:27:50.991 INFO FeatureManager - Using codec BEDCodec to read file file:///data/ref/intervals.bed
00:27:51.016 INFO IntervalArgumentCollection - Processing 6369 bp from intervals
00:27:51.029 INFO HaplotypeCaller - Done initializing engine
00:27:51.040 INFO NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_utils.so
00:27:51.042 INFO NativeLibraryLoader - Loading libgkl_smithwaterman.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_smithwaterman.so
00:27:51.042 INFO SmithWatermanAligner - Using AVX accelerated SmithWaterman implementation
00:27:51.046 INFO HaplotypeCallerEngine - Disabling physical phasing, which is supported only for reference-model confidence output
00:27:51.063 INFO NativeLibraryLoader - Loading libgkl_pairhmm_omp.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_pairhmm_omp.so
00:27:51.085 INFO IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM
00:27:51.086 INFO IntelPairHmm - Available threads: 10
00:27:51.086 INFO IntelPairHmm - Requested threads: 4
00:27:51.086 INFO PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
00:27:51.128 INFO ProgressMeter - Starting traversal
00:27:51.136 INFO ProgressMeter - Current Locus Elapsed Minutes Regions Processed Regions/Minute
00:27:51.882 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position 20_10037292_10066351:3480 and possibly subsequent; at least 10 samples must have called genotypes
00:27:52.969 INFO HaplotypeCaller - 7 read(s) filtered by: MappingQualityReadFilter
0 read(s) filtered by: MappingQualityAvailableReadFilter
0 read(s) filtered by: MappedReadFilter
0 read(s) filtered by: NotSecondaryAlignmentReadFilter
0 read(s) filtered by: NotDuplicateReadFilter
0 read(s) filtered by: PassesVendorQualityCheckReadFilter
0 read(s) filtered by: NonZeroReferenceLengthAlignmentReadFilter
0 read(s) filtered by: GoodCigarReadFilter
0 read(s) filtered by: WellformedReadFilter
7 total reads filtered out of 1867 reads processed
00:27:52.971 INFO ProgressMeter - 20_10037292_10066351:13499 0.0 35 1145.7
00:27:52.971 INFO ProgressMeter - Traversal complete. Processed 35 total regions in 0.0 minutes.
00:27:52.976 INFO VectorLoglessPairHMM - Time spent in setup for JNI call : 0.003346916
00:27:52.976 INFO PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.045731709
00:27:52.977 INFO SmithWatermanAligner - Total compute time in native Smith-Waterman : 0.02 sec
00:27:52.981 INFO HaplotypeCaller - Shutting down engine
[February 8, 2026 at 12:27:52 AM GMT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.04 minutes.
Runtime.totalMemory()=203423744
लॉग आउटपुट बहुत verbose है, इसलिए हमने ऊपर दिए गए उदाहरण में सबसे प्रासंगिक लाइनों को हाइलाइट किया है।
आउटपुट फ़ाइलें, reads_mother.vcf और इसकी इंडेक्स फ़ाइल, reads_mother.vcf.idx, कंटेनर में तुम्हारी कार्यशील डायरेक्टरी के अंदर बनाई गई हैं।
VCF फ़ाइल में वेरिएंट कॉल होते हैं, जैसा कि हम एक मिनट में देखेंगे, और इंडेक्स फ़ाइल का वही कार्य है जो BAM इंडेक्स फ़ाइल का है, टूल्स को संपूर्ण फ़ाइल लोड किए बिना डेटा के सबसेट को खोजने और पुनः प्राप्त करने की अनुमति देना।
चूंकि VCF एक टेक्स्ट फॉर्मेट है और यह एक छोटी परीक्षण फ़ाइल है, तुम इसे खोलने और सामग्री देखने के लिए cat reads_mother.vcf चला सकते हो।
अगर तुम फ़ाइल की शुरुआत तक स्क्रॉल करते हो, तो तुम्हें मेटाडेटा की कई लाइनों से बना एक हेडर मिलेगा, उसके बाद वेरिएंट कॉल की एक सूची होगी, प्रति लाइन एक।
फ़ाइल की सामग्री (संक्षिप्त)
ऊपर दिए गए उदाहरण आउटपुट में, हमने अंतिम हेडर लाइन को हाइलाइट किया है, जो उसके बाद आने वाले टेबुलर डेटा के लिए कॉलम के नाम देती है। डेटा की प्रत्येक लाइन नमूने के सीक्वेंसिंग डेटा में पहचाने गए संभावित वेरिएंट का वर्णन करती है। VCF फॉर्मेट की व्याख्या करने के लिए मार्गदर्शन के लिए, यह उपयोगी लेख देखें।
1.2.4. आउटपुट फ़ाइलों को स्थानांतरित करें¶
कंटेनर के अंदर जो कुछ भी रहता है वह भविष्य के काम के लिए दुर्गम होगा।
BAM इंडेक्स फ़ाइल सीधे माउंट की गई फ़ाइलसिस्टम पर /data/bam डायरेक्टरी में बनाई गई थी, लेकिन VCF फ़ाइल और इसकी इंडेक्स नहीं, इसलिए हमें उन दोनों को मैनुअल रूप से स्थानांतरित करने की आवश्यकता है।
डायरेक्टरी सामग्री
एक बार यह हो जाने के बाद, सभी फ़ाइलें अब तुम्हारी सामान्य फ़ाइलसिस्टम में सुलभ हैं।
1.2.5. GATK कंटेनर से बाहर निकलें¶
कंटेनर से बाहर निकलने के लिए, exit टाइप करो।
तुम्हारा प्रॉम्प्ट वापस सामान्य होना चाहिए। यह प्रति-नमूना वेरिएंट कॉलिंग परीक्षण समाप्त करता है।
इसे वर्कफ़्लो के रूप में लिखो!
अगर तुम इस विश्लेषण को Nextflow वर्कफ़्लो के रूप में लागू करना शुरू करना चाहते हो तो तुम तुरंत भाग 2 पर जा सकते हो। तुम्हें बस भाग 3 पर जाने से पहले परीक्षण के दूसरे दौर को पूरा करने के लिए वापस आना होगा।
2. कोहॉर्ट पर संयुक्त कॉलिंग¶
जिस वेरिएंट कॉलिंग दृष्टिकोण का हमने अभी उपयोग किया, वह प्रति नमूना वेरिएंट कॉल जेनरेट करता है। यह प्रत्येक नमूने से अलग-अलग वेरिएंट को देखने के लिए ठीक है, लेकिन यह सीमित जानकारी देता है। यह अक्सर अधिक दिलचस्प होता है कि वेरिएंट कॉल कई नमूनों में कैसे भिन्न होते हैं। GATK इस उद्देश्य के लिए संयुक्त वेरिएंट कॉलिंग नामक एक वैकल्पिक विधि प्रदान करता है।
संयुक्त वेरिएंट कॉलिंग में प्रत्येक नमूने के लिए GVCF (Genomic VCF के लिए) नामक एक विशेष प्रकार का वेरिएंट आउटपुट जेनरेट करना शामिल है, फिर सभी नमूनों से GVCF डेटा को जोड़ना और एक 'संयुक्त जीनोटाइपिंग' सांख्यिकीय विश्लेषण चलाना।

एक नमूने के GVCF में जो विशेष है वह यह है कि इसमें जीनोम के लक्षित क्षेत्र में सभी पोजीशन्स के बारे में अनुक्रम डेटा सांख्यिकी को सारांशित करने वाले रिकॉर्ड होते हैं, न कि केवल उन पोजीशन्स के जहां प्रोग्राम को भिन्नता के प्रमाण मिले। यह संयुक्त जीनोटाइपिंग गणना के लिए महत्वपूर्ण है (आगे पढ़ना)।
GVCF को GATK HaplotypeCaller द्वारा उत्पादित किया जाता है, वही टूल जिसे हमने अभी परीक्षण किया, एक अतिरिक्त पैरामीटर (-ERC GVCF) के साथ।
GVCFs को जोड़ना GATK GenomicsDBImport के साथ किया जाता है, जो प्रति-नमूना कॉल को GenomicsDB डेटा स्टोर में जोड़ता है।
GenomicsDB डेटा स्टोर एक प्रकार का डेटाबेस फॉर्मेट है जो वेरिएंट जानकारी के लिए एक मध्यवर्ती भंडारण के रूप में कार्य करता है।
वास्तविक 'संयुक्त जीनोटाइपिंग' विश्लेषण फिर GATK GenotypeGVCFs के साथ किया जाता है।
यहां हम GVCFs जेनरेट करने और संयुक्त जीनोटाइपिंग चलाने के लिए आवश्यक कमांड्स का परीक्षण करते हैं। ये वे कमांड्स हैं जिन्हें हम इस कोर्स के भाग 3 में Nextflow वर्कफ़्लो में wrap करेंगे।
- Samtools का उपयोग करके प्रत्येक BAM इनपुट फ़ाइल के लिए एक इंडेक्स फ़ाइल जेनरेट करें
- प्रति-नमूना जीनोमिक वेरिएंट कॉल का GVCF जेनरेट करने के लिए प्रत्येक BAM इनपुट फ़ाइल पर GATK HaplotypeCaller चलाएं
- सभी GVCFs को एकत्र करें और उन्हें GenomicsDB डेटा स्टोर में जोड़ें
- कोहॉर्ट-स्तर VCF उत्पन्न करने के लिए संयुक्त GVCF डेटा स्टोर पर संयुक्त जीनोटाइपिंग चलाएं
हमें अब इन सभी कमांड्स का परीक्षण करने की आवश्यकता है, सभी तीन BAM फ़ाइलों को इंडेक्सिंग करने से शुरुआत करते हुए।
2.1. सभी तीन नमूनों के लिए BAM फ़ाइलों को इंडेक्स करें¶
ऊपर पहले सेक्शन में, हमने केवल एक BAM फ़ाइल को इंडेक्स किया था। अब हमें सभी तीन नमूनों को इंडेक्स करने की आवश्यकता है ताकि GATK HaplotypeCaller उन्हें प्रोसेस कर सके।
2.1.1. Samtools कंटेनर को इंटरएक्टिव रूप से स्पिन अप करें¶
हमने पहले ही Samtools कंटेनर इमेज को पुल कर लिया है, इसलिए हम इसे सीधे स्पिन अप कर सकते हैं:
तुम्हारा प्रॉम्प्ट बदलता है जो दर्शाता है कि तुम कंटेनर के अंदर हो, डेटा डायरेक्टरी पहले की तरह माउंट है।
2.1.2. सभी तीन नमूनों पर इंडेक्सिंग कमांड चलाएं¶
तीनों BAM फ़ाइलों में से प्रत्येक पर इंडेक्सिंग कमांड चलाएं:
samtools index /data/bam/reads_mother.bam
samtools index /data/bam/reads_father.bam
samtools index /data/bam/reads_son.bam
डायरेक्टरी सामग्री
यह संबंधित BAM फ़ाइलों के समान डायरेक्टरी में इंडेक्स फ़ाइलें उत्पन्न करना चाहिए।
2.1.3. Samtools कंटेनर से बाहर निकलें¶
कंटेनर से बाहर निकलने के लिए, exit टाइप करो।
तुम्हारा प्रॉम्प्ट वापस सामान्य होना चाहिए।
2.2. सभी तीन नमूनों के लिए GVCFs जेनरेट करें¶
संयुक्त जीनोटाइपिंग चरण चलाने के लिए, हमें सभी तीन नमूनों के लिए GVCFs की आवश्यकता है।
2.2.1. GATK कंटेनर को इंटरएक्टिव रूप से स्पिन अप करें¶
हमने पहले ही GATK कंटेनर इमेज को पुल कर लिया है, इसलिए हम इसे सीधे स्पिन अप कर सकते हैं:
तुम्हारा प्रॉम्प्ट बदलता है जो दर्शाता है कि तुम GATK कंटेनर के अंदर हो।
2.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
कमांड आउटपुट
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar HaplotypeCaller -R /data/ref/ref.fasta -I /data/bam/reads_mother.bam -O reads_mother.g.vcf -L /data/ref/intervals.bed -ERC GVCF
16:51:00.620 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
16:51:00.749 INFO HaplotypeCaller - ------------------------------------------------------------
16:51:00.751 INFO HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.5.0.0
16:51:00.751 INFO HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
16:51:00.751 INFO HaplotypeCaller - Executing as root@be1a0302f6c7 on Linux v6.8.0-1030-azure amd64
16:51:00.751 INFO HaplotypeCaller - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
16:51:00.752 INFO HaplotypeCaller - Start Date/Time: February 11, 2026 at 4:51:00 PM GMT
16:51:00.752 INFO HaplotypeCaller - ------------------------------------------------------------
16:51:00.752 INFO HaplotypeCaller - ------------------------------------------------------------
16:51:00.752 INFO HaplotypeCaller - HTSJDK Version: 4.1.0
16:51:00.753 INFO HaplotypeCaller - Picard Version: 3.1.1
16:51:00.753 INFO HaplotypeCaller - Built for Spark Version: 3.5.0
16:51:00.753 INFO HaplotypeCaller - HTSJDK Defaults.COMPRESSION_LEVEL : 2
16:51:00.753 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
16:51:00.753 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
16:51:00.754 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
16:51:00.754 INFO HaplotypeCaller - Deflater: IntelDeflater
16:51:00.754 INFO HaplotypeCaller - Inflater: IntelInflater
16:51:00.754 INFO HaplotypeCaller - GCS max retries/reopens: 20
16:51:00.754 INFO HaplotypeCaller - Requester pays: disabled
16:51:00.755 INFO HaplotypeCaller - Initializing engine
16:51:00.893 INFO FeatureManager - Using codec BEDCodec to read file file:///data/ref/intervals.bed
16:51:00.905 INFO IntervalArgumentCollection - Processing 6369 bp from intervals
16:51:00.910 INFO HaplotypeCaller - Done initializing engine
16:51:00.912 INFO HaplotypeCallerEngine - Tool is in reference confidence mode and the annotation, the following changes will be made to any specified annotations: 'StrandBiasBySample' will be enabled. 'ChromosomeCounts', 'FisherStrand', 'StrandOddsRatio' and 'QualByDepth' annotations have been disabled
16:51:00.917 INFO NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_utils.so
16:51:00.919 INFO NativeLibraryLoader - Loading libgkl_smithwaterman.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_smithwaterman.so
16:51:00.919 INFO SmithWatermanAligner - Using AVX accelerated SmithWaterman implementation
16:51:00.923 INFO HaplotypeCallerEngine - Standard Emitting and Calling confidence set to -0.0 for reference-model confidence output
16:51:00.923 INFO HaplotypeCallerEngine - All sites annotated with PLs forced to true for reference-model confidence output
16:51:00.933 INFO NativeLibraryLoader - Loading libgkl_pairhmm_omp.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_pairhmm_omp.so
16:51:00.945 INFO IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM
16:51:00.945 INFO IntelPairHmm - Available threads: 4
16:51:00.945 INFO IntelPairHmm - Requested threads: 4
16:51:00.945 INFO PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
16:51:00.984 INFO ProgressMeter - Starting traversal
16:51:00.985 INFO ProgressMeter - Current Locus Elapsed Minutes Regions Processed Regions/Minute
16:51:01.452 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position 20_10037292_10066351:3480 and possibly subsequent; at least 10 samples must have called genotypes
16:51:02.358 INFO HaplotypeCaller - 7 read(s) filtered by: MappingQualityReadFilter
0 read(s) filtered by: MappingQualityAvailableReadFilter
0 read(s) filtered by: MappedReadFilter
0 read(s) filtered by: NotSecondaryAlignmentReadFilter
0 read(s) filtered by: NotDuplicateReadFilter
0 read(s) filtered by: PassesVendorQualityCheckReadFilter
0 read(s) filtered by: NonZeroReferenceLengthAlignmentReadFilter
0 read(s) filtered by: GoodCigarReadFilter
0 read(s) filtered by: WellformedReadFilter
7 total reads filtered out of 1867 reads processed
16:51:02.359 INFO ProgressMeter - 20_10037292_10066351:13499 0.0 35 1529.5
16:51:02.359 INFO ProgressMeter - Traversal complete. Processed 35 total regions in 0.0 minutes.
16:51:02.361 INFO VectorLoglessPairHMM - Time spent in setup for JNI call : 0.0022800000000000003
16:51:02.361 INFO PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.061637120000000004
16:51:02.361 INFO SmithWatermanAligner - Total compute time in native Smith-Waterman : 0.02 sec
16:51:02.362 INFO HaplotypeCaller - Shutting down engine
[February 11, 2026 at 4:51:02 PM GMT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.03 minutes.
Runtime.totalMemory()=257949696
यह कंटेनर में वर्तमान कार्यशील डायरेक्टरी में GVCF आउटपुट फ़ाइल reads_mother.g.vcf बनाता है, साथ ही इसकी इंडेक्स फ़ाइल, reads_mother.g.vcf.idx।
अगर तुम फ़ाइल की सामग्री की पहली 200 लाइनें देखने के लिए head -200 reads_mother.g.vcf चलाते हो, तो तुम देखोगे कि यह सेक्शन 1 में जेनरेट किए गए समकक्ष VCF की तुलना में बहुत लंबा है, और अधिकांश लाइनें VCF में देखी गई चीज़ों से काफी अलग दिखती हैं।
फ़ाइल की सामग्री (संक्षिप्त)
| reads_mother.g.vcf | |
|---|---|
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 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 | |
हमने एक बार फिर अंतिम हेडर लाइन को हाइलाइट किया है, साथ ही फ़ाइल में पहले तीन 'उचित' वेरिएंट कॉल को भी।
तुम देखोगे कि वेरिएंट कॉल लाइनें कई गैर-वेरिएंट लाइनों के बीच बिखरी हुई हैं, जो गैर-वेरिएंट क्षेत्रों का प्रतिनिधित्व करती हैं जहां वेरिएंट कॉलर को भिन्नता का कोई प्रमाण नहीं मिला। जैसा कि ऊपर संक्षेप में उल्लेख किया गया है, यह वेरिएंट कॉलिंग के GVCF मोड में विशेष है: वेरिएंट कॉलर कुछ सांख्यिकी को कैप्चर करता है जो भिन्नता की अनुपस्थिति में इसके विश्वास के स्तर का वर्णन करती हैं। यह दो बहुत अलग मामलों के बीच अंतर करना संभव बनाता है: (1) अच्छी गुणवत्ता वाला डेटा है जो दिखाता है कि नमूना होमोजाइगस-रेफरेंस है, और (2) किसी भी तरह से निर्धारण करने के लिए पर्याप्त अच्छा डेटा उपलब्ध नहीं है।
इस तरह के GVCF में, आमतौर पर ऐसी कई गैर-वेरिएंट लाइनें होती हैं, जिनके बीच कम संख्या में वेरिएंट रिकॉर्ड बिखरे होते हैं।
2.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
कमांड आउटपुट
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar HaplotypeCaller -R /data/ref/ref.fasta -I /data/bam/reads_father.bam -O reads_father.g.vcf -L /data/ref/intervals.bed -ERC GVCF
17:28:30.677 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
17:28:30.801 INFO HaplotypeCaller - ------------------------------------------------------------
17:28:30.803 INFO HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.5.0.0
17:28:30.804 INFO HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
17:28:30.804 INFO HaplotypeCaller - Executing as root@be1a0302f6c7 on Linux v6.8.0-1030-azure amd64
17:28:30.804 INFO HaplotypeCaller - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
17:28:30.804 INFO HaplotypeCaller - Start Date/Time: February 11, 2026 at 5:28:30 PM GMT
17:28:30.804 INFO HaplotypeCaller - ------------------------------------------------------------
17:28:30.804 INFO HaplotypeCaller - ------------------------------------------------------------
17:28:30.805 INFO HaplotypeCaller - HTSJDK Version: 4.1.0
17:28:30.805 INFO HaplotypeCaller - Picard Version: 3.1.1
17:28:30.805 INFO HaplotypeCaller - Built for Spark Version: 3.5.0
17:28:30.806 INFO HaplotypeCaller - HTSJDK Defaults.COMPRESSION_LEVEL : 2
17:28:30.806 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:28:30.806 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:28:30.806 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:28:30.806 INFO HaplotypeCaller - Deflater: IntelDeflater
17:28:30.807 INFO HaplotypeCaller - Inflater: IntelInflater
17:28:30.807 INFO HaplotypeCaller - GCS max retries/reopens: 20
17:28:30.807 INFO HaplotypeCaller - Requester pays: disabled
17:28:30.807 INFO HaplotypeCaller - Initializing engine
17:28:30.933 INFO FeatureManager - Using codec BEDCodec to read file file:///data/ref/intervals.bed
17:28:30.946 INFO IntervalArgumentCollection - Processing 6369 bp from intervals
17:28:30.951 INFO HaplotypeCaller - Done initializing engine
17:28:30.953 INFO HaplotypeCallerEngine - Tool is in reference confidence mode and the annotation, the following changes will be made to any specified annotations: 'StrandBiasBySample' will be enabled. 'ChromosomeCounts', 'FisherStrand', 'StrandOddsRatio' and 'QualByDepth' annotations have been disabled
17:28:30.957 INFO NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_utils.so
17:28:30.959 INFO NativeLibraryLoader - Loading libgkl_smithwaterman.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_smithwaterman.so
17:28:30.960 INFO SmithWatermanAligner - Using AVX accelerated SmithWaterman implementation
17:28:30.963 INFO HaplotypeCallerEngine - Standard Emitting and Calling confidence set to -0.0 for reference-model confidence output
17:28:30.963 INFO HaplotypeCallerEngine - All sites annotated with PLs forced to true for reference-model confidence output
17:28:30.972 INFO NativeLibraryLoader - Loading libgkl_pairhmm_omp.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_pairhmm_omp.so
17:28:30.987 INFO IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM
17:28:30.987 INFO IntelPairHmm - Available threads: 4
17:28:30.987 INFO IntelPairHmm - Requested threads: 4
17:28:30.987 INFO PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
17:28:31.034 INFO ProgressMeter - Starting traversal
17:28:31.034 INFO ProgressMeter - Current Locus Elapsed Minutes Regions Processed Regions/Minute
17:28:31.570 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position 20_10037292_10066351:3480 and possibly subsequent; at least 10 samples must have called genotypes
17:28:32.865 INFO HaplotypeCaller - 9 read(s) filtered by: MappingQualityReadFilter
0 read(s) filtered by: MappingQualityAvailableReadFilter
0 read(s) filtered by: MappedReadFilter
0 read(s) filtered by: NotSecondaryAlignmentReadFilter
0 read(s) filtered by: NotDuplicateReadFilter
0 read(s) filtered by: PassesVendorQualityCheckReadFilter
0 read(s) filtered by: NonZeroReferenceLengthAlignmentReadFilter
0 read(s) filtered by: GoodCigarReadFilter
0 read(s) filtered by: WellformedReadFilter
9 total reads filtered out of 2064 reads processed
17:28:32.866 INFO ProgressMeter - 20_10037292_10066351:13338 0.0 38 1245.2
17:28:32.866 INFO ProgressMeter - Traversal complete. Processed 38 total regions in 0.0 minutes.
17:28:32.868 INFO VectorLoglessPairHMM - Time spent in setup for JNI call : 0.0035923200000000004
17:28:32.868 INFO PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.10765202500000001
17:28:32.868 INFO SmithWatermanAligner - Total compute time in native Smith-Waterman : 0.03 sec
17:28:32.869 INFO HaplotypeCaller - Shutting down engine
[February 11, 2026 at 5:28:32 PM GMT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.04 minutes.
Runtime.totalMemory()=299892736
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
कमांड आउटपुट
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar HaplotypeCaller -R /data/ref/ref.fasta -I /data/bam/reads_son.bam -O reads_son.g.vcf -L /data/ref/intervals.bed -ERC GVCF
17:30:10.017 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
17:30:10.156 INFO HaplotypeCaller - ------------------------------------------------------------
17:30:10.159 INFO HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.5.0.0
17:30:10.159 INFO HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
17:30:10.159 INFO HaplotypeCaller - Executing as root@be1a0302f6c7 on Linux v6.8.0-1030-azure amd64
17:30:10.159 INFO HaplotypeCaller - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
17:30:10.159 INFO HaplotypeCaller - Start Date/Time: February 11, 2026 at 5:30:09 PM GMT
17:30:10.159 INFO HaplotypeCaller - ------------------------------------------------------------
17:30:10.160 INFO HaplotypeCaller - ------------------------------------------------------------
17:30:10.160 INFO HaplotypeCaller - HTSJDK Version: 4.1.0
17:30:10.160 INFO HaplotypeCaller - Picard Version: 3.1.1
17:30:10.161 INFO HaplotypeCaller - Built for Spark Version: 3.5.0
17:30:10.161 INFO HaplotypeCaller - HTSJDK Defaults.COMPRESSION_LEVEL : 2
17:30:10.161 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:30:10.161 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:30:10.161 INFO HaplotypeCaller -HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:30:10.161 INFO HaplotypeCaller - Deflater: IntelDeflater
17:30:10.162 INFO HaplotypeCaller - Inflater: IntelInflater
17:30:10.162 INFO HaplotypeCaller - GCS max retries/reopens: 20
17:30:10.162 INFO HaplotypeCaller - Requester pays: disabled
17:30:10.162 INFO HaplotypeCaller - Initializing engine
17:30:10.277 INFO FeatureManager - Using codec BEDCodec to read file file:///data/ref/intervals.bed
17:30:10.290 INFO IntervalArgumentCollection - Processing 6369 bp from intervals
17:30:10.296 INFO HaplotypeCaller - Done initializing engine
17:30:10.298 INFO HaplotypeCallerEngine - Tool is in reference confidence mode and the annotation, the following changes will be made to any specified annotations: 'StrandBiasBySample' will be enabled. 'ChromosomeCounts', 'FisherStrand', 'StrandOddsRatio' and 'QualByDepth' annotations have been disabled
17:30:10.302 INFO NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_utils.so
17:30:10.303 INFO NativeLibraryLoader - Loading libgkl_smithwaterman.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_smithwaterman.so
17:30:10.304 INFO SmithWatermanAligner - Using AVX accelerated SmithWaterman implementation
17:30:10.307 INFO HaplotypeCallerEngine - Standard Emitting and Calling confidence set to -0.0 for reference-model confidence output
17:30:10.307 INFO HaplotypeCallerEngine - All sites annotated with PLs forced to true for reference-model confidence output
17:30:10.315 INFO NativeLibraryLoader - Loading libgkl_pairhmm_omp.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_pairhmm_omp.so
17:30:10.328 INFO IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM
17:30:10.329 INFO IntelPairHmm - Available threads: 4
17:30:10.329 INFO IntelPairHmm - Requested threads: 4
17:30:10.329 INFO PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
17:30:10.368 INFO ProgressMeter - Starting traversal
17:30:10.369 INFO ProgressMeter - Current Locus Elapsed Minutes Regions Processed Regions/Minute
17:30:10.875 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position 20_10037292_10066351:3480 and possibly subsequent; at least 10 samples must have called genotypes
17:30:11.980 INFO HaplotypeCaller - 14 read(s) filtered by: MappingQualityReadFilter
0 read(s) filtered by: MappingQualityAvailableReadFilter
0 read(s) filtered by: MappedReadFilter
0 read(s) filtered by: NotSecondaryAlignmentReadFilter
0 read(s) filtered by: NotDuplicateReadFilter
0 read(s) filtered by: PassesVendorQualityCheckReadFilter
0 read(s) filtered by: NonZeroReferenceLengthAlignmentReadFilter
0 read(s) filtered by: GoodCigarReadFilter
0 read(s) filtered by: WellformedReadFilter
14 total reads filtered out of 1981 reads processed
17:30:11.981 INFO ProgressMeter - 20_10037292_10066351:13223 0.0 35 1302.7
17:30:11.981 INFO ProgressMeter - Traversal complete. Processed 35 total regions in 0.0 minutes.
17:30:11.983 INFO VectorLoglessPairHMM - Time spent in setup for JNI call : 0.0034843710000000004
17:30:11.983 INFO PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.048108363
17:30:11.983 INFO SmithWatermanAligner - Total compute time in native Smith-Waterman : 0.02 sec
17:30:11.984 INFO HaplotypeCaller - Shutting down engine
[February 11, 2026 at 5:30:11 PM GMT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.03 minutes.
Runtime.totalMemory()=226492416
एक बार यह पूरा हो जाने के बाद, तुम्हारे पास अपनी वर्तमान डायरेक्टरी में .g.vcf में समाप्त होने वाली तीन फ़ाइलें (प्रति नमूना एक) और उनकी संबंधित इंडेक्स फ़ाइलें .g.vcf.idx में समाप्त होनी चाहिए।
डायरेक्टरी सामग्री
इस बिंदु पर, हमने अपने प्रत्येक इनपुट नमूने के लिए GVCF मोड में वेरिएंट कॉल किए हैं। अब संयुक्त कॉलिंग पर आगे बढ़ने का समय है।
लेकिन कंटेनर से बाहर मत निकलो! हम अगले चरण में उसी का उपयोग करने जा रहे हैं।
2.3. संयुक्त जीनोटाइपिंग चलाएं¶
अब जब हमारे पास सभी GVCFs हैं, तो हम नमूनों के कोहॉर्ट के लिए वेरिएंट कॉल जेनरेट करने के लिए संयुक्त जीनोटाइपिंग दृष्टिकोण को आज़मा सकते हैं। यह एक दो-चरणीय विधि है जिसमें सभी GVCFs से डेटा को डेटा स्टोर में जोड़ना शामिल है, फिर संयुक्त-कॉल किए गए वेरिएंट के अंतिम VCF को जेनरेट करने के लिए वास्तविक संयुक्त जीनोटाइपिंग विश्लेषण चलाना।
2.3.1. सभी प्रति-नमूना GVCFs को जोड़ें¶
यह पहला चरण एक अन्य GATK टूल का उपयोग करता है, जिसे GenomicsDBImport कहा जाता है, सभी GVCFs से डेटा को GenomicsDB डेटा स्टोर में जोड़ने के लिए। 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
कमांड आउटपुट
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar 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
17:37:07.569 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
17:37:07.699 INFO GenomicsDBImport - ------------------------------------------------------------
17:37:07.702 INFO GenomicsDBImport - The Genome Analysis Toolkit (GATK) v4.5.0.0
17:37:07.702 INFO GenomicsDBImport - For support and documentation go to https://software.broadinstitute.org/gatk/
17:37:07.703 INFO GenomicsDBImport - Executing as root@be1a0302f6c7 on Linux v6.8.0-1030-azure amd64
17:37:07.703 INFO GenomicsDBImport - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
17:37:07.704 INFO GenomicsDBImport - Start Date/Time: February 11, 2026 at 5:37:07 PM GMT
17:37:07.704 INFO GenomicsDBImport - ------------------------------------------------------------
17:37:07.704 INFO GenomicsDBImport - ------------------------------------------------------------
17:37:07.706 INFO GenomicsDBImport - HTSJDK Version: 4.1.0
17:37:07.706 INFO GenomicsDBImport - Picard Version: 3.1.1
17:37:07.707 INFO GenomicsDBImport - Built for Spark Version: 3.5.0
17:37:07.709 INFO GenomicsDBImport - HTSJDK Defaults.COMPRESSION_LEVEL : 2
17:37:07.709 INFO GenomicsDBImport - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:37:07.709 INFO GenomicsDBImport - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:37:07.710 INFO GenomicsDBImport - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:37:07.710 INFO GenomicsDBImport - Deflater: IntelDeflater
17:37:07.711 INFO GenomicsDBImport - Inflater: IntelInflater
17:37:07.711 INFO GenomicsDBImport - GCS max retries/reopens: 20
17:37:07.711 INFO GenomicsDBImport - Requester pays: disabled
17:37:07.712 INFO GenomicsDBImport - Initializing engine
17:37:07.883 INFO FeatureManager - Using codec BEDCodec to read file file:///data/ref/intervals.bed
17:37:07.886 INFO IntervalArgumentCollection - Processing 6369 bp from intervals
17:37:07.889 INFO GenomicsDBImport - Done initializing engine
17:37:08.560 INFO GenomicsDBLibLoader - GenomicsDB native library version : 1.5.1-84e800e
17:37:08.561 INFO GenomicsDBImport - Vid Map JSON file will be written to /tmp/family_trio_gdb/vidmap.json
17:37:08.561 INFO GenomicsDBImport - Callset Map JSON file will be written to /tmp/family_trio_gdb/callset.json
17:37:08.561 INFO GenomicsDBImport - Complete VCF Header will be written to /tmp/family_trio_gdb/vcfheader.vcf
17:37:08.561 INFO GenomicsDBImport - Importing to workspace - /tmp/family_trio_gdb
17:37:08.878 INFO GenomicsDBImport - Importing batch 1 with 3 samples
17:37:09.359 INFO GenomicsDBImport - Importing batch 1 with 3 samples
17:37:09.487 INFO GenomicsDBImport - Importing batch 1 with 3 samples
17:37:09.591 INFO GenomicsDBImport - Done importing batch 1/1
17:37:09.592 INFO GenomicsDBImport - Import completed!
17:37:09.592 INFO GenomicsDBImport - Shutting down engine
[February 11, 2026 at 5:37:09 PM GMT] org.broadinstitute.hellbender.tools.genomicsdb.GenomicsDBImport done. Elapsed time: 0.03 minutes.
Runtime.totalMemory()=113246208
Tool returned:
true
इस चरण का आउटपुट प्रभावी रूप से एक डायरेक्टरी है जिसमें कई अलग-अलग फ़ाइलों के रूप में संयुक्त वेरिएंट डेटा रखने वाली आगे नेस्टेड डायरेक्टरियों का एक सेट होता है। तुम इसके चारों ओर घूम सकते हो लेकिन तुम जल्दी ही देखोगे कि यह डेटा स्टोर फॉर्मेट मनुष्यों द्वारा सीधे पढ़े जाने के लिए नहीं है।
सुझाव
GATK में ऐसे टूल शामिल हैं जो आवश्यकतानुसार डेटा स्टोर से वेरिएंट कॉल डेटा का निरीक्षण और निष्कर्षण संभव बनाते हैं।
2.3.2. वास्तविक संयुक्त जीनोटाइपिंग विश्लेषण चलाएं¶
यह दूसरा चरण एक और GATK टूल का उपयोग करता है, जिसे GenotypeGVCFs कहा जाता है, कोहॉर्ट में सभी नमूनों में उपलब्ध डेटा के आलोक में वेरिएंट सांख्यिकी और व्यक्तिगत जीनोटाइप्स की पुनर्गणना करने के लिए।
कमांड आउटपुट
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar GenotypeGVCFs -R /data/ref/ref.fasta -V gendb://family_trio_gdb -O family_trio.vcf
17:38:45.084 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
17:38:45.217 INFO GenotypeGVCFs - ------------------------------------------------------------
17:38:45.220 INFO GenotypeGVCFs - The Genome Analysis Toolkit (GATK) v4.5.0.0
17:38:45.220 INFO GenotypeGVCFs - For support and documentation go to https://software.broadinstitute.org/gatk/
17:38:45.220 INFO GenotypeGVCFs - Executing as root@be1a0302f6c7 on Linux v6.8.0-1030-azure amd64
17:38:45.220 INFO GenotypeGVCFs - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
17:38:45.221 INFO GenotypeGVCFs - Start Date/Time: February 11, 2026 at 5:38:45 PM GMT
17:38:45.221 INFO GenotypeGVCFs - ------------------------------------------------------------
17:38:45.221 INFO GenotypeGVCFs - ------------------------------------------------------------
17:38:45.221 INFO GenotypeGVCFs - HTSJDK Version: 4.1.0
17:38:45.222 INFO GenotypeGVCFs - Picard Version: 3.1.1
17:38:45.222 INFO GenotypeGVCFs - Built for Spark Version: 3.5.0
17:38:45.222 INFO GenotypeGVCFs - HTSJDK Defaults.COMPRESSION_LEVEL : 2
17:38:45.222 INFO GenotypeGVCFs - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:38:45.222 INFO GenotypeGVCFs - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:38:45.222 INFO GenotypeGVCFs - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:38:45.223 INFO GenotypeGVCFs - Deflater: IntelDeflater
17:38:45.223 INFO GenotypeGVCFs - Inflater: IntelInflater
17:38:45.223 INFO GenotypeGVCFs - GCS max retries/reopens: 20
17:38:45.223 INFO GenotypeGVCFs - Requester pays: disabled
17:38:45.223 INFO GenotypeGVCFs - Initializing engine
17:38:45.544 INFO GenomicsDBLibLoader - GenomicsDB native library version : 1.5.1-84e800e
17:38:45.561 INFO NativeGenomicsDB - pid=221 tid=222 No valid combination operation found for INFO field InbreedingCoeff - the field will NOT be part of INFO fields in the generated VCF records
17:38:45.561 INFO NativeGenomicsDB - pid=221 tid=222 No valid combination operation found for INFO field MLEAC - the field will NOT be part of INFO fields in the generated VCF records
17:38:45.561 INFO NativeGenomicsDB - pid=221 tid=222 No valid combination operation found for INFO field MLEAF - the field will NOT be part of INFO fields in the generated VCF records
17:38:45.577 INFO GenotypeGVCFs - Done initializing engine
17:38:45.615 INFO ProgressMeter - Starting traversal
17:38:45.615 INFO ProgressMeter - Current Locus Elapsed Minutes Variants Processed Variants/Minute
17:38:45.903 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position 20_10037292_10066351:3480 and possibly subsequent; at least 10 samples must have called genotypes
GENOMICSDB_TIMER,GenomicsDB iterator next() timer,Wall-clock time(s),0.07757032800000006,Cpu time(s),0.07253379200000037
17:38:46.421 INFO ProgressMeter - 20_10037292_10066351:13953 0.0 3390 252357.3
17:38:46.422 INFO ProgressMeter - Traversal complete. Processed 3390 total variants in 0.0 minutes.
17:38:46.423 INFO GenotypeGVCFs - Shutting down engine
[February 11, 2026 at 5:38:46 PM GMT] org.broadinstitute.hellbender.tools.walkers.GenotypeGVCFs done. Elapsed time: 0.02 minutes.
Runtime.totalMemory()=203423744
यह कंटेनर में वर्तमान कार्यशील डायरेक्टरी में VCF आउटपुट फ़ाइल family_trio.vcf बनाता है, साथ ही इसकी इंडेक्स, family_trio.vcf.idx।
यह एक और उचित रूप से छोटी फ़ाइल है, इसलिए तुम फ़ाइल की सामग्री देखने के लिए cat family_trio.vcf चला सकते हो, और पहली कुछ वेरिएंट लाइनों को खोजने के लिए नीचे स्क्रॉल कर सकते हो।
फ़ाइल की सामग्री (संक्षिप्त)
हमने एक बार फिर अंतिम हेडर लाइन को हाइलाइट किया है, जो वेरिएंट कॉल डेटा की शुरुआत को चिह्नित करती है।
यह पहले जेनरेट किए गए VCF के समान दिखता है, सिवाय इस बार हमारे पास सभी तीन नमूनों के लिए जीनोटाइप-स्तरीय जानकारी है। फ़ाइल में अंतिम तीन कॉलम नमूनों के लिए जीनोटाइप ब्लॉक हैं, जो हाइलाइट की गई हेडर लाइन में दिखाए अनुसार उनके ID फ़ील्ड के वर्णानुक्रम में सूचीबद्ध हैं।
अगर हम बहुत पहले वेरिएंट के लिए हमारे परीक्षण परिवार ट्रायो के लिए कॉल किए गए जीनोटाइप्स को देखते हैं, तो हम देखते हैं कि पिता हेटेरोजाइगस-वेरिएंट (0/1) है, और मां और बेटा दोनों होमोजाइगस-वेरिएंट (1/1) हैं।
अंततः यही वह जानकारी है जिसे हम डेटासेट से निकालना चाहते हैं!
2.3.3. आउटपुट फ़ाइलों को स्थानांतरित करें¶
जैसा कि पहले उल्लेख किया गया है, कंटेनर के अंदर जो कुछ भी रहता है वह भविष्य के काम के लिए दुर्गम होगा। कंटेनर से बाहर निकलने से पहले, हम GVCF फ़ाइलों, अंतिम मल्टी-सैंपल VCF और उनकी सभी इंडेक्स फ़ाइलों को कंटेनर के बाहर फ़ाइलसिस्टम में मैनुअल रूप से स्थानांतरित करने जा रहे हैं। इस तरह, हमारे पास तुलना करने के लिए कुछ होगा जब हम इस सभी काम को स्वचालित करने के लिए अपना वर्कफ़्लो बनाएंगे।
डायरेक्टरी की सामग्री" hl_lines="14-19 22-23
data
├── bam
│ ├── reads_father.bam
│ ├── reads_father.bam.bai
│ ├── reads_mother.bam
│ ├── reads_mother.bam.bai
│ ├── reads_son.bam
│ └── reads_son.bam.bai
├── ref
│ ├── intervals.bed
│ ├── ref.dict
│ ├── ref.fasta
│ └── ref.fasta.fai
├── samplesheet.csv
└── vcf
├── family_trio.vcf
├── family_trio.vcf.idx
├── reads_father.g.vcf
├── reads_father.g.vcf.idx
├── reads_mother.g.vcf
├── reads_mother.g.vcf.idx
├── reads_mother.vcf
├── reads_mother.vcf.idx
├── reads_son.g.vcf
└── reads_son.g.vcf.idx
एक बार यह हो जाने के बाद, सभी फ़ाइलें अब तुम्हारी सामान्य फ़ाइलसिस्टम में सुलभ हैं।
2.3.4. GATK कंटेनर से बाहर निकलें¶
कंटेनर से बाहर निकलने के लिए, exit टाइप करो।
तुम्हारा प्रॉम्प्ट वापस सामान्य होना चाहिए। यह संयुक्त वेरिएंट कॉलिंग कमांड्स के मैनुअल परीक्षण को समाप्त करता है।
सारांश¶
तुम जानते हो कि Samtools इंडेक्सिंग और GATK वेरिएंट कॉलिंग कमांड्स को उनके संबंधित कंटेनरों में कैसे परीक्षण करना है, जिसमें GVCFs जेनरेट करना और कई नमूनों पर संयुक्त जीनोटाइपिंग चलाना शामिल है।
आगे क्या है?¶
एक ब्रेक लो, फिर भाग 2 पर जाओ ताकि सीख सको कि उन्हीं कमांड्स को वर्कफ़्लो में कैसे wrap करें जो काम को निष्पादित करने के लिए कंटेनरों का उपयोग करते हैं।