#!/usr/bin/env python

import os,re,sys
import csv
import logging
import gzip
from cyvcf2 import VCF, Writer
import subprocess


csv.field_size_limit(500 * 1024 * 1024)
threeLettertoOneLetterAA = {'Ala':'A','Arg':'R','Asn':'N','Asp':'D','Cys':'C','Glu':'E','Gln':'Q','Gly':'G','His':'H','Ile':'I','Leu':'L','Lys':'K', 'Met':'M','Phe':'F','Pro':'P','Ser':'S','Thr':'T','Trp':'W','Tyr':'Y','Val':'V','Ter':'X'}


def read_infotag_file(vcf_info_tags_tsv):
   """
   Function that reads a file that lists VCF INFO tags produced by PCGR/CPSR/gvanno.
   An example line of the VCF info tag file is the following:
   
   tag	number	type	description category
   Consequence	.	String	"Impact modifier for the consequence type (picked by VEP's --flag_pick_allele option)."   vep
   
   A dictionary is returned, with the tag as the key, and the full dictionary record as the value
   """
   info_tag_xref = {} ##dictionary returned
   if not os.path.exists(vcf_info_tags_tsv):
      return info_tag_xref
   tsvfile = open(vcf_info_tags_tsv, 'r')
   reader = csv.DictReader(tsvfile, delimiter='\t')
   for row in reader:
      if not row['tag'] in info_tag_xref:
         info_tag_xref[row['tag']] = row
   
   return info_tag_xref

def check_subprocess(command):
   #if debug:
      #logger.info(command)
   try:
      output = subprocess.check_output(str(command), stderr=subprocess.STDOUT, shell=True)
      if len(output) > 0:
         print (str(output.decode()).rstrip())
   except subprocess.CalledProcessError as e:
      print (e.output.decode())
      exit(0)


def read_genexref_namemap(gene_xref_namemap_tsv):
  
   """
   Function that reads a file that lists names of tags in PCGR_ONCO_XREF annotation.
   """
  
   namemap_xref = {} ##dictionary returned
   if not os.path.exists(gene_xref_namemap_tsv):
      return namemap_xref
   tsvfile = open(gene_xref_namemap_tsv, 'r')
   reader = csv.DictReader(tsvfile, delimiter='\t')
   for row in reader:
      namemap_xref[row['name']] = int(row['index'])

   return namemap_xref

def is_integer(n):
    try:
        float(n)
    except ValueError:
        return False
    else:
        return float(n).is_integer()

def error_message(message, logger):
   logger.error('')
   logger.error(message)
   logger.error('')
   exit(1)

def write_pass_vcf(annotated_vcf, logger):
   """
   Function prints all PASS variants from a VCF file. New VCF file appends '.pass.' to the filename.
   
   """
   #out_vcf = re.sub(r'\.annotated\.vcf\.gz$','.annotated.pass.vcf',annotated_vcf)
   out_vcf = re.sub(r'\.vcf\.gz$','.pass.vcf',annotated_vcf)
   vcf = VCF(annotated_vcf)
   w = Writer(out_vcf, vcf)

   num_rejected = 0
   num_pass = 0
   for rec in vcf:
      if rec.FILTER is None or rec.FILTER == 'None':
         w.write_record(rec)
         num_pass += 1
      else:
         num_rejected +=1

   vcf.close()
   w.close()
   
   logger.info('Number of non-PASS/REJECTED variant calls: ' + str(num_rejected))
   logger.info('Number of PASSed variant calls: ' + str(num_pass))
   if num_pass == 0:
      logger.warning('There are zero variants with a \'PASS\' filter in the VCF file')
      os.system('bgzip -dc ' + str(annotated_vcf) + ' egrep \'^#\' > ' + str(out_vcf))
   #else:
   os.system('bgzip -f ' + str(out_vcf))
   os.system('tabix -f -p vcf ' + str(out_vcf) + '.gz')

   return

def warn_message(message, logger):
   logger.warning(message)

def getlogger(logger_name):
	logger = logging.getLogger(logger_name)
	logger.setLevel(logging.DEBUG)

	# create console handler and set level to debug
	ch = logging.StreamHandler(sys.stdout)
	ch.setLevel(logging.DEBUG)

	# add ch to logger
	logger.addHandler(ch)
	
	# create formatter
	formatter = logging.Formatter("%(asctime)s - %(name)s - %(levelname)s - %(message)s", "20%y-%m-%d %H:%M:%S")
	
	#add formatter to ch
	ch.setFormatter(formatter)
	
	return logger


def is_valid_vcf(input_vcf, output_dir, logger, debug):
   """
   Function that reads the output file of EBIvariation/vcf-validator and reports potential errors and validation status
   """

   logger.info('Validating VCF file with EBIvariation/vcf-validator (v0.9.3)')

   ## Filename for output file from VCF validation
   vcf_validation_output_file = os.path.join(output_dir, re.sub(r'(\.vcf$|\.vcf\.gz$)', '.vcf_validator_output', os.path.basename(input_vcf)))
   
   ## Command for running vcf-validator
   command_v42 = 'vcf_validator -i ' + str(input_vcf) + ' -l warning -r text -o ' + str(output_dir) + ' > ' + str(vcf_validation_output_file)  + ' 2>&1'
   if debug:
      logger.info(command_v42)
   try:
      subprocess.run(command_v42, stderr=subprocess.STDOUT, shell=True)
   except subprocess.CalledProcessError as e:
      print (e.output.decode())
      exit(0)

   ## read output file from vcf-validator and get name of log/error file
   validation_results = {}
   validation_results['validation_status'] = 0
   validation_results['error_messages'] = []
   validation_log_fname = None
   if os.path.exists(vcf_validation_output_file):
      f = open(vcf_validation_output_file, 'r')
      for line in f:
         if 'Text report written to' in line.rstrip():
            validation_log_fname = re.split(r' : ',line.rstrip(),maxsplit=2)[1]
      f.close()
      if not debug:
         check_subprocess('rm -f ' + str(vcf_validation_output_file))

   if validation_log_fname is None:
      err_msg = 'Cannot find file with error messages from vcf_validator'
      return error_message(err_msg, logger)

   if os.path.exists(validation_log_fname):
      f = open(validation_log_fname, 'r')
      for line in f:
         if not re.search(r' \(warning\)$|^Reading from ',line.rstrip()): ## ignore warnings
            if line.startswith('Line '):
               validation_results['error_messages'].append('ERROR: ' + line.rstrip())
            if 'the input file is valid' in line.rstrip(): ## valid VCF
               validation_results['validation_status'] = 1
            if 'the input file is not valid' in line.rstrip():  ## non-valid VCF
               validation_results['validation_status'] = 0
      f.close()
      if not debug:
         check_subprocess('rm -f ' + str(validation_log_fname))
   else:
      err_msg = str(validation_log_fname) + ' does not exist'
      return error_message(err_msg, logger)

   if validation_results['validation_status'] == 0:
      error_string_42 = '\n'.join(validation_results['error_messages'])
      validation_status = 'According to the VCF specification, the VCF file (' + str(input_vcf) + ') is NOT valid'
      err_msg = validation_status + ':\n' + str(error_string_42)
      return error_message(err_msg, logger)
   else:
      validation_status = 'According to the VCF specification, the VCF file ' + str(input_vcf) + ' is valid'
      logger.info(validation_status)
   return 0


def map_regulatory_variant_annotations(vep_csq_records):
   """
   Function that considers an array of VEP CSQ records and appends all regulatory variant consequent annotations (open chromatin, TF_binding_site, 
   CTCF_binding_site, promoter (flanks), enhancers ) into a single comma-separated string. Each individual regulatory annotation is formatted as:
   <Consequence>|<Feature_type>|<Feature>|<BIOTYPE>|<MOTIF_NAME>|<MOTIF_POS>|<HIGH_INF_POS>|<MOTIF_SCORE_CHANGE>|<TRANSCRIPTION_FACTORS>
   """

   regulatory_annotation = '.'
   if len(vep_csq_records) == 1:
      return regulatory_annotation

   j = 0
   regulatory_annotations = []
   while j < len(vep_csq_records):
      missing_key_annotation = False
      for k in ['Feature_type','Consequence','Feature']:
         if not k in vep_csq_records[j].keys():
            missing_key_annotation = True

      if missing_key_annotation is False:

         ## RegulatoryFeature annotations - open chromatin, promoters (flanks), enhancers, CTCF binding sites
         if vep_csq_records[j]['Feature_type'] == 'RegulatoryFeature':
            biotype = ""
            if re.match(r"^(enhancer|promoter|open|CTCF|TF_)", vep_csq_records[j]['BIOTYPE']):
               biotype = vep_csq_records[j]['BIOTYPE']
            
            annotation = str(vep_csq_records[j]['Consequence']) + '|' + \
               str(vep_csq_records[j]['Feature_type']) + '|' + \
               str(vep_csq_records[j]['Feature']) + '|' + \
               str(biotype) + '|||||'

            regulatory_annotations.append(annotation)

         ## MotifFeature annotations (TF)
         if vep_csq_records[j]['Feature_type'] == 'MotifFeature':
            missing_motif_annotation = False
            for annotation in ['MOTIF_NAME','MOTIF_POS','HIGH_INF_POS','MOTIF_SCORE_CHANGE','TRANSCRIPTION_FACTORS']:
               if not annotation in vep_csq_records[j].keys():
                  missing_motif_annotation = True
            
            if missing_motif_annotation is False:
               annotation = str(vep_csq_records[j]['Consequence']) + '|' + \
                  str(vep_csq_records[j]['Feature_type']) + '|' + \
                  str(vep_csq_records[j]['Feature']) + '|TF_binding_site|' + \
                  str(vep_csq_records[j]['MOTIF_NAME']) + '|' + \
                  str(vep_csq_records[j]['MOTIF_POS']) + '|' + \
                  str(vep_csq_records[j]['HIGH_INF_POS']) + '|' + \
                  str(vep_csq_records[j]['MOTIF_SCORE_CHANGE']) + '|' + \
                  str(vep_csq_records[j]['TRANSCRIPTION_FACTORS'])

               regulatory_annotations.append(annotation)
      
      j = j + 1
   
   if len(regulatory_annotations) > 0:
      regulatory_annotation = ','.join(regulatory_annotations)

   return(regulatory_annotation)
               
def get_correct_cpg_transcript(vep_csq_records):
   """
   Function that considers an array of VEP CSQ records and picks most relevant consequence (and gene) from
   neighbouring genes/transcripts of relevance for cancer predisposition (cpg = cancer predisposition gene)
   """

  
   csq_idx = 0
   if len(vep_csq_records) == 1:
      return csq_idx
   
   
   ## some variants iare assigned multiple transcript consequences
   ## if cancer predisposition genes are in the vicinity of other genes, choose the cancer predisposition gene
   ## if there are neighbouring cancer-predispositon genes, choose custom gene, preferring coding change (see below, KLLN/PTEN, XPC/TMEM43, NTHL1/TSC2)
   csq_idx_dict = {}
   for g in ['KLLN','PTEN','XPC','TMEM43','NTHL1','TSC2']:
      csq_idx_dict[g] = {}
      csq_idx_dict[g]['idx'] = -1
      csq_idx_dict[g]['coding'] = False

   j = 0
   while j < len(vep_csq_records):
      if 'CANCER_PREDISPOSITION_SOURCE' in vep_csq_records[j].keys() or 'GE_PANEL_ID' in vep_csq_records[j].keys():
         csq_idx = j
         if 'SYMBOL' in vep_csq_records[j].keys():
            if vep_csq_records[j]['SYMBOL'] in csq_idx_dict.keys():
               csq_idx_dict[str(vep_csq_records[j]['SYMBOL'])]['idx'] = j
               if vep_csq_records[j]['CODING_STATUS'] == 'coding':
                  csq_idx = j  # prefer coding on over anything else
                  csq_idx_dict[str(vep_csq_records[j]['SYMBOL'])]['coding'] = True
      j = j + 1
   
   if csq_idx_dict['KLLN']['idx'] != -1 and csq_idx_dict['PTEN']['idx'] != -1:
      csq_idx = csq_idx_dict['PTEN']['idx']
      if csq_idx_dict['KLLN']['coding'] is True:
         csq_idx = csq_idx_dict['KLLN']['idx']
   
   if csq_idx_dict['XPC']['idx'] != -1 and csq_idx_dict['TMEM43']['idx'] != -1:
      csq_idx = csq_idx_dict['XPC']['idx']
      if csq_idx_dict['TMEM43']['coding'] is True:
         csq_idx = csq_idx_dict['TMEM43']['idx']
   
   if csq_idx_dict['TSC2']['idx'] != -1 and csq_idx_dict['NTHL1']['idx'] != -1:
      csq_idx = csq_idx_dict['TSC2']['idx']
      if csq_idx_dict['NTHL1']['coding'] is True:
         csq_idx = csq_idx_dict['NTHL1']['idx']

   if csq_idx is None:
      csq_idx = 0
   return csq_idx


def threeToOneAA(aa_change):
	
   for three_letter_aa in threeLettertoOneLetterAA.keys():
      aa_change = aa_change.replace(three_letter_aa,threeLettertoOneLetterAA[three_letter_aa])

   aa_change = re.sub(r'[A-Z]{1}fsX([0-9]{1,}|\?)','fs',aa_change)
   return aa_change

def map_variant_effect_predictors(rec, algorithms):
   
   dbnsfp_predictions = map_dbnsfp_predictions(str(rec.INFO.get('DBNSFP')), algorithms)
   if rec.INFO.get('Gene') is None or rec.INFO.get('Consequence') is None:
      return
   gene_id = str(rec.INFO.get('Gene'))
   consequence = str(rec.INFO.get('Consequence'))
     
   dbnsfp_key = ''

   found_key = 0
   if not rec.INFO.get('HGVSp_short') is None and not rec.INFO.get('HGVSp_short') == '.':
      aa_change = str(rec.INFO.get('HGVSp_short'))
      dbnsfp_key = gene_id + ':' + str(aa_change)
      if dbnsfp_key in dbnsfp_predictions:
         found_key = 1
   
   if found_key == 0 and re.search('splice_',consequence):
      dbnsfp_key = gene_id

   if dbnsfp_key != '':
      if dbnsfp_key in dbnsfp_predictions:
         rec.INFO['EFFECT_PREDICTIONS'] = dbnsfp_predictions[dbnsfp_key]
         for algo_pred in rec.INFO['EFFECT_PREDICTIONS'].split('&'):
            if algo_pred.startswith('sift:'):
               rec.INFO['DBNSFP_SIFT'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('provean:'):
               rec.INFO['DBNSFP_PROVEAN'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('m-cap:'):
               rec.INFO['DBNSFP_M_CAP'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('mutpred:'):
               rec.INFO['DBNSFP_MUTPRED'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('metarnn:'):
               rec.INFO['DBNSFP_META_RNN'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('fathmm:'):
               rec.INFO['DBNSFP_FATHMM'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('fathmm_mkl_coding:'):
               rec.INFO['DBNSFP_FATHMM_MKL'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('mutationtaster:'):
               rec.INFO['DBNSFP_MUTATIONTASTER'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('mutationassessor:'):
               rec.INFO['DBNSFP_MUTATIONASSESSOR'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('deogen2:'):
               rec.INFO['DBNSFP_DEOGEN2'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('primateai:'):
               rec.INFO['DBNSFP_PRIMATEAI'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('list_s2:'):
               rec.INFO['DBNSFP_LIST_S2'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('gerp_rs:'):
               rec.INFO['DBNSFP_GERP'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('bayesdel_addaf:'):
               rec.INFO['DBNSFP_BAYESDEL_ADDAF'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('aloft:'):
               rec.INFO['DBNSFP_ALOFTPRED'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('splice_site_rf:'):
               rec.INFO['DBNSFP_SPLICE_SITE_RF'] = str(algo_pred.split(':')[1])
            if algo_pred.startswith('splice_site_ada:'):
               rec.INFO['DBNSFP_SPLICE_SITE_ADA'] = str(algo_pred.split(':')[1])



def detect_reserved_info_tag(tag, tag_name, logger):
   reserved_tags = ['AA','AC','AF','AN','BQ','CIGAR','DB','DP','END','H2','H3','MQ','MQ0','NS','SB','SOMATIC','VALIDATED','1000G']
   if tag in reserved_tags:
      err_msg = 'Custom INFO tag (' + str(tag_name) + ') needs another name - ' + str(tag) + ' is a reserved field in the VCF specification (INFO)'
      return error_message(err_msg, logger)
   
   reserved_format_tags = ['GT','DP','FT','GL','GLE','GQ','PL','HQ','PS','PQ','EC','MQ']
   if tag in reserved_format_tags:
      err_msg = 'Custom INFO tag (' + str(tag_name) + ') needs another name - ' + str(tag) + ' is a reserved field in the VCF specification (FORMAT)'
      return error_message(err_msg, logger)


def assign_cds_exon_intron_annotations(csq_record):
   csq_record['CODING_STATUS'] = 'noncoding'
   csq_record['EXONIC_STATUS'] = 'nonexonic'
   csq_record['SPLICE_DONOR_RELEVANT'] = False
   csq_record['NULL_VARIANT'] = False
   csq_record['INTRON_POSITION'] = 0
   csq_record['EXON_POSITION'] = 0
   
   coding_csq_pattern = r"^(stop_|start_lost|frameshift_|missense_|splice_donor|splice_acceptor|protein_altering|inframe_)"
   wes_csq_pattern = r"^(stop_|start_lost|frameshift_|missense_|splice_donor|splice_acceptor|inframe_|protein_altering|synonymous)"
   null_pattern = r"^(stop_|frameshift_)"
   if re.match(coding_csq_pattern, str(csq_record['Consequence'])):
      csq_record['CODING_STATUS'] = 'coding'
   
   if re.match(wes_csq_pattern, str(csq_record['Consequence'])):
      csq_record['EXONIC_STATUS'] = 'exonic'
   
   if re.match(null_pattern, str(csq_record['Consequence'])):
      csq_record['NULL_VARIANT'] = True

   if re.match(r"^splice_region_variant", str(csq_record['Consequence'])) and re.search(r'(\+3(A|G)>|\+4G>|\+5G>)', str(csq_record['HGVSc'])):
      csq_record['SPLICE_DONOR_RELEVANT'] = True

   if re.match(r"splice_region_variant|intron_variant", str(csq_record['Consequence'])):
      match = re.search(r"((-|\+)[0-9]{1,}(dup|del|inv|((ins|del|dup|inv|delins)(A|G|C|T){1,})|(A|C|T|G){1,}>(A|G|C|T){1,}))$", str(csq_record['HGVSc']))
      if match is not None:
         pos = re.sub(r"(\+|dup|del|delins|ins|inv|(A|G|C|T){1,}|>)","",match.group(0))
         if is_integer(pos):
            csq_record['INTRON_POSITION'] = int(pos)
   
   if 'NearestExonJB' in csq_record.keys():
      if not csq_record['NearestExonJB'] is None:
         if re.match(r"synonymous_|missense_|stop_|inframe_|start_", str(csq_record['Consequence'])) and str(csq_record['NearestExonJB']) != "":
            exon_pos_info = csq_record['NearestExonJB'].split("+")
            if len(exon_pos_info) == 4:
               if is_integer(exon_pos_info[1]) and str(exon_pos_info[2]) == "end":
                  csq_record['EXON_POSITION'] = -int(exon_pos_info[1])
               if is_integer(exon_pos_info[1]) and str(exon_pos_info[2]) == "start":
                  csq_record['EXON_POSITION'] = int(exon_pos_info[1])


   for m in ['HGVSp_short','CDS_CHANGE']:
      csq_record[m] = '.'
   if not csq_record['HGVSc'] is None:
      if csq_record['HGVSc'] != '.':
         if 'splice_acceptor_variant' in csq_record['Consequence'] or 'splice_donor_variant' in csq_record['Consequence']:
            key = str(csq_record['Consequence']) + ':' + str(csq_record['HGVSc'])
            csq_record['CDS_CHANGE'] = key
   if csq_record['Amino_acids'] is None or csq_record['Protein_position'] is None or csq_record['Consequence'] is None:
      return
   if not csq_record['Protein_position'] is None:
      if csq_record['Protein_position'].startswith('-'):
         return

   protein_change = '.'
   if '/' in csq_record['Protein_position']:
      protein_position = str(csq_record['Protein_position'].split('/')[0])
      if '-' in protein_position:
         if protein_position.split('-')[0].isdigit():
            csq_record['AMINO_ACID_START'] = protein_position.split('-')[0]
         if protein_position.split('-')[1].isdigit():
            csq_record['AMINO_ACID_END'] = protein_position.split('-')[1]
      else:
         if protein_position.isdigit():
            csq_record['AMINO_ACID_START'] = protein_position
            csq_record['AMINO_ACID_END'] = protein_position
   
   if not csq_record['HGVSp'] is None:
      if csq_record['HGVSp'] != '.':
         if ':' in csq_record['HGVSp']:
            protein_identifier = str(csq_record['HGVSp'].split(':')[0])
            if protein_identifier.startswith('ENSP'):
               protein_change_VEP = str(csq_record['HGVSp'].split(':')[1])
               protein_change = threeToOneAA(protein_change_VEP)
  
   if 'synonymous_variant' in csq_record['Consequence']:
      protein_change = 'p.' + str(csq_record['Amino_acids']) + str(protein_position) + str(csq_record['Amino_acids'])
      if 'stop_lost' in str(csq_record['Consequence']) and '/' in str(csq_record['Amino_acids']):
         protein_change = 'p.X' + str(protein_position) + str(csq_record['Amino_acids']).split('/')[1]
    
   csq_record['HGVSp_short'] = protein_change
   exon_number = 'NA'
   if not csq_record['EXON'] is None:
      if csq_record['EXON'] != '.':
         if '/' in csq_record['EXON']:
            exon_number = str(csq_record['EXON']).split('/')[0]
            num_exons = str(csq_record['EXON']).split('/')[1]
            if exon_number == num_exons:
               csq_record['LAST_EXON'] = True

   if not csq_record['INTRON'] is None:
      if csq_record['INTRON'] != '.':
         if '/' in csq_record['INTRON']:
            intron_number = str(csq_record['INTRON']).split('/')[0]
            num_introns = str(csq_record['INTRON']).split('/')[1]
            if intron_number == num_introns:
               csq_record['LAST_INTRON'] = True

   if not csq_record['HGVSc'] is None:
      if csq_record['HGVSc'] != '.':
         if protein_change != '.':
            key = str(csq_record['Consequence']) + ':' + str(csq_record['HGVSc']) + ':exon' + str(exon_number) + ':' + str(protein_change)
            csq_record['CDS_CHANGE'] = key

   return

def map_dbnsfp_predictions(dbnsfp_tag, algorithms):
   
   effect_predictions = {}
   
   for v in dbnsfp_tag.split(','):
   
      dbnsfp_info = v.split('|')
      if len(dbnsfp_info) == 1:
         return effect_predictions
      ref_aa = dbnsfp_info[0]
      alt_aa = dbnsfp_info[1]
      all_ids = dbnsfp_info[4].split('&')
      unique_ids = {}
      for s in all_ids:
         unique_ids[s] = 1
         
      isoform_aa_keys = []
      if ref_aa != '.' and alt_aa != '.' and ref_aa != '' and alt_aa != '':
         aa_pos = dbnsfp_info[6].split('&')
         for pos in aa_pos:
            for gene_id in unique_ids:
               k = str(gene_id) + ':p.' + str(ref_aa) + pos + str(alt_aa)
               isoform_aa_keys.append(k)
      else:
         #continue
         for gene_id in unique_ids:
            isoform_aa_keys.append(gene_id)
   
      algorithm_raw_predictions = {}
   
      i = 7
      v = 0
      
      if len(algorithms) != len(dbnsfp_info[7:]):
         return effect_predictions
      
      while i < len(dbnsfp_info):
         algorithm_raw_predictions[str(algorithms[v]).lower()] = dbnsfp_info[i].split('&')
         i = i + 1
         v = v + 1
      dbnsfp_predictions = {}
      
      for k in isoform_aa_keys:
         if not k in dbnsfp_predictions:
            dbnsfp_predictions[k] = {}
         all_preds = []
         for algo in algorithm_raw_predictions.keys():
            unique_algo_predictions = {}
            for pred in algorithm_raw_predictions[algo]: