diff --git a/README.md b/README.md index 8aedafd..03abbc8 100644 --- a/README.md +++ b/README.md @@ -19,11 +19,9 @@ source activate vcf_parse ``` -usage: vcf_parse.py [-h] [-v] [-O OUTPUT] [-t TRANSCRIPTS] - [-T TRANSCRIPT_STRICTNESS] [-b BED | -B BED_FOLDER] - [-k KNOWN_VARIANTS] [-c CONFIG] [-l] - input - +usage: vcf_parse.py [-h] [-v] [-i INPUT VCF] [-O OUTPUT] + [-t TRANSCRIPTS] [-n NTC] + summary: Takes a VCF file and parses the variants to produce a tab delimited variant report. @@ -34,109 +32,17 @@ positional arguments: optional arguments: -h, --help show this help message and exit -v, --version show program's version number and exit - -O OUTPUT, --output OUTPUT + -i INPUT VCF --inp INPUT VCF + VCF file and filepath + + -o OUTPUT, --outfile OUTPUT Filepath to folder where output reports will be saved. If missing, defaults to current directory. - -t TRANSCRIPTS, --transcripts TRANSCRIPTS - - Filepath to preferred transcripts file. - - Must be a tab seperated file with preferred transcripts in the second - column. If missing, all entries in the preferred transcript column - will be labelled as 'Unknown'. - - -T TRANSCRIPT_STRICTNESS, --transcript_strictness TRANSCRIPT_STRICTNESS - - Strictness of matching while annotating preferred transcripts. - Default setting is low. - - Options: - - high - Transcripts must be an exact match. - e.g. NM_001007553.2 and NM_001007553.1 won't match, - NM_001007553.1 and NM_001007553.1 will. - - low - Transcripts will match regardless of the version number. The - version number is after the . at the end of a transcript - e.g. NM_001007553.2 and NM_001007553.1 will match. - - -b BED, --bed BED - Filepath to a single BED file. - - The BED file will be applied to the variant report and a seperate - report saved with the BED file applied. This report will be saved in - the same output folder as the original variant report, with the BED - file name added to it. - Cannot be used together with -B flag. - - -B BED_FOLDER, --bed_folder BED_FOLDER - - Filepath to folder containing BED files. - - Each BED file will be applied to the variant report and a seperate - report saved with the BED file applied. These reports will be saved in - a new folder within the output folder, named the same as the input BED - folder. - The file names will be the same as the original variant report, with - the BED file name added to them. - Cannot be used together with -b flag. - - -k KNOWN_VARIANTS, --known_variants KNOWN_VARIANTS - - Filepath to known variants file. - - This is a VCF file containing any known variants and an associated - classification. The classification will be added to the variant - report. The VCF must have an annotation named 'Classification' within - the INFO field for each variant. - - Key: - 0 - Artifact - 1 - Benign - 2 - Likely benign - 3 - VUS - 4 - Likely pathogenic - 5 - Pathogenic - - -c CONFIG, --config CONFIG - - Filepath to config file. - - This is a tab seperated text file containing a number of rows, where - each row specifies an annotation to be included in the variant report. - Only annotations included in the config file will be included in the - variant report. - The columns in the variant report will be in the same order as the - order in which the annotations appear in the config file. - - Each row contains: - - Column 1 - Required. Annotation headers, these must match up with how - they appear in the VCF (case sensitive). - - Column 2 - Required. Location where to find the data within the VCF, - used to select the correct parsing function. - options: info, format, vep, filter or pref. - - Column 3 - Optional. Alternative name for column header. - - To make a config file with all available options from a VCF, run: - vcf_parse -l path_to_input_vcf > config.txt - - -l, --config_list - Return a list of all availabile config to the screen, then exit. - See CONFIG section for usage. +-t PREFERRED TRANSCRIPT --pref_trans PREFERRED TRANSCRIPT + Preferred transcript file and filepath ``` -### Special formatting for report fields: -- % Allele frequency has been calculated from the AD -- Genotypes are reformatted from 0/0, 0/1 and 1/1 to HOM_REF, HET and HOM_ALT, respectively, to prevent them from appearing as dates in Excel -- dbSNP, COSMIC and HGMD are parsed from the 'Existing_variation' VEP field -- All population frequencies from ExAC and 1KG are re-formatted to a percentage -- Intron and exon numberings are renamed from x/y to x|y (where x and y are numbers), to prevent them appearing as dates in Excel -- HGVS coding (HGVSc) and protein (HGVSp) sequences have had the transcript name trimmed off - ### Testing To run unit tests run `python -m unittest test` diff --git a/config/somatic_amplicon_config.txt b/config/somatic_amplicon_config.txt deleted file mode 100644 index 0dd2eb0..0000000 --- a/config/somatic_amplicon_config.txt +++ /dev/null @@ -1,31 +0,0 @@ -Frequency format -DP info Depth -GT format Genotype -GQ format Quality -Classification custom -Preferred pref -dbSNP vep -Cosmic vep -HGMD vep -ExAC_AFR_MAF vep ExAC_African -ExAC_AMR_MAF vep ExAC_American -ExAC_NFE_MAF vep ExAC_EuropeanNonFinnish -ExAC_FIN_MAF vep ExAC_Finnish -ExAC_EAS_MAF vep ExAC_EastAsian -ExAC_SAS_MAF vep ExAC_SouthAsian -ExAC_OTH_MAF vep ExAC_Other -AFR_MAF vep 1KG_African -AMR_MAF vep 1KG_American -EUR_MAF vep 1KG_European -EAS_MAF vep 1KG_EastAsian -SAS_MAF vep 1KG_SouthAsian -SYMBOL vep Gene -Feature vep Transcript -HGVSc vep -HGVSp vep -Consequence vep -INTRON vep -EXON vep -SIFT vep -PolyPhen vep -CANONICAL vep diff --git a/env/vcf_parse.yaml b/env/vcf_parse.yaml index f8ae6cd..5e38d43 100644 --- a/env/vcf_parse.yaml +++ b/env/vcf_parse.yaml @@ -1,8 +1,9 @@ -name: vcf_parse +name: vcf_parse_somamp channels: - bioconda dependencies: -- python=2.7 +- python=3.9.16 - bedtools=2.27.1 - pip: - - pyvcf==0.6.8 \ No newline at end of file + - pyvariantfilter==1.6.2 + - pysam==0.20.0 diff --git a/scripts/vcf_report.py b/scripts/vcf_report.py index 8f04c57..c106cb4 100644 --- a/scripts/vcf_report.py +++ b/scripts/vcf_report.py @@ -8,8 +8,8 @@ Author: Erik Waskiewicz Created: 31 Aug 2018 -Version: 0.1.0 -Updated: 31 Oct 2018 +Edited: Niamh Teague - 25 May 2023 +Version: 0.1.1 """ @@ -37,6 +37,7 @@ def load_data(self, inp, out): 'loading VCF file from {}'.format(os.path.abspath(inp))) with open(inp, 'r') as vcf_input: vcf_reader = vcf.Reader(vcf_input) + print(vcf_reader) vcf_records = [] for var in vcf_reader: vcf_records.append(var) @@ -46,7 +47,7 @@ def load_data(self, inp, out): # load sample name from vcf self.sample = vcf_reader.samples[0] - # load info and format fields + # load info, variant and format fields self.info_fields = vcf_reader.infos self.format_fields = vcf_reader.formats @@ -68,7 +69,6 @@ def load_data(self, inp, out): # make empty vep variable self.config = None - def load_config(self, config_file): """ Load in config file that defines what annotations to include @@ -93,51 +93,76 @@ def list_config(self): Returns a list to screen containing all possible column headers and the source that the data comes from. """ - # preferred transcript and filter fields - print('Preferred\tpref') - print('Filter\tfilter') - print('Classification\tcustom') - # info fields - for record in self.info_fields: - print(record + '\tinfo') + # Chrom, pos, alt, ref fields + print('SYMBOL\tvep\tgene') + print('CHROM\tvar\tchr') + print('POS\tvar\tpos') + print('REF\tvar\tref') + print('ALT\tvar\talt') + + # preferred transcript and filter fields # format fields - print('Frequency\tformat') - for record in self.format_fields: - print(record + '\tformat') + print('Frequency\tformat\tvaf') + + # info fields - only want depth + for record in self.info_fields: + if record == 'DP': + print('DP' + '\tinfo' + '\tdepth') # vep fields - print('dbSNP\tvep') - print('Cosmic\tvep') - print('HGMD\tvep') - for record in self.vep_fields: - print(record + '\tvep') + print('Feature\tvep\ttranscript') + print('HGVSp\tvep\thgvs_p') + print('HGVSc\tvep\thgvs_c') + print('Consequence\tvep\tconsequence') + print('EXON\tvep\texon') + + def parse_chrom_field(self, variant): - def make_variant_name(self, variant): - var_name = '{}:{}{}>{}'.format( - str(variant.CHROM), - str(variant.POS), - str(variant.REF), - str(variant.ALT).strip('[]').replace(' ', '') - ) - return(var_name) + for variant in self.data: + out = [variant.CHROM] + else: + out = [''] + #print(out) + return(out) - def parse_filter_field(self, variant): - try: - out = [str(variant.FILTER)] - except: + def parse_pos_field(self, variant): + for variant in self.data: + out = [variant.POS] + else: + out = [''] + + return(out) + + + def parse_ref_field(self, variant): + for variant in self.data: + out = [variant.REF] + else: out = [''] + return(out) - def parse_info_field(self, variant, field): + def parse_alt_field(self, variant): + for variant in self.data: + out = variant.ALT + else: + out = [''] + + return(out) + + + def parse_info_field(self, var, field): try: - out = [str(variant.INFO[field])] + out = [str(var.INFO[field])] + except: out = [''] + return(out) @@ -148,7 +173,7 @@ def parse_format_field(self, variant, field): out = [ sample[field] ] except: out = [''] - + # custom setting for allele freq if field == 'Frequency': out = [ sample['AD'] ] @@ -156,28 +181,17 @@ def parse_format_field(self, variant, field): alt = float(out[0][1]) freq = float((alt / (ref + alt)) * 100) out = ['{}%'.format(round(freq, 2))] - - # custom setting for genotype - if field == 'GT': - gt = out[0] - if gt == '0/1': - out = ['HET'] - if gt == '1/1': - out = ['HOM_VAR'] - if gt == '0/0': - out = ['HOM_REF'] return(out) - def parse_vep_field(self, field, vep): if vep: try: - pos = self.vep_fields.index(field) - out = str(vep[pos]) + pos = self.vep_fields.index(field) #chrom? + out = str(vep[pos]) #gene except: - out = '' - + out = '' + # custom exon/intron tweak if field in ('EXON', 'INTRON'): out = out.replace('/', '|') @@ -189,73 +203,24 @@ def parse_vep_field(self, field, vep): out = split[1] except: pass - - # custom existing variantion field tweak - if field in ('dbSNP', 'Cosmic', 'HGMD'): - # split existing variation field by & sign - existing_variation_pos = self.vep_fields.index('Existing_variation') - existing_variation = str(vep[existing_variation_pos]).split('&') - - # define identifier for each different annotation type - if field == 'dbSNP': - id = 'rs' - if field == 'Cosmic': - id = 'COSM' - if field == 'HGMD': - id = 'CM' - - # make output string containing only records of the desired type - out_list = '' - for item in existing_variation: - if item.startswith(id): - out_list += '{},'.format(str(item)) - out = out_list.rstrip(',') - - # custom exac/1kg tweak - if field in ('ExAC_AFR_MAF', 'ExAC_AMR_MAF', 'ExAC_EAS_MAF', 'ExAC_FIN_MAF', - 'ExAC_NFE_MAF', 'ExAC_SAS_MAF', 'ExAC_OTH_MAF', 'AFR_MAF', 'AMR_MAF', - 'EAS_MAF', 'EUR_MAF', 'SAS_MAF'): - - out_string = '' - - try: - for record in out.split('&'): - split = record.split(':') - percent = float(split[1]) * 100 - out_string += '{}:{}%,'.format(split[0], str(percent)) - out = out_string.rstrip(',') - except: - pass - + else: out = 'No VEP output' - + return([out]) def make_header(self): - # Sample and variant are always the first two columns - header = 'SampleID\tVariant' - + # gene to be first column # config file provided - if self.config: - for annotation in self.config: - if annotation[2] != '': - header += '\t' + annotation[2] - else: - header += '\t' + annotation[0] - - # config file not provided - all headers - else: - header += '\tPreferred\tClassification\tFilter' - for annotation in self.info_fields: - if annotation != 'CSQ': - header += '\t' + annotation - for annotation in self.format_fields: - header += '\t' + annotation - for annotation in self.vep_fields: - header += '\t' + annotation - + header = '' + #if self.config: + for annotation in self.config: + if annotation[2] != '': + header += '\t' + annotation[2] + else: + header += '\t' + annotation[0] + # add newline and return header += '\n' return(header) @@ -267,6 +232,22 @@ def make_record_config(self, setting, variant, vep=None): """ out = [''] + # chrom + if setting[2] == 'chr': + out = self.parse_chrom_field(setting[0]) + + # pos + if setting[2] == 'pos': + out = self.parse_pos_field(setting[0]) + + # ref + if setting[2] == 'ref': + out = self.parse_ref_field(setting[0]) + + # alt + if setting[2] == 'alt': + out = self.parse_alt_field(setting[0]) + # preferred if setting[1] == 'pref': out = ['Unknown'] @@ -349,23 +330,24 @@ def make_report(self, filter_setting): # loop through variants for var in self.data: - + # PASS filter - pass will be empty - [], anything else will be filtered out if filter_setting and var.FILTER : pass else: # make variant name - variant = self.make_variant_name(var) - + #variant = self.make_variant_name(var) + # if VEP annotation exists, loop through each transcript try: vep = var.INFO['CSQ'] + for record in range(len(vep)): out = [] vep_split = vep[record].split('|') - # filter out any transcripts that dont begin with NM + # filter out any transcripts that dont begin with NM transcript_col = self.vep_fields.index('Feature') if vep_split[transcript_col].startswith('NM'): @@ -378,9 +360,9 @@ def make_report(self, filter_setting): # filter and preferred must be first, in that order else: out = self.make_record_no_config(var, vep=vep_split) - + # save to file then repeat for all transcripts - report_writer.writerow([self.sample] + [variant] + out) + report_writer.writerow(out) # if variant has no vep annotations except: @@ -390,18 +372,18 @@ def make_report(self, filter_setting): if self.config: for annotation in self.config: out += self.make_record_config(annotation, var) - + # if no config file - include all annotations # filter and preferred must be first, in that order else: out = self.make_record_no_config(var) # save to file then repeat for next variant - report_writer.writerow([self.sample] + [variant] + out) + report_writer.writerow(out) # once loop has finished, close the output file outfile.close() - + # remove duplicates in output file uniq = os.popen('cat {} | uniq'.format(self.report_path)).read() diff --git a/vcf_parse.py b/vcf_parse.py index 82f0731..3e349ee 100755 --- a/vcf_parse.py +++ b/vcf_parse.py @@ -1,313 +1,329 @@ -#!/anaconda3/envs/python2/bin/python +#!/usr/bin/env python """ -vcf_parse.py - -Takes a VCF file and parses the variants to produce a tab delimited -variant report. - -Usage: vcf_parse.py [-h] [-v] - [-O OUTPUT] - [-t TRANSCRIPTS] [-T TRANSCRIPT_STRICTNESS] - [-b BED | -B BED_FOLDER] - [-k KNOWN_VARIANTS] - [-c CONFIG] [-l] - input - vcf_parse.py -h for full description of options. - -Author: Erik Waskiewicz -Created: 31 Aug 2018 -Version: 0.1.0 -Updated: 31 Oct 2018 +Convert annotated vcf from 1_SomaticAmplicon.sh into readable variant report file """ -__version__ = '0.1.0' -__updated__ = '31 Oct 2018' - - -import argparse +import re +import os +import csv import logging +import argparse import textwrap +from pysam import VariantFile +from pyvariantfilter.utils import parse_csq_field -from scripts.vcf_report import vcf_report -from scripts.preferred_transcripts import preferred_transcripts -from scripts.bed_object import bed_object -from scripts.known_variants import known_variants +# Main functions to create variant report file +# """ +# Logger deals with all status messages +# """ +logger = logging.getLogger('vcf_parse.vcf') -## -- PARSE INPUT ARGUMENTS ------------------------------------------- +def parse_transcripts(pref_trans): + """ + Get list of preferred transcripts + """ -def get_args(): - """ - Use argparse package to take arguments from the command line. - See descriptions for full detail of each argument. - """ - - # Make argparse object, add description - parser = argparse.ArgumentParser( - formatter_class=argparse.RawTextHelpFormatter, - description=textwrap.dedent( - ''' - summary: - Takes a VCF file and parses the variants to produce a tab delimited - variant report. - ''' - )) - - - # Version info - parser.add_argument( - '-v', '--version', action='version', - version= - '%(prog)s\nversion:\t{}\nlast updated:\t{}'.format( - __version__, __updated__ - )) - - - # Arguments (see help string for full descriptions): - # REQUIRED: VCF file input - parser.add_argument( - 'input', action='store', - help='Filepath to input VCF file. REQUIRED.' - ) - - - # OPTIONAL: Output folder, defaults to current directory if empty - parser.add_argument( - '-O', '--output', action='store', - help=textwrap.dedent( - ''' - Filepath to folder where output reports will be saved. - If missing, defaults to current directory. - \n''' - )) - - - # OPTIONAL: List of preferred transcripts - parser.add_argument( - '-t', '--transcripts', action='store', - help=textwrap.dedent( - ''' - Filepath to preferred transcripts file. - - Must be a tab seperated file with preferred transcripts in the second - column. If missing, all entries in the preferred transcript column - will be labelled as 'Unknown'. - \n''' - )) - - - # OPTIONAL: Preferred transcripts strictness - parser.add_argument( - '-T', '--transcript_strictness', action='store', default='low', - help=textwrap.dedent( - ''' - Strictness of matching while annotating preferred transcripts. - Default setting is low. - - Options: - - high - Transcripts must be an exact match. - e.g. NM_001007553.2 and NM_001007553.1 won't match, - NM_001007553.1 and NM_001007553.1 will. - - low - Transcripts will match regardless of the version number. The - version number is after the . at the end of a transcript - e.g. NM_001007553.2 and NM_001007553.1 will match. - \n''' - )) - - - # OPTIONAL: either a single BED file or a folder containing BED - # files, only one of these can be used - bed_files = parser.add_mutually_exclusive_group() - - # Single BED file - bed_files.add_argument( - '-b', '--bed', action='store', - help=textwrap.dedent( - ''' - Filepath to a single BED file. - - The BED file will be applied to the variant report and a seperate - report saved with the BED file applied. This report will be saved in - the same output folder as the original variant report, with the BED - file name added to it. - Cannot be used together with -B flag. - \n''' - )) - - # Multiple BED files - bed_files.add_argument( - '-B', '--bed_folder', action='store', - help=textwrap.dedent( - ''' - Filepath to folder containing BED files. - - Each BED file will be applied to the variant report and a seperate - report saved with the BED file applied. These reports will be saved in - a new folder within the output folder, named the same as the input BED - folder. - The file names will be the same as the original variant report, with - the BED file name added to them. - Cannot be used together with -b flag. - \n''' - )) - - - # OPTIONAL: File containing known variants - parser.add_argument( - '-k', '--known_variants', action='store', - help=textwrap.dedent( - ''' - Filepath to known variants file. - - This is a VCF file containing any known variants and an associated - classification. The classification will be added to the variant - report. The VCF must have an annotation named 'Classification' within - the INFO field for each variant. - - Key: - 0 - Artifact - 1 - Benign - 2 - Likely benign - 3 - VUS - 4 - Likely pathogenic - 5 - Pathogenic - \n''' - )) - - - # OPTIONAL: File containing the headers for the report - parser.add_argument( - '-c', '--config', action='store', - help=textwrap.dedent( - ''' - Filepath to config file. - - This is a tab seperated text file containing a number of rows, where - each row specifies an annotation to be included in the variant report. - Only annotations included in the config file will be included in the - variant report. - The columns in the variant report will be in the same order as the - order in which the annotations appear in the config file. - - Each row contains: - - Column 1 - Required. Annotation headers, these must match up with how - they appear in the VCF (case sensitive). - - Column 2 - Required. Location where to find the data within the VCF, - used to select the correct parsing function. - options: info, format, vep, filter or pref. - - Column 3 - Optional. Alternative name for column header. - - To make a config file with all available options from a VCF, run: - vcf_parse -l path_to_input_vcf > config.txt - \n''' - )) - - - # OPTIONAL: Lists all headers in a vcf then exits - parser.add_argument( - '-l', '--config_list', action='store_true', - help=textwrap.dedent( - ''' - Return a list of all availabile config to the screen, then exit. - See CONFIG section for usage. - \n''' - )) - - - # OPTIONAL: Filter out any variants where FILTER column is not PASS - parser.add_argument( - '-F', '--filter_non_pass', action='store_true', - help=textwrap.dedent( - ''' - Filters out any variants where the FILTER annotation is not - PASS. If missing then there will be no fitering based on the - FILTER annotation. - \n''' - )) - - return parser.parse_args() - - -# -- MAIN FUNCTION ---------------------------------------------------- + #Empty variable + pref_transcript = [] + + #Make list of transcripts + trans = open(pref_trans, 'r') + + for line in trans: + + #Skip header + if line.startswith('Transcript'): + next + else: + transcript = line.strip().split('\t')[0] + pref_transcript.append(transcript) + trans.close() + + return(pref_transcript) + + logger.info('Loaded preferred transcripts') + +def parse_ntc(ntc_file): + """ + Create dictionary of all NTC variants + """ + + # Empty ntc list + ntc_list = [] + + # Loop over the ntc file + ntc = VariantFile(ntc_file) + + for line in ntc.fetch(): + + chrom = line.chrom + pos = line.pos + ref = line.ref + alt = line.alts[0] + var = chrom+":"+str(pos)+ref+">"+alt + vaf = line.info['AF'] + # Convert vaf tuple to string and remove parentheses + vaf = str(vaf).replace('(','').replace(')','').replace(',','') + dp = line.samples[0]['DP'] + + #For alt read count, get last value from AD + alt = line.samples[0]['AD'][1] + + ntc_var = [var, vaf, dp, alt] + ntc_list.append(ntc_var) + + return(ntc_list) + + logger.info('created list of NTC variants') + +def parse_sample(inp, outfile, ntc_list, pref_transcript): + """ + Loads in data from vcf + Returns list of variants in correct format + Only gives ntc_vaf and ntc_depth if in_ntc = True + """ + + # Read filtered meta annotated vcf produced in 1_SomaticAmplicon.sh + logger.info( + 'loading VCF file from {}'.format(os.path.abspath(inp))) + + # Empty list to hold variants + variant_list = [] + + # Read vcf file + inp = VariantFile(inp) + # CSQ fields from header to be parsed with pyvariantfilter + csq_fields = str(inp.header.info['CSQ'].record) + csq_fields = csq_fields.strip() + index = csq_fields.index('Format') + 8 + csq_fields = csq_fields[index:len(csq_fields)-2].split('|') + + for record in inp.fetch(): + + if 'CSQ' in record.info: + + # Load chrom, pos, ref & alt fields + chrom = record.chrom + pos = record.pos + ref = record.ref + alt = record.alts[0] + + # Load vaf and depth + vaf = record.info['AF'] + # Convert vaf tuple to string and remove parentheses + vaf = str(vaf).replace('(','').replace(')','').replace(',','') + depth = record.info['DP'] + # Alt-reads - take second value from AD + alt_reads = record.samples[0]['AD'][1] + + # Sample for output filepath + + # Get annotations and go through these + csq = record.info['CSQ'] + annotations = parse_csq_field(csq, csq_fields) + + # Empty list of annotations + transcripts = [] + # If any of the 'Feature' values match the preferred transcript + for entry in annotations: + + # Loop over annotations and if preferred transcript add to transcript list + if entry['Feature'] in pref_transcript: + gene = entry['SYMBOL'] + + # If transcripts for variant still empty, add MANE select + if len(transcripts) == 0: + + for entry in annotations: + + if entry['MANE']: + + gene = entry['SYMBOL'] + hgvs_p = entry['HGVSp'] + hgvs_c = entry['HGVSc'] + consequence = entry['Consequence'] + exon = entry['EXON'] + + transcript_info = (gene + , hgvs_p + , hgvs_c + , consequence + , exon) + + transcripts.append(transcript.info) + + # If transcripts still empty, use PICK annotation + if len(transcripts) == 0: + + for entry in annotations: + + if 'PICK': + + gene = entry['SYMBOL'] + hgvs_p = entry['HGVSp'] + hgvs_c = entry['HGVSc'] + consequence = entry['Consequence'] + exon = entry['EXON'] + + transcript_info = (gene + , hgvs_p + , hgvs_c + , consequence + , exon) + + transcripts.append(transcript_info) + + # If still empty leave empty + if len(transcripts) == 0: + + gene = "" + hgvs_p = "" + hgvs_c = "" + consequence = "" + exon = "" + + transcript_info = (gene + , hgvs_p + , hgvs_c + , consequence + , exon) + + transcripts.append(transcript.info) + + #Check if NTC variant list empty + if len(ntc_list) == 0: + in_ntc = "False" + ntc_vaf = "" + ntc_depth = "" + ntc_alt_reads = "" + else: + #Check if variant in NTC + for entry in ntc_list: + variant_full = chrom+":"+str(pos)+ref+">"+alt + if variant_full == entry[0]: + in_ntc = "True" + ntc_vaf = entry[1] + ntc_depth = entry[2] + ntc_alt_reads = entry[3] + else: + in_ntc = "False" + ntc_vaf = "" + ntc_depth = "" + ntc_alt_reads = "" + + #Go through all annotations we're keeping for this variant and add them + for keep in transcripts: + var = f'{keep[0]}\t{chrom}\t{pos}\t{ref}\t{alt}\t{vaf}\t{depth}\t{keep[1]}\t{keep[2]}\t{keep[3]}\t{keep[4]}\t{alt_reads}\t{in_ntc}\t{ntc_vaf}\t{ntc_depth}\t{ntc_alt_reads}' + variant_list.append(var) + + + return(variant_list) + + +## Input arguments +def get_args(): + """ + Uses argparse to take arguments from command line + """ + + parser = argparse.ArgumentParser( + formatter_class=argparse.RawTextHelpFormatter, + description=textwrap.dedent( + ''' + Summary: + Takes VCF file created by 1_SomaticAmplicon.sh and parses + the variants to produce a tab delimited variant report + ''' + )) + + # Required: VCF input file + parser.add_argument('-i', '--inp', action='store', + help ='Filepath to input VCF file. REQUIRED.' + ) + + # OPTIONAL: Output folder, defaults to current directory if empty + parser.add_argument( + '-o', '--outfile', action='store', + help=textwrap.dedent( + ''' + Filepath to folder where output reports will be saved. + If missing, defaults to current directory. + \n''' + )) + + # Required: Preferred transcript file + parser.add_argument( + '-t', '--pref_trans', action='store', + help='Filepath to preferred transcrips' + ) + + # Required: Psth to ntc vcf file + parser.add_argument( + '-n', '--ntc_vcf', action='store', + help='Filepath to NTC VCF' + ) + + return parser.parse_args() + +# Function to run main functions def main(args): - # setup logger - logger = logging.getLogger('vcf_parse') - logger.setLevel(logging.DEBUG) - handler = logging.StreamHandler() - handler.setLevel(logging.DEBUG) - formatter = logging.Formatter( - '%(levelname)s\t%(asctime)s\t%(name)s\t%(message)s' - ) - handler.setFormatter(formatter) - logger.addHandler(handler) - logger.info('running vcf_parse.py...') - - # Load arguments, make vcf report object and load data - report = vcf_report() - report.load_data(args.input, args.output) - - # If -l flag called, print headers and exit - if args.config_list: - report.list_config() - exit() - - # If config file provided, load config - if args.config: - report.load_config(args.config) - else: - logger.info('no config file found -- outputting all data from VCF.') - - # Make variant report of whole VCF - report.make_report(args.filter_non_pass) - - # If preferred transcripts provided, apply to variant report - if args.transcripts: - pt = preferred_transcripts() - pt.load(args.transcripts) - pt.apply(report, args.transcript_strictness) - else: - logger.info('no preferred transcripts file provided -- preferred ' + - 'transcripts column will all be labelled as "Unknown"') - - # If known variants provided, apply to variant report - if args.known_variants: - known = known_variants() - known.load_known_variants(args.known_variants) - known.apply_known_variants(report) - - else: - logger.info('no known variants file provided -- Classification ' + - 'column will be empty') - - # If single BED file provided, make variant report with BED file - # applied - if args.bed: - bed = bed_object() - bed.apply_single(args.bed, report) - - # If folder of BED file provided, make a seperate variant report - # for each BED file. Output will be saved in a folder named the - # same as the BED file folder, within the output directory. - elif args.bed_folder: - bed = bed_object() - bed.apply_multiple(args.bed_folder, report) - - # If no BED files provided, pass - else: - logger.info('no BED files provided') - - # Finish - logger.info('vcf_parse.py completed\n{}'.format('---'*30)) - - -# -- CALL FUNCTIONS --------------------------------------------------- + # Setup logger + logger = logging.getLogger('vcf_parse') + logger.setLevel(logging.DEBUG) + + handler = logging.StreamHandler() + handler.setLevel(logging.DEBUG) + + formatter = logging.Formatter( + '%(levelname)s\t%(asctime)s\t%(name)s\t%(message)s' + ) + + handler.setFormatter(formatter) + + logger.addHandler(handler) + logger.info('running vcf_parse.py...') + + #Load arguments, make vcf report object and load data + + # Parse preferred transcripts + if args.pref_trans: + pref_transc = parse_transcripts(args.pref_trans) + else: + logger.info('no preffered transcript file provided - please provide') + + # Parse NTC + if args.ntc_vcf: + ntc_list = parse_ntc(args.ntc_vcf) + else: + logger.info('no ntc vcf provided - please provide') + + # Parse sample + variant_list = parse_sample(args.inp, args.outfile, ntc_list, pref_transc) + + # Finish + logger.info('steps completed\n{}'.format('---'*30)) + + # Export data + # Header lines + logger.info('creating output file') + + outfile_headers = ['gene', 'chr', 'pos', 'ref', 'alt', 'vaf', 'depth', 'hgvs_p', 'hgvs_c', 'consequence', 'exon', 'alt_reads', 'in_ntc', 'ntc_vaf', 'ntc_depth', 'ntc_alt_reads'] + outfile_headers_string = '\t'.join(outfile_headers) + + #Info + with open(args.outfile,'w',newline='') as file: + file.write(outfile_headers_string + '\n') + + for line in variant_list: + file.write(str(line) + '\n') + + # Final file + logger.info('vcf_parse.py completed\n{}'.format('---'*30)) + +# Call functions if __name__ == '__main__': - args = get_args() - main(args) + args = get_args() + main(args)