import argparse
from Bio import SeqIO
import json

if __name__ == '__main__':
    parser = argparse.ArgumentParser(
        description="Append newly sequenced genomes to current genome dataset, and export metadata",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter
    )
    parser.add_argument("--genomes", required=True, help="FASTA or JSON file with latest genomes from GISAID")
    parser.add_argument("--new-genomes", required=True, help="FASTA file with newly sequenced genomes")
    parser.add_argument("--keep", required=True, help="TXT file with accession number of genomes to be included")
    parser.add_argument("--remove", required=True, help="TXT file with accession number of genomes to be removed")
    parser.add_argument("--output", required=True, help="FASTA file containing filtered sequences")
    args = parser.parse_args()

    genomes = args.genomes
    new_genomes = args.new_genomes
    keep = args.keep
    remove = args.remove
    outfile = args.output

    # genomes = path + "pre-analyses/provision.json"
    # # genomes = path + "pre-analyses/gisaid_hcov-19.fasta"
    # new_genomes = path + "pre-analyses/new_genomes.fasta"
    # keep = path + 'config/keep.txt'
    # remove = path + "config/remove.txt"
    # outfile = path + "pre-analyses/temp_sequences.fasta"

    # inspect genome coverage
    genome_size = 29903
    max_gaps = 30
    min_size = genome_size - int(genome_size * (max_gaps / 100))

    # store only new sequences in a dictionary, ignoring existing ones
    print('\n### Loading new sequences, and reporting genome coverage\n')
    newly_sequenced = {}
    low_coverage = {}
    for fasta in SeqIO.parse(open(new_genomes), 'fasta'):
        id, seq = fasta.description, fasta.seq
        size = len(str(seq).replace('N', '').replace('-', ''))
        if size > min_size:
            coverage = str(round(size / genome_size, 3))
            print(id + ', coverage = ' + coverage + ' (PASS)')
            if id not in newly_sequenced.keys():  # avoid potential duplicates
                newly_sequenced[id] = str(seq)
        else:
            coverage = str(round(size / genome_size, 3))
            print(id + ', coverage = ' + coverage + ' (FAIL)')
            low_coverage[id] = coverage
    print('\nDone!\n')

    # create a list of sequences to be added in all instances
    keep_sequences = []
    for id in open(keep, "r").readlines():
        if id[0] not in ["#", "\n"]:
            id = id.strip().replace('hCoV-19/', '')
            if id not in newly_sequenced:
                if id not in keep_sequences:
                    keep_sequences.append(id)

    # create a list of sequences to be ignored in all instances
    remove_sequences = []
    for id in open(remove, "r").readlines():
        if id[0] not in ["#", "\n"]:
            id = id.strip()
            remove_sequences.append(id)

    # export only sequences to be used in the nextstrain build
    c = 1
    print('\n### Exporting sequences\n')
    exported = []
    ignored = []
    all_sequences = []
    with open(outfile, 'w') as output:
        if genomes.split('.')[1] == 'json':
            with open(genomes) as infile:
                c = 0
                for line in infile:
                    try:
                        entry = json.loads(line)
                        id = entry['covv_virus_name']
                        id = id.split('|')[0].replace('hCoV-19/', '')
                        if len(id.split('/')) == 3:
                            country, index, year = id.split('/')
                        elif len(id.split('/')) == 4:
                            host, country, index, year = id.split('/')
                        country = country.replace(' ', '').replace('\'', '-').replace('_', '')
                        id = '/'.join([country, index, year])

                        all_sequences.append(id)
                        if id not in remove_sequences:
                            if id in keep_sequences:  # filter out unwanted sequences
                                seq = entry['sequence'].replace('\n', '')
                                if int(len(seq)) >= min_size:
                                    entry = ">" + id + "\n" + seq.upper() + "\n"
                                    exported.append(id)
                                    output.write(entry)
                                    print(str(c) + '. ' + id)
                                    c += 1
                        else:
                            ignored.append(id)
                    except:
                        print('Corrupted JSON row detected. Skipping...')

        else:
            for fasta in SeqIO.parse(open(genomes), 'fasta'):
                id, seq = fasta.description, fasta.seq
                id = id.split('|')[0].replace('hCoV-19/', '')
                if len(id.split('/')) == 3:
                    country, index, year = id.split('/')
                elif len(id.split('/')) == 4:
                    host, country, index, year = id.split('/')
                country = country.replace(' ', '').replace('\'', '-').replace('_', '')
                id = '/'.join([country, index, year])
                all_sequences.append(id)

                if id not in remove_sequences:
                    if id in keep_sequences:  # filter out unwanted sequences
                        entry = ">" + id + "\n" + str(seq).upper() + "\n"
                        exported.append(id)
                        output.write(entry)
                        print(str(c) + '. ' + id)
                        c += 1
                else:
                    ignored.append(id)

        for id, seq in newly_sequenced.items():
            print('* ' + str(c) + '. ' + id)
            entry = ">" + id + "\n" + seq.upper() + "\n"
            exported.append(id)
            output.write(entry)
            c += 1
    print('\n- Done!\n')

    # mismatched sequence headers
    mismatch = [genome for genome in keep_sequences if genome not in all_sequences]
    if len(mismatch) + len(low_coverage) > 0:
        print('\n### WARNINGS!')

    if len(mismatch) > 0:
        print('\n### Possible sequence header mismatches\n')
        m = 1
        for id in mismatch:
            print(str(m) + '. ' + id)
            m += 1
    else:
        print('\nNo sequence name mismatches found...')

    if len(low_coverage) > 0:
        print('\n- Low quality sequences were ignored.\n')
        l = 1
        for id, coverage in low_coverage.items():
            print('\t' + str(l) + '. ' + id + ', coverage = ' + coverage + ' (FAIL)')
            l += 1

    print('\n\n\n### Final result\n')

    print('Lab file contains ' + str(len(newly_sequenced)) + ' high coverage sequences')
    print('Lab file contains ' + str(len(low_coverage)) + ' low coverage sequences, which were ignored')
    print('GISAID file contains ' + str(len(all_sequences)) + ' sequences\n')

    print(str(len(mismatch)) + ' genomes in keep.txt were NOT FOUND on GISAID database')
    print(str(len(keep_sequences)) + ' genomes ADDED from GISAID dataset')
    print(str(len(newly_sequenced)) + ' newly sequenced genomes were added')
    print(str(len(low_coverage)) + ' low coverage genomes were ignored')
    print(str(len(ignored)) + ' genomes were REMOVED according to remove.txt\n')
    print(str(len(exported)) + ' genomes included in FINAL dataset\n')

import pycountry_convert as pyCountry
import pandas as pd
import argparse
from bs4 import BeautifulSoup as BS
import pycountry
from matplotlib import cm
import numpy as np


if __name__ == '__main__':
    parser = argparse.ArgumentParser(
        description="Generate ordered colour file for nextstrain build",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter
    )
    parser.add_argument("--metadata", required=True, help="Reformatted nextstrain metadata file")
    parser.add_argument("--coordinates", required=True,  help="TSV coordinates file being used in the build")
    parser.add_argument("--geoscheme", required=True, help="XML file with geographic scheme")
    parser.add_argument("--grid", required=True, help="HTML file with HEX colour matrices")
    parser.add_argument("--columns", nargs='+', type=str,   help="list of columns with geographic information")
    parser.add_argument("--output", required=True, help="TSV file containing ordered HEX colours based on locations")
    parser.add_argument("--filter", required=False, nargs='+', type=str,  help="List of filters for tagged rows in lab metadata")
    args = parser.parse_args()

    metadata = args.metadata
    geoscheme = args.geoscheme
    coordinates = args.coordinates
    grid = args.grid
    columns = args.columns
    output = args.output
    filt = args.filter

    # path = '/Users/anderson/GLab Dropbox/Anderson Brito/projects/ncov/ncov_variants/nextstrain/runX_ncov_20210407_vocvoicolours/'
    # metadata = path + 'data/metadata.tsv'
    # coordinates = path + 'config/latlongs.tsv'
    # geoscheme = path + 'config/geoscheme.tsv'
    # grid = path + 'config/colour_grid.html'
    # columns = ['region', 'country', 'division', 'location']
    # output = path + 'data/colors.tsv'


    # pre-determined HEX colours and hues
    print(filt)
    if filt==['connecticut']:
        force_colour = {'Connecticut': '#8FEECB', 'New York': '#19ae77', 'Canada': '#663300'}
        force_hue = {'North America': 0}
    if filt==['caribbean']:
        force_colour = {'Puerto Rico': '#8FEECB', 'US Virgin Islands': '#19ae77', 'Dominican Republic': '#663300'}
        force_hue = {'Caribbean': 0}

    # content to be exported as final result
    latlongs = {trait: {} for trait in columns}

    # get ISO alpha3 country codes
    isos =  {'US Virgin Islands':'VIR','Virgin Islands':'VIR','British Virgin Islands':'VGB','Curacao':'CUW','Northern Mariana Islands':'MNP',
            'Sint Maarten':'MAF','St Eustatius':'BES'}
    def get_iso(country):
        global isos
        if country not in isos.keys():
            try:
                isoCode = pyCountry.country_name_to_country_alpha3(country, cn_name_format="default")
                isos[country] = isoCode
            except:
                try:
                    isoCode = pycountry.countries.search_fuzzy(country)[0].alpha_3
                    isos[country] = isoCode
                except:
                    isos[country] = ''
        return isos[country]


    # extract coordinates from latlongs file for sorting places by latitude
    for line in open(coordinates).readlines():
        if not line.startswith('\n'):
            try:
                trait, place, lat, long = line.strip().split('\t')
                if trait in latlongs.keys():
                    entry = {place: (str(lat), str(long))}
                    latlongs[trait].update(entry) # save as pre-existing result
                else:
                    print('### WARNING! ' + trait + ' is not among the pre-selected columns with geographic information!')
            except:
                pass



    ''' REORDER LOCATIONS FOR LEGEND FORMATTING '''

    # open metadata file as dataframe
    dfN = pd.read_csv(metadata, encoding='utf-8', sep='\t', dtype=str)
    df = dfN
    dfN = dfN[['region', 'country', 'division', 'location']]

    ordered_regions = {}
    dcountries = {}
    places = []
    pinpoints = [dfN[trait].values.tolist() for trait in columns]
    for region_index in latlongs[columns[0]]:
        for address in zip(*pinpoints):
            address = list(address)
            region, country, division, location = [trait for trait in columns]
            region_name = address[0]
            country_name = address[1]
            division_name = address[2]
            location_name = address[3]
            places.append(address)

            if region_index == region_name:
                if 'subcontinent' not in ordered_regions.keys():
                    ordered_regions['subcontinent'] = {}
                if region_name not in ordered_regions['subcontinent'].keys():
                    ordered_regions['subcontinent'].update({region_name: latlongs[region][region_name]})

                ### STORE COUNTRIES WITH REGIONS AS KEYS
                if country_name in latlongs[country].keys():
                    if region_name not in dcountries.keys():
                        dcountries[region_name] = {country_name: latlongs[country][country_name]}
                    else:
                        if country_name not in dcountries[region_name].keys():
                            dcountries[region_name].update({country_name: latlongs[country][country_name]})



    # sort division entries based on sorted country entries
    ordered_countries = {}
    for region, countries in dcountries.items():
        ordered_countries[region] = {k: v for k, v in sorted(countries.items(), key=lambda item: item[1])}

    ddivisions = {}
    for country_index in [key for dict_ in ordered_countries.values() for key in dict_]:
        for address in places:
            region, country, division, location = [trait for trait in columns]
            country_name = address[1]
            division_name = address[2]
            location_name = address[3]

            if country_index == country_name:
                ### STORE DIVISIONS WITH COUNTRIES AS KEYS
                if division_name in latlongs[division].keys():
                    if country_name not in ddivisions.keys():
                        ddivisions[country_name] = {division_name: latlongs[division][division_name]}
                    else:
                        if division_name not in ddivisions[country_name].keys():
                            ddivisions[country_name].update({division_name: latlongs[division][division_name]})


    # sort division entries based on sorted country entries
    ordered_divisions = {}
    for country, divisions in ddivisions.items():
        ordered_divisions[country] = {k: v for k, v in sorted(divisions.items(), key=lambda item: item[1])}



    dlocations = {}
    for division_index in [key for dict_ in ordered_divisions.values() for key in dict_]:
        for address in places:
            address = list(address)
            region, country, division, location = [trait for trait in columns]
            division_name = address[2]
            location_name = address[3]

            if division_index == division_name:
                ### STORE LOCATIONS WITH DIVISIONS AS KEYS
                if location_name in latlongs[location].keys():
                    if division_name not in dlocations.keys():
                        dlocations[division_name] = {location_name: latlongs[location][location_name]}
                    else:
                        if location_name not in dlocations[division_name].keys():
                            dlocations[division_name].update({location_name: latlongs[location][location_name]})

    # sort locations entries based on sorted division entries
    ordered_locations = {}
    for division, locations in dlocations.items():
        ordered_locations[division] = {k: v for k, v in sorted(locations.items(), key=lambda item: item[1])}



    ''' CONVERT RGB TO HEX COLOURS '''

    # original source: https://bsou.io/posts/color-gradients-with-python

    # convert colour codes
    def RGB_to_hex(RGB):
        ''' [255,255,255] -> "#FFFFFF" '''
        # Components need to be integers for hex to make sense
        RGB = [int(x) for x in RGB]
        return "#" + "".join(["0{0:x}".format(v) if v < 16 else
                              "{0:x}".format(v) for v in RGB])

    def hex_to_RGB(hex):
        ''' "#FFFFFF" -> [255,255,255] '''
        # Pass 16 to the integer function for change of base
        return [int(hex[i:i + 2], 16) for i in range(1, 6, 2)]


    def color_dict(gradient):
        ''' Takes in a list of RGB sub-lists and returns dictionary of
          colors in RGB and hex form for use in a graphing function
          defined later on '''
        return [RGB_to_hex(RGB) for RGB in gradient]

    # create gradient
    def linear_gradient(start_hex, finish_hex, n):
        ''' returns a gradient list of (n) colors between
          two hex colors. start_hex and finish_hex
          should be the full six-digit color string,
          ("#FFFFFF") '''
        # Starting and ending colors in RGB form
        s = hex_to_RGB(start_hex)
        f = hex_to_RGB(finish_hex)
        # Initilize a list of the output colors with the starting color
        RGB_list = [s]
        # Calcuate a color at each evenly spaced value of t from 1 to n
        n = n
        for t in range(1, n):
            # Interpolate RGB vector for color at the current value of t
            curr_vector = [
                int(s[j] + (float(t) / (n - 1)) * (f[j] - s[j]))
                for j in range(3)]
            # Add it to list of output colors
            RGB_list.append(curr_vector)
        RGB_list = RGB_list
        return color_dict(RGB_list)




    ''' IMPORT GEOSCHEME '''

    # xml = BS(open(geoscheme, "r").read(), 'xml')
    # levels = xml.find('levels')

    scheme_list = open(geoscheme, "r").readlines()[1:]
    sampled_region = [key for dict_ in ordered_regions.values() for key in dict_]
    geodata = {}
    for line in scheme_list:
        if not line.startswith('\n'):
            type = line.split('\t')[0]
            if type == 'region':
                continent = line.split('\t')[1]
                region = line.split('\t')[2]
                if region in sampled_region:
                    if continent not in geodata.keys():
                        geodata[continent] = [region]
                    else:
                        geodata[continent] += [region]


    ''' IMPORT COLOUR SCHEME '''

    print('\nGenerating colour scheme...\n')
    html = BS(open(grid, "r").read(), 'html.parser')
    tables = html.find_all('table')

    # for colour_name in colour_scale.keys():
    limits = {'dark': (60, 30), 'light': (100, 80)}  # define saturation and luminance, max and min, respectively
    hue_to_hex = {}
    for table in tables:
        string_head = str(table.caption)
        if string_head == 'None':
            continue
        else:
            hue_value = int(string_head.split('"')[1])

        lux = table.tbody.find_all('tr')
        lum_value = 90  # brightest colour
        hexdark = ''
        hexligth = ''
        for row in lux:
            sat_value = 10  # unsaturated colour
            for cell in row.find_all('td'):
                # pass
                if sat_value == limits['dark'][0] and lum_value == limits['dark'][1]:
                    hexdark = cell.text.strip()
                if sat_value == limits['light'][0] and lum_value == limits['light'][1]:
                    hexligth = cell.text.strip()
                hex = (hexdark, hexligth)
                sat_value += 10
            lum_value -= 10
        hex = (hexdark, hexligth)
        hue_to_hex[hue_value] = hex


    colour_scale = {'magenta': [320], 'purple': [310, 300, 290, 280, 270, 260],
                    'blue': [250, 240, 230, 220], 'cyan': [210, 200, 190, 180], 'turquoise': [170, 160, 150],
                    'green': [140, 130, 120], 'yellowgreen': [110, 100, 90, 80, 70],
                    'yellow': [60, 50, 40], 'orange': [30, 20], 'red': [10, 0]}

    continent_hues = {'Oceania': colour_scale['magenta'], 'Asia': colour_scale['purple'],
                      'Europe': colour_scale['blue'] + colour_scale['cyan'], 'Africa': colour_scale['yellowgreen'],
                      'America': colour_scale['yellow'] + colour_scale['orange'] + colour_scale['red']}

    colour_wheel = {}
    palette = {}
    for area, subareas in geodata.items():
        num_subareas = len(subareas)
        hues = len(continent_hues[area])
        print(area, subareas)
        for position, subarea in zip([int(x) for x in np.linspace(0, int(hues), num_subareas, endpoint=False)], subareas):
            if subarea not in palette.keys():
                hue = continent_hues[area][position] # colour picker
                palette[subarea] = hue
                # print(subarea, hue)
    print(colour_wheel)
    # print(sampled_region)
    # print(palette)
    for region in sampled_region:
        if region in force_hue:
            colour_wheel[region] = hue_to_hex[force_hue[region]]
        else:
            colour_wheel[region] = hue_to_hex[palette[region]]


    ''' SET COLOUR SCHEME FOR UPDATES '''

    # convert a hue value into an rgb colour, and then hex colour
    def hue_to_rgb(hue):
        colour = int(hue / 240 * 255)
        # print(colour)
        luminance = 0.7
        rgb = [c * 255 * luminance for c in list(cm.jet(colour))[:-1]]
        # print(rgb)
        return RGB_to_hex(rgb)

    results = {trait: {} for trait in columns}


    ''' APPLY SAME HUE FOR MEMBERS OF THE SAME SUB-CONTINENT '''

    # assign countries to regions
    country_colours = {}
    reference_countries = {}
    for region, members in ordered_countries.items():
        countries = list(members.keys())
        for country in countries:
            reference_countries[country] = region
            country_colours[region] = countries


    # assign divisions to countries
    division_colours = {}
    reference_divisions = {}
    for country, members in ordered_divisions.items():
        divisions = list(members.keys())
        for division in divisions:
            reference_divisions[division] = reference_countries[country]
            if reference_countries[country] not in division_colours.keys():
                division_colours[reference_countries[country]] = divisions
            else:
                if division not in division_colours[reference_countries[country]]:
                    division_colours[reference_countries[country]].append(division)


    # assign locations to divisions
    location_colours = {}
    reference_locations = {}
    for division, members in ordered_locations.items():
        locations = list(members.keys())
        for location in locations:
            if reference_divisions[division] not in location_colours.keys():
                location_colours[reference_divisions[division]] = locations
            else:
                if location not in location_colours[reference_divisions[division]]:
                    location_colours[reference_divisions[division]].append(location)



    ''' CREATE COLOUR GRADIENT '''
    # define gradients for regions
    for continent, regions in geodata.items():
        hex_limits = []
        for region in sampled_region:
            if region in regions:
                if continent == 'America':
                    hex_limits += [list(colour_wheel[region])[0]]
                else:
                    hex_limits += list(colour_wheel[region])

        start, end = hex_limits[0], hex_limits[-1]
        if len(regions) == 1:
            gradient = linear_gradient(start, end, 4)
            gradient = [list(gradient)[2]]
        else:
            gradient = linear_gradient(start, end, len(regions))

        for region, colour in zip(regions, gradient):
            print('region', region, colour)
            results['region'].update({region: colour})

    # define gradients for country
    for hue, countries in country_colours.items():
        start, end = colour_wheel[hue]
        if len(countries) == 1:
            gradient = linear_gradient(start, end, 4)
            gradient = [list(gradient)[2]]
        else:
            gradient = linear_gradient(start, end, len(countries))
        for country, colour in zip(countries, gradient):
            print('country', country, colour)
            results['country'].update({country: colour})

    # define gradients for divisions
    for hue, divisions in division_colours.items():
        start, end = colour_wheel[hue]
        if len(divisions) == 1:
            gradient = linear_gradient(start, end, 3)
            gradient = [list(gradient)[1]]
        else:
            gradient = linear_gradient(start, end, len(divisions))
        for division, colour in zip(divisions, gradient):
            print('division', division, colour)
            results['division'].update({division: colour})

    # define gradients for locations
    for hue, locations in location_colours.items():
        start, end = colour_wheel[hue]
        if len(locations) == 1:
            gradient = linear_gradient(start, end, 3)
            gradient = [list(gradient)[1]]
        else:
            gradient = linear_gradient(start, end, len(locations))
        for location, colour in zip(locations, gradient):
            print('location', location, colour)
            results['location'].update({location: colour})


    # special colouring
    geoLevels = {}
    for line in scheme_list:
        if not line.startswith('\n'):
            line = line.strip()
            id = line.split('\t')[2]
            type = line.split('\t')[0]

            # parse subnational regions for countries in geoscheme
            if type == 'country':
                members = [item.strip() for item in line.split('\t')[5].split(',')] # elements inside the subarea
                if id not in geoLevels:
                    geoLevels[id] = members

    categories = {'Global': '#CCCCCC', 'Europe': '#666666'}
    results['us_region'] = {}
    usregion_hues = {
        'USA-Northeast': colour_scale['purple'][3],
        'USA-Midwest': colour_scale['green'][0],
        'USA-Southwest': colour_scale['yellow'][1],
        'USA-Southeast': colour_scale['cyan'][0],
        'USA-West': colour_scale['red'][0]
        }

    for us_region, hue in usregion_hues.items():
        start, end = hue_to_hex[hue]
        divisions = geoLevels[us_region]
        gradient = linear_gradient(start, end, len(divisions))
        # print(us_region, hue, divisions, gradient)
        for state, colour in zip(divisions, gradient):
            categories[state] = colour

    for reg, hex in categories.items():
        results['us_region'].update({reg: hex})
        print('us_region', reg, hex)

    # VOC / VOI
    list_category = [up_number for up_number in sorted(set(df['category'].to_list())) if up_number != 'Other variants']
    list_hex = list([hue_to_rgb(int(x)) for x in np.linspace(30, 240, len(list_category)*2, endpoint=True)])
    skip_hex = [h for n, h in enumerate(list_hex) if n in range(0, len(list_hex), 2)][::-1]

    results['category'] = {}
    for category, hex in zip(list_category, skip_hex):
        results['category'].update({category: hex})
        print(category, hex)
    results['category'].update({'Other variants': '#808080'})



    ''' EXPORT COLOUR FILE '''

    with open(output, 'w') as outfile:
        for trait, entries in results.items():
            for place, hexcolour in entries.items():
                if place in force_colour and trait not in ['location']:
                    hexcolour = force_colour[place]
                    print('* ' + place + ' is hardcode with the colour ' + hexcolour)
                line = "{}\t{}\t{}\n".format(trait, place, hexcolour.upper())
                outfile.write(line)
            outfile.write('\n')
    print('\nOrdered colour file successfully created!\n')

import pycountry_convert as pyCountry
import pycountry
import pandas as pd
import argparse
from uszipcode import SearchEngine


if __name__ == '__main__':
    parser = argparse.ArgumentParser(
        description="Reformat metadata file by adding column with subcontinental regions based on the UN geo-scheme",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter
    )
    parser.add_argument("--metadata", required=True, help="Nextstrain metadata file")
    parser.add_argument("--geoscheme", required=True, help="XML file with geographic classifications")
    parser.add_argument("--output", required=True, help="Updated metadata file")
    parser.add_argument("--filter", required=True, nargs='+', type=str, help="Filter region to define focus")
    args = parser.parse_args()

    metadata = args.metadata
    geoscheme = args.geoscheme
    output = args.output
    filt = args.filter

    # path = '/Users/anderson/GLab Dropbox/Anderson Brito/projects/ncov/ncov_variants/nextstrain/run6_20210202_b117/ncov/'
    # metadata = path + 'pre-analyses/metadata_filtered.tsv'
    # geoscheme = path + 'config/geoscheme.tsv'
    # output = path + 'pre-analyses/metadata_geo.tsv'
    #

    ##keep location (county) data for focus region of build
    if filt==['caribbean']:
            focus = ['US Virgin Islands','Puerto Rico','Dominican Republic','Anguilla','Antigua and Barbuda',
                'Aruba','Bahamas','Barbados','Bonaire','Sint Eustatius','Saba','British Virgin Islands','Cayman Islands',
                'Cuba','Curacao','Dominica','Grenada','Guadeloupe','Haiti','Jamaica','Martinique','Montserrat',
                'Saint Barthelemy','Saint Kitts and Nevis','Saint Lucia','Saint Martin','Saint Vincent and the Grenadines',
                'Sint Maarten','Trinidad and Tobago','Turks and Caicos', 'USA','Netherlands','United Kingdom','France']
    if filt == ['connecticut']:
            focus = ['USA', 'Canada', 'United Kingdom', 'Maine', 'New Hampshire',
             'Massachusetts', 'Connecticut', 'Vermont', 'New York']
    if filt == ['']:
        focus = []
    print(focus)

    # get ISO alpha3 country codes
    isos = {'British Virgin Islands':'VGB','Cayman Islands':'CYM','Guam':'GUM','US Virgin Islands':'VIR','Virgin Islands':'VIR',
            'Puerto Rico':'PRI','Northern Mariana Islands':'MNP','Sint Maarten':'MAF','Sint Eustatius':'BES','St Eustatius':'BES',
            'Bonaire':'BES','Saba':'BES','Martinique':'MTQ','French Guiana':'GUF','Montserrat':'MSR','Turks and Caicos':'TCA',
            'Anguilla':'AIA','Mayotte':'MYT','Reunion':'RUE','Wallis and Futuna':'WLF'}
    def get_iso(country):
        global isos
        if country not in isos.keys():
            try:
                isoCode = pyCountry.country_name_to_country_alpha3(country, cn_name_format="default")
                isos[country] = isoCode
            except:
                try:
                    isoCode = pycountry.countries.search_fuzzy(country)[0].alpha_3
                    isos[country] = isoCode
                except:
                    isos[country] = ''
        return isos[country]

    #keys to correct ISOs for colonies
    colonies =  {'British Virgin Islands':'VGB','Cayman Islands':'CYM','Guam':'GUM','US Virgin Islands':'VIR','Virgin Islands':'VIR',
                'Puerto Rico':'PRI','Northern Mariana Islands':'MNP','Sint Maarten':'MAF','Sint Eustatius':'BES','St Eustatius':'BES',
                'Bonaire':'BES','Saba':'BES','Martinique':'MTQ','French Guiana':'GUF','Montserrat':'MSR','Turks and Caicos':'TCA',
                'Anguilla':'AIA','Mayotte':'MYT','Reunion':'RUE','Wallis and Futuna':'WLF'}
    #keys to correct regions for colonies
    colony_regions =  {'British Virgin Islands':'Caribbean','Cayman Islands':'Caribbean','Guam':'Oceania','US Virgin Islands':'Caribbean','Virgin Islands':'Caribbean',
                'Puerto Rico':'Caribbean','Northern Mariana Islands':'Oceania','Sint Maarten':'Caribbean','Sint Eustatius':'Caribbean','St Eustatius':'Caribbean',
                'Bonaire':'Caribbean','Saba':'Caribbean','Martinique':'Caribbean','French Guiana':'South America','Montserrat':'Caribbean','Turks and Caicos':'Caribbean',
                'Anguilla':'Caribbean','Mayotte':'Eastern Africa','Reunion':'Eastern Africa','Wallis and Futuna':'Oceania'}
    #keys to correct common misspellings (in GISAID or in our metadata)
    misspelled =    {'Virgin Islands':'US Virgin Islands','St Eustatius':'Sint Eustatius','Turks and Caicos Islands':'Turks and Caicos'}

    # parse subcontinental regions in geoscheme
    scheme_list = open(geoscheme, "r").readlines()[1:]
    geoLevels = {}
    c = 0
    for line in scheme_list:
        if not line.startswith('\n'):
            id = line.split('\t')[2]
            type = line.split('\t')[0]
            if type == 'region':
                members = line.split('\t')[5].split(',') # elements inside the subarea
                for country in members:
                    iso = get_iso(country.strip())
                    geoLevels[iso] = id

            # parse subnational regions for countries in geoscheme
            if type == 'country':
                members = line.split('\t')[5].split(',') # elements inside the subarea
                for state in members:
                    if state.strip() not in geoLevels.keys():
                        geoLevels[state.strip()] = id

            # parse subareas for states in geoscheme
            if type == 'location':
                members = line.split('\t')[5].split(',')  # elements inside the subarea
                for zipcode in members:
                    if zipcode.strip() not in geoLevels.keys():
                        geoLevels[zipcode.strip()] = id


    # open metadata file as dataframe
    dfN = pd.read_csv(metadata, encoding='utf-8', sep='\t')
    try:
        dfN.insert(4, 'region', '')
    except:
        pass
    dfN['region'] = dfN['iso'].map(geoLevels) # add 'column' region in metadata
    dfN['us_region'] = ''
    listA = dfN['iso'] + dfN['iso'].map(geoLevels)
    listA.to_csv('listA.csv', sep='\t', index=False)

    # convert sets of locations into sub-locations
    print('\nApplying geo-schemes...')
    dfN.fillna('', inplace=True)
    for idx, row in dfN.iterrows():

        # flatten divison names as country names, for countries that are not a focus of study
        country = dfN.loc[idx, 'country']
        division = dfN.loc[idx,'division']
        if country not in focus:
            if division not in focus:
                dfN.loc[idx, 'division'] = country

        # assign US region
        if country not in ['USA']:
            if 'Europe' in dfN.loc[idx, 'region']:
                dfN.loc[idx, 'us_region'] = 'Europe'
            else:
                dfN.loc[idx, 'us_region'] = 'Global'
        if country == 'USA' and dfN.loc[idx, 'us_region'] == '':
            dfN.loc[idx, 'us_region'] = dfN.loc[idx, 'division']

        # divide country into subnational regions
        division = dfN.loc[idx, 'division']
        if division not in ['', 'unknown']:
            if division in geoLevels.keys():
                dfN.loc[idx, 'country'] = geoLevels[dfN.loc[idx, 'division']]

        # correct ISO codes and regions for colonies (any insular terrritories)
        island = dfN.loc[idx,'country']
        if island in colonies.keys():
            dfN.loc[idx,'iso'] = colonies[island]
            dfN.loc[idx,'region'] = colony_regions[island]


        # flatten location names as division names for divisions that are not a focus of study
        if division not in focus:
            dfN.loc[idx, 'location'] = division
        #print('Processing metadata for... ' + row['strain'])

        #rename some commonly misspelled or mismatched data
        if country in misspelled:
            dfN.loc[idx,'country'] = misspelled[country]
        if division in misspelled:
            dfN.loc[idx,'division'] = misspelled[division]
        location = dfN.loc[idx,'location']
        if location in misspelled:
            dfN.loc[idx,'location'] = misspelled[location]

        #rename Florida, PR
        if dfN.loc[idx,'country'] in ['Puerto Rico']:
            if dfN.loc[idx,'location'] in ['Florida']:
                dfN.loc[idx,'location'] = 'Florida (PR)'



    dfN = dfN.drop_duplicates(subset=['strain'])
    dfN.to_csv(output, sep='\t', index=False)

print('\nMetadata file successfully reformatted applying geo-scheme!\n')

import pycountry_convert as pyCountry
import pycountry
from Bio import SeqIO
import pandas as pd
from epiweeks import Week
import time
import argparse

if __name__ == '__main__':
    parser = argparse.ArgumentParser(
        description="Filter nextstrain metadata files re-formmating and exporting only selected lines",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter
    )
    parser.add_argument("--genomes", required=True, help="FASTA file genomes to be used")
    parser.add_argument("--metadata1", required=True, help="Metadata file from NextStrain")
    parser.add_argument("--metadata2", required=False, help="Custom lab metadata file")
    parser.add_argument("--filter", required=False, nargs='+', type=str,  help="List of filters for tagged rows in lab metadata")
    parser.add_argument("--output1", required=True, help="Filtered metadata file")
    parser.add_argument("--output2", required=True, help="Reformatted, final FASTA file")
    args = parser.parse_args()

    genomes = args.genomes
    metadata1 = args.metadata1
    metadata2 = args.metadata2
    filterby = args.filter
    output1 = args.output1
    output2 = args.output2

    # path = '/Users/anderson/GLab Dropbox/Anderson Brito/projects/ncov/ncov_variants/nextstrain/runX_20210617_filter/'
    # genomes = path + 'pre-analyses/temp_sequences.fasta'
    # metadata1 = path + 'pre-analyses/metadata_nextstrain.tsv'
    # metadata2 = path + 'pre-analyses/GLab_SC2_sequencing_data.xlsx'
    # filterby = ['caribe', 'test']
    # output1 = path + 'pre-analyses/metadata_filtered.tsv'
    # output2 = path + 'pre-analyses/sequences.fasta'

    # temporal boundaries
    today = time.strftime('%Y-%m-%d', time.gmtime())
    min_date = '2019-12-15'

    variants = {'B.1.1.7': 'Alpha (B.1.1.7)',
                'B.1.351': 'Beta (B.1.351)',
                'B.1.351.2': 'Beta (B.1.351.2)',
                'B.1.351.3': 'Beta (B.1.351.3)',
                'P.1': 'Gamma (P.1)',
                'P.1.1': 'Gamma (P.1.1)',
                'P.1.2': 'Gamma (P.1.2)',
                'B.1.617.2': 'Delta (B.1.617.2)',
                'AY.1': 'Delta (AY.1)',
                'AY.2': 'Delta (AY.2)',
                'AY.3': 'Delta (AY.3)',
                'AY.3.1': 'Delta (AY.3.1)',
                'AY.4': 'Delta (AY.4)',
                'AY.4.1': 'Delta (AY.4.1)',
                'AY.4.2': 'Delta (AY.4.2)',              
                'AY.4.3': 'Delta (AY.4.3)',
                'AY.4.4': 'Delta (AY.4.4)',
                'AY.4.5': 'Delta (AY.4.5)',
                'AY.5': 'Delta (AY.5)',
                'AY.5.1': 'Delta (AY.5.1)',
                'AY.5.2': 'Delta (AY.5.2)',
                'AY.5.3': 'Delta (AY.5.3)',
                'AY.5.4': 'Delta (AY.5.4)',
                'AY.6': 'Delta (AY.6)',
                'AY.7': 'Delta (AY.7)',
                'AY.7.1': 'Delta (AY.7.1)',
                'AY.7.2': 'Delta (AY.7.2)',
                'AY.8': 'Delta (AY.8)',
                'AY.9': 'Delta (AY.9)',
                'AY.9.1': 'Delta (AY.9.1)',
                'AY.9.2': 'Delta (AY.9.2)',
                'AY.9.2.1': 'Delta (AY.9.2.1)',               
                'AY.10': 'Delta (AY.10)',
                'AY.11': 'Delta (AY.11)',
                'AY.12': 'Delta (AY.12)',
                'AY.13': 'Delta (AY.13)',
                'AY.14': 'Delta (AY.14)',
                'AY.15': 'Delta (AY.15)',
                'AY.16': 'Delta (AY.16)',
                'AY.16.1': 'Delta (AY.16.1)',
                'AY.17': 'Delta (AY.17)',
                'AY.18': 'Delta (AY.18)',
                'AY.19': 'Delta (AY.19)',
                'AY.20': 'Delta (AY.20)',
                'AY.21': 'Delta (AY.21)',
                'AY.22': 'Delta (AY.22)',
                'AY.23': 'Delta (AY.23)',
                'AY.23.1': 'Delta (AY.23.1)',
                'AY.24': 'Delta (AY.24)',
                'AY.25': 'Delta (AY.25)',
                'AY.26': 'Delta (AY.26)',
                'AY.27': 'Delta (AY.27)',
                'AY.28': 'Delta (AY.28)',
                'AY.29': 'Delta (AY.29)',
                'AY.29.1': 'Delta (AY.29.1)',
                'AY.30': 'Delta (AY.30)',
                'AY.31': 'Delta (AY.31)',
                'AY.32': 'Delta (AY.32)',
                'AY.33': 'Delta (AY.33)',
                'AY.34': 'Delta (AY.34)',
                'AY.35': 'Delta (AY.35)',
                'AY.36': 'Delta (AY.36)',
                'AY.37': 'Delta (AY.37)',
                'AY.38': 'Delta (AY.38)',
                'AY.39': 'Delta (AY.39)',
                'AY.39.1': 'Delta (AY.39.1)',
                'AY.39.1.1': 'Delta (AY.39.1.1)',
                'AY.39.2': 'Delta (AY.39.2)',
                'AY.40': 'Delta (AY.40)',
                'AY.41': 'Delta (AY.41)',
                'AY.42': 'Delta (AY.42)',
                'AY.43': 'Delta (AY.43)',
                'AY.44': 'Delta (AY.44)',
                'AY.45': 'Delta (AY.45)',
                'AY.46': 'Delta (AY.46)',
                'AY.46.1': 'Delta (AY.46.1)',
                'AY.46.2': 'Delta (AY.46.2)',
                'AY.46.3': 'Delta (AY.46.3)',
                'AY.46.4': 'Delta (AY.46.4)',
                'AY.46.5': 'Delta (AY.46.5)',
                'AY.46.6': 'Delta (AY.46.6)',
                'AY.47': 'Delta (AY.47)',
                'AY.48': 'Delta (AY.48)',
                'AY.49': 'Delta (AY.49)',
                'AY.50': 'Delta (AY.50)',
                'AY.51': 'Delta (AY.51)',
                'AY.52': 'Delta (AY.52)',
                'AY.53': 'Delta (AY.53)',
                'AY.54': 'Delta (AY.54)',
                'AY.55': 'Delta (AY.55)',
                'AY.56': 'Delta (AY.56)',
                'AY.57': 'Delta (AY.57)',
                'AY.58': 'Delta (AY.58)',
                'AY.59': 'Delta (AY.59)',
                'AY.60': 'Delta (AY.60)',
                'AY.61': 'Delta (AY.61)',
                'AY.62': 'Delta (AY.62)',
                'AY.63': 'Delta (AY.63)',
                'AY.64': 'Delta (AY.64)',
                'AY.65': 'Delta (AY.65)',
                'AY.66': 'Delta (AY.66)',
                'AY.67': 'Delta (AY.67)',
                'AY.68': 'Delta (AY.68)',
                'AY.69': 'Delta (AY.69)',
                'AY.70': 'Delta (AY.70)',
                'AY.71': 'Delta (AY.71)',
                'AY.72': 'Delta (AY.72)',
                'AY.73': 'Delta (AY.73)',
                'AY.74': 'Delta (AY.74)',
                'AY.75': 'Delta (AY.75)',
                'AY.75.1': 'Delta (AY.75.1)',
                'AY.76': 'Delta (AY.76)',
                'AY.77': 'Delta (AY.77)',
                'AY.78': 'Delta (AY.78)',
                'AY.79': 'Delta (AY.79)',
                'AY.80': 'Delta (AY.80)',
                'AY.81': 'Delta (AY.81)',
                'AY.82': 'Delta (AY.82)',
                'AY.83': 'Delta (AY.83)',
                'AY.84': 'Delta (AY.84)',
                'AY.85': 'Delta (AY.85)',
                'AY.86': 'Delta (AY.86)',
                'AY.87': 'Delta (AY.87)',
                'AY.88': 'Delta (AY.88)',
                'AY.89': 'Delta (AY.89)',
                'AY.90': 'Delta (AY.90)',
                'AY.91': 'Delta (AY.91)',
                'AY.91.1': 'Delta (AY.91.1)',
                'AY.92': 'Delta (AY.92)',
                'AY.93': 'Delta (AY.93)',
                'AY.94': 'Delta (AY.94)',
                'AY.95': 'Delta (AY.95)',
                'AY.96': 'Delta (AY.96)',
                'AY.97': 'Delta (AY.97)',
                'AY.98': 'Delta (AY.98)',
                'AY.98.1': 'Delta (AY.98.1)',
                'AY.99': 'Delta (AY.99)',
                'AY.99.1': 'Delta (AY.99.1)',
                'AY.99.2': 'Delta (AY.99.2)',
                'AY.100': 'Delta (AY.100)',
                'AY.101': 'Delta (AY.101)',
                'AY.102': 'Delta (AY.102)',
                'AY.103': 'Delta (AY.103)',
                'AY.104': 'Delta (AY.104)',
                'AY.105': 'Delta (AY.105)',
                'AY.106': 'Delta (AY.106)',
                'AY.107': 'Delta (AY.107)',
                'AY.108': 'Delta (AY.108)',
                'AY.109': 'Delta (AY.109)',
                'AY.110': 'Delta (AY.110)',
                'AY.111': 'Delta (AY.111)',
                'AY.112': 'Delta (AY.112)',
                'AY.113': 'Delta (AY.113)',
                'AY.114': 'Delta (AY.114)',
                'AY.115': 'Delta (AY.115)',
                'AY.116': 'Delta (AY.116)',
                'AY.116.1': 'Delta (AY.116.1)',
                'AY.117': 'Delta (AY.117)', 
                'AY.118': 'Delta (AY.118)',    
                'AY.119': 'Delta (AY.119)',   
                'AY.120': 'Delta (AY.120)',   
                'AY.120.1': 'Delta (AY.120.1)', 
                'AY.120.2': 'Delta (AY.120.2)',
                'AY.120.2.1': 'Delta (AY.120.2.1)',
                'AY.121': 'Delta (AY.121)',
                'AY.121.1': 'Delta (AY.121.1)',
                'AY.122': 'Delta (AY.122)', 
                'AY.122.1': 'Delta (AY.122.1)',
                'AY.123': 'Delta (AY.123)',
                'AY.124': 'Delta (AY.124)', 
                'AY.125': 'Delta (AY.125)',
                'AY.126': 'Delta (AY.126)',
                'AY.127': 'Delta (AY.127)',
                'AY.128': 'Delta (AY.128)',
                'AY.129': 'Delta (AY.129',
                'AY.130': 'Delta (AY.130)',                                                                                                                                          
                'B.1.525': 'Eta (B.1.525)',
                'B.1.526': 'Iota (B.1.526)',
                'B.1.617.1': 'Kappa (B.1.617.1)',
                'C.37': 'Lambda (C.37)',
                'B.1.427': 'Epsilon (B.1.427/B.1.429)',
                'B.1.429': 'Epsilon (B.1.427/B.1.429)',
                'P.2': 'Zeta (P.2)',
                'B.1.621': 'Mu (B.1.621)',
                'B.1.621.1': 'Mu (B.1.621.1)',
                'BA.1': 'Omicron (BA.1)'
                }


    # get ISO alpha3 country codes

    isos =  {'US Virgin Islands':'VIR','Virgin Islands':'VIR','British Virgin Islands':'VGB','Curacao':'CUW','Northern Mariana Islands':'MNP',
            'Sint Maarten':'MAF','St Eustatius':'BES'}
    def get_iso(country):
        global isos
        if country not in isos.keys():
            try:
                isoCode = pyCountry.country_name_to_country_alpha3(country, cn_name_format="default")
                isos[country] = isoCode
            except:
                try:
                    isoCode = pycountry.countries.search_fuzzy(country)[0].alpha_3
                    isos[country] = isoCode
                except:
                    isos[country] = ''
        return isos[country]


    # create epiweek column
    def get_epiweeks(date):
        date = pd.to_datetime(date)
        epiweek = str(Week.fromdate(date, system="cdc"))  # get epiweeks
        epiweek = epiweek[:4] + '_' + 'EW' + epiweek[-2:]
        return epiweek


    # add state code
    us_state_abbrev = {
        'Alabama': 'AL',
        'Alaska': 'AK',
        'American Samoa': 'AS',
        'Arizona': 'AZ',
        'Arkansas': 'AR',
        'California': 'CA',
        'Colorado': 'CO',
        'Connecticut': 'CT',
        'Delaware': 'DE',
        'District of Columbia': 'DC',
        'Washington DC': 'DC',
        'Florida': 'FL',
        'Georgia': 'GA',
        'Guam': 'GU',
        'Hawaii': 'HI',
        'Idaho': 'ID',
        'Illinois': 'IL',
        'Indiana': 'IN',
        'Iowa': 'IA',
        'Kansas': 'KS',
        'Kentucky': 'KY',
        'Louisiana': 'LA',
        'Maine': 'ME',
        'Maryland': 'MD',
        'Massachusetts': 'MA',
        'Michigan': 'MI',
        'Minnesota': 'MN',
        'Mississippi': 'MS',
        'Missouri': 'MO',
        'Montana': 'MT',
        'Nebraska': 'NE',
        'Nevada': 'NV',
        'New Hampshire': 'NH',
        'New Jersey': 'NJ',
        'New Mexico': 'NM',
        'New York': 'NY',
        'North Carolina': 'NC',
        'North Dakota': 'ND',
        'Northern Mariana Islands': 'MP',
        'Ohio': 'OH',
        'Oklahoma': 'OK',
        'Oregon': 'OR',
        'Pennsylvania': 'PA',
        'Puerto Rico': 'PR',
        'Rhode Island': 'RI',
        'South Carolina': 'SC',
        'South Dakota': 'SD',
        'Tennessee': 'TN',
        'Texas': 'TX',
        'Utah': 'UT',
        'Vermont': 'VT',
        'Virgin Islands': 'VI',
        'Virginia': 'VA',
        'Washington': 'WA',
        'West Virginia': 'WV',
        'Wisconsin': 'WI',
        'Wyoming': 'WY'
    }

    # nextstrain metadata
    dfN = pd.read_csv(metadata1, encoding='utf-8', sep='\t', dtype='str')
    dfN['strain'] = dfN['strain'].replace('hCoV-19/', '')
    dfN.insert(4, 'iso', '')
    dfN.insert(1, 'category', '')
    dfN.fillna('', inplace=True)

    # add tag of variant category
    def variant_category(lineage):
        var_category = 'Other variants'
        for name in variants.keys():
            if lineage == name:
                var_category = variants[lineage]
        return var_category

    dfN['category'] = dfN['pango_lineage'].apply(lambda x: variant_category(x))


    list_columns = dfN.columns.values  # list of column in the original metadata file
    # print(dfN)

    # Lab genomes metadata
    dfE = pd.read_excel(metadata2, index_col=None, header=0, sheet_name=0,
                        converters={'Sample-ID': str, 'Collection-date': str})
    dfE.fillna('', inplace=True)

    dfE = dfE.rename(
        columns={'Sample-ID': 'id', 'Collection-date': 'date', 'Country': 'country', 'Division (state)': 'division',
                 'Location (county)': 'location', 'Country of exposure': 'country_exposure',
                 'State of exposure': 'division_exposure', 'Lineage': 'pango_lineage', 'Source': 'originating_lab',
                 'Filter': 'filter'})
    dfE['epiweek'] = ''

    # exclude rows with no ID
    if 'id' in dfE.columns.to_list():
        dfE = dfE[~dfE['id'].isin([''])]

    lab_sequences = dfE['id'].tolist()
    # exclude unwanted lab metadata row
    if len(filterby) > 0:
        print('\nFiltering metadata by category: ' + ', '.join(filterby) + '\n')
    dfL = pd.DataFrame(columns=dfE.columns.to_list())
    for value in filterby:
        dfF = dfE[dfE['filter'].isin([value])]  # batch inclusion of specific rows
        dfL = pd.concat([dfL, dfF]) # add filtered rows to dataframe with lab metadata

    # list of relevant genomes sequenced
    keep_only = dfL['id'].tolist()
    excluded = [id for id in lab_sequences if id not in keep_only]

    # create a dict of existing sequences
    sequences = {}
    for fasta in SeqIO.parse(open(genomes), 'fasta'):  # as fasta:
        id, seq = fasta.description, fasta.seq
        if id not in sequences.keys() and id not in excluded:
            sequences[id] = str(seq)

    # add inexistent columns
    for col in list_columns:
        if col not in dfL.columns:
            dfL[col] = ''

    # output dataframe
    outputDF = pd.DataFrame(columns=list_columns)
    found = []
    lab_label = {}
    metadata_issues = {}
    # process metadata from excel sheet
    for idx, row in dfL.iterrows():
        id = dfL.loc[idx, 'id'].replace('hCoV-19/', '')
        if id in sequences:
            dict_row = {}
            for col in list_columns:
                dict_row[col] = ''
                if col in row:
                    dict_row[col] = dfL.loc[idx, col].strip()  # add values to dictionary

            # check for missing geodata
            geodata = ['country'] # column
            for level in geodata:
                if len(dict_row[level]) < 1:
                    if id not in metadata_issues:
                        metadata_issues[id] = [level]
                    else:
                        metadata_issues[id].append(level)

            if dict_row['location'] in ['', None]:
                dict_row['location'] = dfL.loc[idx, 'location']

            collection_date = ''
            if len(str(dict_row['date'])) > 1 and 'X' not in dict_row['date']:
                collection_date = dict_row['date'].split(' ')[0].replace('.', '-').replace('/', '-')
                dict_row['date'] = collection_date
                # check is date is appropriate: not from the 'future', not older than 'min_date'
                if pd.to_datetime(today) < pd.to_datetime(collection_date) or pd.to_datetime(min_date) > pd.to_datetime(collection_date):
                    if id not in metadata_issues:
                        metadata_issues[id] = ['date']
                    else:
                        metadata_issues[id].append('date')
            else: # missing date
                if id not in metadata_issues:
                    metadata_issues[id] = ['date']
                else:
                    metadata_issues[id].append('date')

            # fix exposure
            columns_exposure = ['country_exposure', 'division_exposure']
            for level_exposure in columns_exposure:
                level = level_exposure.split('_')[0]
                dict_row[level_exposure] = dfL.loc[idx, level_exposure]
                if dict_row[level_exposure] in ['', None]:
                    if level_exposure == 'country_exposure':
                        dict_row[level_exposure] = dict_row[level]
                    else:
                        if dict_row['country_exposure'] != dfL.loc[idx, 'country']:
                            dict_row[level_exposure] = dict_row['country_exposure']
                        else:
                            dict_row[level_exposure] = dict_row[level]

            code = ''
            if dict_row['division'] in us_state_abbrev:
                code = us_state_abbrev[dict_row['division']] + '-'

            strain = dfL.loc[idx, 'country'].replace(' ', '') + '/' + code + dfL.loc[idx, 'id'] + '/' + collection_date.split('-')[0]

            # set the strain name
            dict_row['strain'] = strain
            dict_row['iso'] = get_iso(dict_row['country'])
            dict_row['originating_lab'] = dfL.loc[idx, 'originating_lab']
            dict_row['submitting_lab'] = 'Grubaugh Lab - Yale School of Public Health'
            dict_row['authors'] = 'GLab team'

            # add lineage
            lineage = ''
            if dfL.loc[idx, 'pango_lineage'] != '':
                lineage = dfL.loc[idx, 'pango_lineage']
            dict_row['pango_lineage'] = lineage

            # variant classication (VOI, VOC, VHC)
            dict_row['category'] = variant_category(lineage)

            # assign epiweek
            if len(dict_row['date']) > 0 and 'X' not in dict_row['date']:
                dict_row['epiweek'] = get_epiweeks(collection_date)
            else:
                dict_row['epiweek'] = ''

            # record sequence and metadata as found
            found.append(strain)
            if id not in metadata_issues.keys():
                lab_label[id] = strain
                outputDF = outputDF.append(dict_row, ignore_index=True)


    # process metadata from TSV
    dfN = dfN[dfN['strain'].isin(sequences.keys())]
    for idx, row in dfN.iterrows():
        strain = dfN.loc[idx, 'strain'].replace('hCoV-19/', '')
        if strain in sequences:
            if strain in outputDF['strain'].to_list():
                continue
            dict_row = {}
            date = ''
            for col in list_columns:
                if col == 'date':
                    date = dfN.loc[idx, col]
                dict_row[col] = ''
                if col in row:
                    dict_row[col] = dfN.loc[idx, col]

            # fix exposure
            columns_exposure = ['country_exposure', 'division_exposure']
            for level_exposure in columns_exposure:
                level = level_exposure.split('_')[0]
                if dict_row[level_exposure] in ['', None]:
                    dict_row[level_exposure] = dict_row[level]

            dict_row['iso'] = get_iso(dict_row['country'])
            dict_row['epiweek'] = get_epiweeks(date)
            found.append(strain)

            outputDF = outputDF.append(dict_row, ignore_index=True)

    # write new metadata files
    outputDF = outputDF.drop(columns=['region'])
    outputDF.to_csv(output1, sep='\t', index=False)

    # write sequence file
    exported = []
    with open(output2, 'w') as outfile2:
        # export new metadata lines
        for id, sequence in sequences.items():
            if id in lab_label and id not in metadata_issues.keys(): # export lab generated sequences
                if lab_label[id] not in exported:
                    entry = '>' + lab_label[id] + '\n' + sequence + '\n'
                    outfile2.write(entry)
                    print('* Exporting newly sequenced genome and metadata for ' + id)
                    exported.append(lab_label[id])
            else:  # export publicly available sequences
                if id not in exported and id in outputDF['strain'].tolist():
                    entry = '>' + id + '\n' + sequence + '\n'
                    outfile2.write(entry)
                    exported.append(id)

    if len(metadata_issues) > 0:
        print('\n\n### WARNINGS!\n')
        print('\nPlease check for metadata issues related to these samples and column (which will be otherwise ignored)\n')
        for id, columns in metadata_issues.items():
            print('\t- ' + id + ' (issues found at: ' + ', '.join(columns) + ')')

    print('\nMetadata file successfully reformatted and exported!\n')