-
Notifications
You must be signed in to change notification settings - Fork 7
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
58 test variant extractor for svs #60
Merged
Merged
Changes from all commits
Commits
Show all changes
9 commits
Select commit
Hold shift + click to select a range
3429248
add variant_extractor
kubranarci 28e3231
add vcf_to_csv.py
kubranarci ee9e870
remove truvari phab
kubranarci e22485b
benchmark params structure
kubranarci 2413c0b
remove unneccsary modules
kubranarci 76aecda
fix params in schema
kubranarci 2264f12
add new samples and fix r plot
kubranarci 81ff2a1
bgzip added to variant_extractor
kubranarci 3deb09d
vcf_to_csv.py changed
kubranarci File filter
Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,5 +1,9 @@ | ||
| test_vcf,caller,vartype | ||
| "https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg37/test/HG002_delly_SV_hg19.chr21.vcf.gz",delly,sv | ||
| "https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg37/test/HG002_lumpy_SV_hg19.sorted.chr21.vcf.gz",lumpy,sv | ||
| "https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg37/test/HG002_manta_SV_hg19_genotype.chr21.vcf",manta,sv | ||
| "https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg37/test/full.svaba.germline.sv.chr21.vcf",svaba,sv | ||
| test_vcf,caller,vartype,refdist,chunksize,normshift,normdist,normsizediff,maxdist,typeignore,dup_to_ins,pctsize,pctseq,pctovl,evaluationmode | ||
| "https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg37/test/HG002_delly_SV_hg19.chr21.vcf.gz",delly,sv,100000,100000,0.3,0.3,0.3,100000,true,true,0.3,0,0,cts | ||
| "https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg37/test/HG002_lumpy_SV_hg19.sorted.vcf.gz",lumpy,sv,100000,100000,0.3,0.3,0.3,100000,true,true,0.3,0,0,cts | ||
| "https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg37/test/HG002_manta_SV_hg19_genotype.chr21.vcf.gz",manta,sv,100000,100000,0.3,0.3,0.3,100000,true,true,0.3,0,0,cts | ||
| "https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg37/test/HG002.svaba.germline.sv.chr21.vcf.gz",svaba,sv,100000,100000,0.3,0.3,0.3,100000,true,true,0.3,0,0,cts | ||
| "https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg37/test/HG002-NA24385-50x.union_170414.split.chr21.vcf.gz",sbg_graph,sv,100000,100000,0.3,0.3,0.3,100000,true,true,0.3,0,0,cts | ||
| "https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg37/test/HG002.jointVC.filter.chr21.vcf.gz",gatk_joint,sv,100000,100000,0.3,0.3,0.3,100000,true,true,0.3,0,0,cts | ||
| "https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg37/test/HG002_ALLCHROM_hs37d5_novoalign_Ilmn150bp300X_FB.chr21.vcf.gz",freebayes,sv,100000,100000,0.3,0.3,0.3,100000,true,true,0.3,0,0,cts | ||
| "https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg37/test/hg002.Assemblytics_structural_variants.sorted.chr21.vcf.gz",assemblytics,sv,100000,100000,0.3,0.3,0.3,100000,true,true,0.3,0,0,cts |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,7 +1,6 @@ | ||
| test_vcf,caller,vartype,pctsize,pctseq,pctovl,refdist,chunksize,normshift,normdist,normsizediff,maxdist | ||
| https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg38/test/ajtrio.lumpy.svtyper.HG002.md.sorted.recal.chr21.vcf.gz,lumpy,sv,0.3,0,0,100000,100000,0.3,0.3,0.3,100000 | ||
| https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg38/test/manta.HG002.chr21.vcf.gz,manta,sv,0.3,0,0,100000,100000,0.3,0.3,0.3,100000 | ||
| https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg38/test/Ashkenazim_HG002.filtered.sv.chr21.vcf.gz,merged,sv,0.3,0,0,100000,100000,0.3,0.3,0.3,100000 | ||
| https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg38/test/HG002_DRAGEN_SV_hg19.chr21.vcf.gz,dragen,sv,0.3,0,0,100000,100000,0.3,0.3,0.3,100000 | ||
| https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/sarek/HG002.strelka.variants.chr21.vcf.gz,strelka,small,,,,,,,,,, | ||
| https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/sarek/HG002.bcftools.chr21.vcf.gz,bcftools,small,,,,,,,,,, | ||
| test_vcf,caller,vartype,pctsize,pctseq,pctovl,refdist,chunksize,normshift,normdist,normsizediff,maxdist,typeignore,evaluationmode | ||
| https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg38/test/manta.HG002.chr21.vcf.gz,manta,sv,0.3,0,0,100000,100000,0.3,0.3,0.3,100000,true,cts | ||
| https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg38/test/Ashkenazim_HG002.filtered.sv.chr21.vcf.gz,merged,sv,0.3,0,0,100000,100000,0.3,0.3,0.3,100000,true,cts | ||
| https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/SV_testdata/hg38/test/HG002_DRAGEN_SV_hg19.chr21.vcf.gz,dragen,sv,0.3,0,0,100000,100000,0.3,0.3,0.3,100000,true,cts | ||
| https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/sarek/HG002.strelka.variants.chr21.vcf.gz,strelka,small,,,,,,,,,,,, | ||
| https://raw.githubusercontent.com/kubranarci/benchmark_datasets/main/sarek/HG002.bcftools.chr21.vcf.gz,bcftools,small,,,,,,,,,,,, |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,29 @@ | ||
| #!/usr/bin/env python | ||
|
|
||
| import pysam | ||
| from argparse import ArgumentParser | ||
|
|
||
| def compare_vcf(file1, file2): | ||
| vcf1 = pysam.VariantFile(file1) | ||
| vcf2 = pysam.VariantFile(file2) | ||
|
|
||
| variants1 = {rec.id for rec in vcf1.fetch()} | ||
| variants2 = {rec.id for rec in vcf2.fetch()} | ||
|
|
||
| unique_to_file1 = variants1 - variants2 | ||
| unique_to_file2 = variants2 - variants1 | ||
|
|
||
| print(f"Unique to {file1}: {unique_to_file1}") | ||
| print(f"Unique to {file2}: {unique_to_file2}") | ||
|
|
||
| if __name__ == '__main__': | ||
| import os | ||
| import sys | ||
|
|
||
| # Parse arguments | ||
| parser = ArgumentParser(description='Check differences') | ||
| parser.add_argument('file1', help='VCF file') | ||
| parser.add_argument('file2', help='Output VCF file') | ||
|
|
||
| args = parser.parse_args() | ||
| compare_vcf(args.file1, args.file2) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,58 @@ | ||
| #!/usr/bin/env python | ||
|
|
||
| # Copyright 2022 - Barcelona Supercomputing Center | ||
| # Author: Rodrigo Martin | ||
| # BSC Dual License | ||
| ''' | ||
| Generates a normalized VCF file from a VCF | ||
| Expected usage: | ||
| $ python normalize.py <vcf_file> <output_vcf_file> | ||
| Use --help for more information. | ||
| ''' | ||
| import sys | ||
| from argparse import ArgumentParser | ||
|
|
||
| import pysam | ||
|
|
||
|
|
||
| def _contig_to_int(contig): | ||
| contig = contig.lower().replace('chr', '') | ||
| if contig.isdigit(): | ||
| return int(contig) | ||
| else: | ||
| return 22 + ord(contig[0]) | ||
|
|
||
|
|
||
| if __name__ == '__main__': | ||
| import os | ||
| import sys | ||
| sys.path.insert(0, os.path.abspath(os.path.dirname(__file__)) + '/../src/') | ||
| from variant_extractor import VariantExtractor | ||
| from variant_extractor.variants import VariantType | ||
|
|
||
| # Parse arguments | ||
| parser = ArgumentParser(description='Generate normalized VCF file from a VCF file') | ||
| parser.add_argument('vcf_file', help='VCF file') | ||
| parser.add_argument('output_vcf_file', help='Output VCF file') | ||
| parser.add_argument('fasta', help='Reference fasta file coupled with fai') | ||
|
|
||
| args = parser.parse_args() | ||
|
|
||
| # Extract header from original file | ||
| with open(args.vcf_file, 'r') as f: | ||
| with pysam.VariantFile(f) as input_file: | ||
| input_file.header.add_meta('cmdline', ' '.join(sys.argv)) | ||
| header_str = str(input_file.header) | ||
|
|
||
| # Open output file, write as stream | ||
| with open(args.output_vcf_file, 'w') as output_vcf: | ||
| # Write header | ||
| output_vcf.write(header_str) | ||
| print(f'Reading {args.vcf_file}...') | ||
| # Open input file, read with variant_extractor | ||
| extractor = VariantExtractor(args.vcf_file,ensure_pairs=False,fasta_ref=args.fasta) | ||
| records = list(extractor) | ||
| # Sort record by chromosome and position | ||
| records.sort(key=lambda x: (_contig_to_int(x.contig), x.pos)) | ||
| for variant_record in records: | ||
| output_vcf.write(str(variant_record)+'\n') |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,117 @@ | ||
| #!/usr/bin/env python | ||
|
|
||
| # Copyright 2024 - GHGA | ||
| # Author: Kuebra Narci | ||
| ''' | ||
| Generates a CSV file from a VCF | ||
| Expected usage: | ||
| $ python vcf_to_csv.py <vcf_file> <output> | ||
| Use --help for more information. | ||
| ''' | ||
| import sys | ||
| from argparse import ArgumentParser | ||
| import re | ||
|
|
||
| import csv | ||
|
|
||
| def parse_info_field(info): | ||
| """Parse the INFO field of a VCF line into a dictionary.""" | ||
| info_dict = {} | ||
| for entry in info.split(';'): | ||
| key_value = entry.split('=') | ||
| if len(key_value) == 2: | ||
| info_dict[key_value[0]] = key_value[1] | ||
| else: | ||
| info_dict[key_value[0]] = True | ||
| return info_dict | ||
|
|
||
| def extract_gt_from_sample(sample, format_field): | ||
| """Extract GT value from the sample field.""" | ||
| format_fields = format_field.split(':') | ||
| sample_values = sample.split(':') | ||
| if 'GT' in format_fields: | ||
| gt_index = format_fields.index('GT') | ||
| return sample_values[gt_index] | ||
| return './.' # Default GT value if not present | ||
|
|
||
| def vcf_to_csv(vcf_file, csv_file): | ||
| """Convert a VCF file to a CSV file.""" | ||
| with open(vcf_file, 'r') as vcf: | ||
| headers = [] | ||
| sample_headers = [] | ||
| records = [] | ||
| include_supp_vec = False | ||
| include_supp = False | ||
| include_type_inferred = False | ||
| include_svtype = False | ||
| include_svlen = False | ||
|
|
||
| for line in vcf: | ||
| if line.startswith('##'): | ||
| continue # Skip meta-information lines | ||
| elif line.startswith('#'): | ||
| headers = line[1:].strip().split('\t') | ||
| sample_headers = headers[9:] # The sample headers start from the 10th column | ||
| else: | ||
| row = line.strip().split('\t') | ||
| info_dict = parse_info_field(row[7]) | ||
|
|
||
| # Check for SUPP_VEC, SUPP, type_inferred, SVTYPE, and SVLEN in the INFO field | ||
| if 'SUPP_VEC' in info_dict: | ||
| include_supp_vec = True | ||
| if 'SUPP' in info_dict: | ||
| include_supp = True | ||
| if 'type_inferred' in info_dict: | ||
| include_type_inferred = True | ||
| if 'SVTYPE' in info_dict: | ||
| include_svtype = True | ||
| if 'SVLEN' in info_dict: | ||
| include_svlen = True | ||
|
|
||
| records.append((row, info_dict)) | ||
|
|
||
| # Write the header with optional fields | ||
| headers_to_write = headers[:7] # Only keep CHROM, POS, ID, REF, ALT, QUAL, FILTER | ||
| if include_supp_vec: | ||
| headers_to_write.append("SUPP_VEC") | ||
| if include_supp: | ||
| headers_to_write.append("SUPP") | ||
| if include_type_inferred: | ||
| headers_to_write.append("type_inferred") | ||
| if include_svtype: | ||
| headers_to_write.append("SVTYPE") | ||
| if include_svlen: | ||
| headers_to_write.append("SVLEN") | ||
| headers_to_write.extend([f'{sample}_GT' for sample in sample_headers]) | ||
|
|
||
| with open(csv_file, 'w', newline='') as csvf: | ||
| csv_writer = csv.writer(csvf) | ||
| csv_writer.writerow(headers_to_write) | ||
|
|
||
| for row, info_dict in records: | ||
| row_to_write = row[:7] # Only keep CHROM, POS, ID, REF, ALT, QUAL, FILTER | ||
| if include_supp_vec: | ||
| row_to_write.append(info_dict.get('SUPP_VEC', '')) | ||
| if include_supp: | ||
| row_to_write.append(info_dict.get('SUPP', '')) | ||
| if include_type_inferred: | ||
| row_to_write.append(info_dict.get('type_inferred', '')) | ||
| if include_svtype: | ||
| row_to_write.append(info_dict.get('SVTYPE', '')) | ||
| if include_svlen: | ||
| row_to_write.append(info_dict.get('SVLEN', '')) | ||
| format_field = row[8] | ||
| gt_values = [extract_gt_from_sample(sample, format_field) for sample in row[9:]] | ||
| row_to_write.extend(gt_values) | ||
|
|
||
| csv_writer.writerow(row_to_write) | ||
|
|
||
| if __name__ == '__main__': | ||
|
|
||
| # Parse arguments | ||
| parser = ArgumentParser(description='Generates a CSV file from a VCF') | ||
| parser.add_argument('vcf_file', help='VCF file') | ||
| parser.add_argument('output', help='Output CSV file') | ||
| args = parser.parse_args() | ||
|
|
||
| vcf_to_csv(args.vcf_file, args.output) |
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You can keep these. they won't be validated if the value has not been given
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I wanted to be safe 🗡️ ahaha
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Don't worry it won't break anything :)