<reponame>cherrytrees-kpu/sturdy-journey
import argparse
import pathlib
import json
import pandas
# data.loc[data['acc'].isin(accs)]
def import_assays(assay_path): 
    """
    Import assays into list
    Parameters: 
    assay_path - pathlib.Path - path to assay.xlsx file
    Output: 
    list: 
    [0] - assay1
    [1] - assay2
    ...
    [x] - assayx
    """
    data = pandas.read_excel(assay_path, engine='openpyxl')
    return data
def import_archive(archive_path):
    """
    Import archive data
    Parameters:
    archive_path - path to archive files
    Output: 
    fw_dict, probes_dict, rev_dict, key = sequence of oligo
    """
    def import_data(file_path):
        input_file = open(file_path, 'r')
        data = json.loads(input_file.readline())
        oligo_data = dict()
        for item in data: 
            oligo_data[item['seq']] = item
        return oligo_data
    for oligo_path in archive_path.glob('*.txt'):
        file_oligo = oligo_path.stem.split('_')[1]
        if file_oligo == 'fw':
            fw_dict = import_data(oligo_path)
        elif file_oligo == 'rev': 
            rev_dict = import_data(oligo_path)
        elif file_oligo =='probes':
            probes_dict = import_data(oligo_path)
    return(fw_dict, probes_dict, rev_dict)
def import_database(database_path):
    data = pandas.read_csv(database_path)
    return data
def import_genus_blastdb_data(genus_regions_path, blastdb_regions_path, genus_blastdb_map_path):
    genus_regions_file = open(genus_regions_path, 'r')
    blastdb_regions_file = open(blastdb_regions_path, 'r')
    genus_blastdb_map_file = open(genus_blastdb_map_path, 'r')
    genus_regions = json.loads(genus_regions_file.readline())
    blastdb_regions = json.loads(blastdb_regions_file.readline())
    genus_blastdb_map = json.loads(genus_blastdb_map_file.readline())
    genus_regions_file.close()
    blastdb_regions_file.close()
    genus_blastdb_map_file.close()
    return (genus_regions, blastdb_regions, genus_blastdb_map)
def calc_assay_sens_spec(fw, probe, rev, genus_regions, blastdb_regions, genus_blastdb_map): 
    def check_within_regions(sequence_regions, start, end): 
        """Identify which sequences contain the target region.
        Note that the target region is considered continuous from 
        the specified start and end coordinates.
        Parameters: 
        sequence_regions - dictionary returned from get_sequence_regions
        [accession_id]:(seq_start, seq_end)
        start - start of target regoin
        end - end of target region
        Output: 
        traversed_regions - list of accession IDs that contain
        the target region
        """
        traversed_regions = []
        curr_index = 0
        for accession in sequence_regions: 
            region_start = sequence_regions[accession][0]
            region_end = sequence_regions[accession][1]
            if (
                region_start <= start
                and region_end >= end
            ): 
                traversed_regions.append(accession)
            curr_index = curr_index + 1
            if curr_index%1000 == 0: 
                print(f'Current progress: {curr_index}')
        return traversed_regions
    def calc_assay_sens(fw, probe, rev, genus_regions):
        #Calculate total possible targets
        assay_start = fw['start']
        assay_end = rev['end']
        possible_assay_targets = check_within_regions(genus_regions, assay_start, assay_end)
        #Identify accessions that are in all three
        assay_amplified_accs = []
        for acc in probe['targets']:
            if acc in fw['targets'] and acc in rev['targets']:
                assay_amplified_accs.append(acc)
        sens = len(assay_amplified_accs)/len(possible_assay_targets)
        return (sens, assay_amplified_accs)
    def calc_assay_spec(fw, probe, rev, genus_regions, blastdb_regions, genus_blastdb_map):
        #Calculate total possible non-targets
        assay_start = fw['start']
        assay_end = rev['end']
        mapped_assay_start = genus_blastdb_map[assay_start]
        mapped_assay_end = genus_blastdb_map[assay_end]
        possible_assay_targets = check_within_regions(blastdb_regions, mapped_assay_start, mapped_assay_end)
        #Remove accs that belong to the genus from possible_assay_targets
        for acc in genus_regions.keys():
            if acc in possible_assay_targets: 
                possible_assay_targets.remove(acc)
        #Identify accessions that are in all three
        assay_amplified_accs = []
        for acc in probe['non_targets']:
            if acc in fw['non_targets'] and acc in rev['non_targets']:
                assay_amplified_accs.append(acc)
        spec = (len(possible_assay_targets) - len (assay_amplified_accs))/len(possible_assay_targets)
        return (spec, assay_amplified_accs)
    sens, targets = calc_assay_sens(fw, probe, rev, genus_regions)
    spec, non_targets = calc_assay_spec(fw, probe, rev, genus_regions, blastdb_regions, genus_blastdb_map)
    return (sens, targets, spec, non_targets)
def get_metadata(accessions, metadata_database):
    #Quirk - metadata accession has .1 as a part of the accession id
    accessions_fixed = []
    for accession in accessions: 
        accessions_fixed.append(accession+'.1')
    accessions_data = metadata_database.loc[metadata_database['acc'].isin(accessions_fixed)]
    return accessions_data
    
def parse_args():
    parser = argparse.ArgumentParser(description='MetaData analyzer - output information for an assay')
    parser.add_argument('assay_path', 
        metavar='path_to_assay', 
        action='store', 
        type=str, 
        help = 'Path to assay.tab file'
    )
    parser.add_argument('archive_path', 
        metavar='path_to_archive', 
        action='store', 
        type=str, 
        help = 'Path to folder containing all three archive files (fw, pb, rev). Ensure that the files are suffixed _fw.txt, _probes.txt, _rev.txt '
    )
    parser.add_argument('database_path',
        metavar='path_to_database_metadata', 
        action='store', 
        type=str, 
        help = 'Path to database metadata file'
    )
    parser.add_argument('blastdb_regions',
        metavar='path_to_blastdb_regions', 
        action='store', 
        type=str, 
        help = 'Path to blast db region files'
    )
    parser.add_argument('genus_regions',
        metavar='path_to_genus_regions', 
        action='store', 
        type=str, 
        help = 'Path to genus region files'
    )
    parser.add_argument('genus_blastdb_map',
        metavar='path_to_genus_blastdb_map', 
        action='store', 
        type=str, 
        help = 'Path to the genus-blastdb map'
    )
    args = parser.parse_args()
    assay_path = pathlib.Path(args.assay_path)
    archive_path = pathlib.Path(args.archive_path)
    database_path = pathlib.Path(args.database_path)
    genus_regions_path = pathlib.Path(args.genus_regions)
    blastdb_regions_path = pathlib.Path(args.blastdb_regions)
    genus_blastdb_map_path = pathlib.Path(args.genus_blastdb_map)
    return(assay_path, archive_path, database_path, genus_regions_path, blastdb_regions_path, genus_blastdb_map_path)
def main(assay_path, archive_path, database_path, genus_regions_path, blastdb_regions_path, genus_blastdb_map_path):
    output_path = assay_path.parent
    assays = import_assays(assay_path)
    fw_dict, probe_dict, rev_dict = import_archive(archive_path)
    metadata_database = import_database(database_path)
    genus_regions, blastdb_regions, genus_blastdb_map = import_genus_blastdb_data(genus_regions_path, blastdb_regions_path, genus_blastdb_map_path)
    list_sens = []
    list_spec = []
    i = 0
    for assay in assays.iterrows(): 
        fw_seq = assay[1]['FW']
        rev_seq = assay[1]['REV']
        probe_seq = assay[1]['PB']
        fw_data = fw_dict[fw_seq]
        rev_data = rev_dict[rev_seq]
        probe_data = probe_dict[probe_seq]
        sens, targets, spec, non_targets = calc_assay_sens_spec(fw_data, probe_data, rev_data, genus_regions, blastdb_regions, genus_blastdb_map)
        list_sens.append(sens)
        list_spec.append(spec)
        target_metadata = get_metadata(targets, metadata_database)
        non_target_metadata = get_metadata(non_targets, metadata_database)
        target_metadata_file = open(output_path.joinpath(f'assay_{i}_target_metadata.csv'), 'w', newline='')
        non_target_metadata_file = open(output_path.joinpath(f'assay_{i}_non_target_metadata.csv'), 'w', newline='')
        target_metadata.to_csv(target_metadata_file)
        non_target_metadata.to_csv(non_target_metadata_file)
        i = i + 1
    
    assays = assays.assign(sens=pandas.Series(list_sens).values)
    assays = assays.assign(spec=pandas.Series(list_spec).values)
    output_assays = open(output_path.joinpath('assay_complete.csv'), 'w', newline='')
    assays.to_csv(output_assays)

if __name__ == '__main__': 
    assay_path, archive_path, database_path, genus_regions_path, blastdb_regions_path, genus_blastdb_map_path = parse_args()
    main(assay_path, archive_path, database_path, genus_regions_path, blastdb_regions_path, genus_blastdb_map_path)

import argparse
import pathlib
import json
import pandas
# data.loc[data['acc'].isin(accs)]
def import_assays(assay_path): 
    """
    Import assays into list
    Parameters: 
    assay_path - pathlib.Path - path to assay.xlsx file
    Output: 
    list: 
    [0] - assay1
    [1] - assay2
    ...
    [x] - assayx
    """
    data = pandas.read_excel(assay_path, engine='openpyxl')
    return data
def import_archive(archive_path):
    """
    Import archive data
    Parameters:
    archive_path - path to archive files
    Output: 
    fw_dict, probes_dict, rev_dict, key = sequence of oligo
    """
    def import_data(file_path):
        input_file = open(file_path, 'r')
        data = json.loads(input_file.readline())
        oligo_data = dict()
        for item in data: 
            oligo_data[item['seq']] = item
        return oligo_data
    for oligo_path in archive_path.glob('*.txt'):
        file_oligo = oligo_path.stem.split('_')[1]
        if file_oligo == 'fw':
            fw_dict = import_data(oligo_path)
        elif file_oligo == 'rev': 
            rev_dict = import_data(oligo_path)
        elif file_oligo =='probes':
            probes_dict = import_data(oligo_path)
    return(fw_dict, probes_dict, rev_dict)
def import_database(database_path):
    data = pandas.read_csv(database_path)
    return data
def import_genus_blastdb_data(genus_regions_path, blastdb_regions_path, genus_blastdb_map_path):
    genus_regions_file = open(genus_regions_path, 'r')
    blastdb_regions_file = open(blastdb_regions_path, 'r')
    genus_blastdb_map_file = open(genus_blastdb_map_path, 'r')
    genus_regions = json.loads(genus_regions_file.readline())
    blastdb_regions = json.loads(blastdb_regions_file.readline())
    genus_blastdb_map = json.loads(genus_blastdb_map_file.readline())
    genus_regions_file.close()
    blastdb_regions_file.close()
    genus_blastdb_map_file.close()
    return (genus_regions, blastdb_regions, genus_blastdb_map)
def calc_assay_sens_spec(fw, probe, rev, genus_regions, blastdb_regions, genus_blastdb_map): 
    def check_within_regions(sequence_regions, start, end): 
        """Identify which sequences contain the target region.
        Note that the target region is considered continuous from 
        the specified start and end coordinates.
        Parameters: 
        sequence_regions - dictionary returned from get_sequence_regions
        [accession_id]:(seq_start, seq_end)
        start - start of target regoin
        end - end of target region
        Output: 
        traversed_regions - list of accession IDs that contain
        the target region
        """
        traversed_regions = []
        curr_index = 0
        for accession in sequence_regions: 
            region_start = sequence_regions[accession][0]
            region_end = sequence_regions[accession][1]
            if (
                region_start <= start
                and region_end >= end
            ): 
                traversed_regions.append(accession)
            curr_index = curr_index + 1
            if curr_index%1000 == 0: 
                print(f'Current progress: {curr_index}')
        return traversed_regions
    def calc_assay_sens(fw, probe, rev, genus_regions):
        #Calculate total possible targets
        assay_start = fw['start']
        assay_end = rev['end']
        possible_assay_targets = check_within_regions(genus_regions, assay_start, assay_end)
        #Identify accessions that are in all three
        assay_amplified_accs = []
        for acc in probe['targets']:
            if acc in fw['targets'] and acc in rev['targets']:
                assay_amplified_accs.append(acc)
        sens = len(assay_amplified_accs)/len(possible_assay_targets)
        return (sens, assay_amplified_accs)
    def calc_assay_spec(fw, probe, rev, genus_regions, blastdb_regions, genus_blastdb_map):
        #Calculate total possible non-targets
        assay_start = fw['start']
        assay_end = rev['end']
        mapped_assay_start = genus_blastdb_map[assay_start]
        mapped_assay_end = genus_blastdb_map[assay_end]
        possible_assay_targets = check_within_regions(blastdb_regions, mapped_assay_start, mapped_assay_end)
        #Remove accs that belong to the genus from possible_assay_targets
        for acc in genus_regions.keys():
            if acc in possible_assay_targets: 
                possible_assay_targets.remove(acc)
        #Identify accessions that are in all three
        assay_amplified_accs = []
        for acc in probe['non_targets']:
            if acc in fw['non_targets'] and acc in rev['non_targets']:
                assay_amplified_accs.append(acc)
        spec = (len(possible_assay_targets) - len (assay_amplified_accs))/len(possible_assay_targets)
        return (spec, assay_amplified_accs)
    sens, targets = calc_assay_sens(fw, probe, rev, genus_regions)
    spec, non_targets = calc_assay_spec(fw, probe, rev, genus_regions, blastdb_regions, genus_blastdb_map)
    return (sens, targets, spec, non_targets)
def get_metadata(accessions, metadata_database):
    #Quirk - metadata accession has .1 as a part of the accession id
    accessions_fixed = []
    for accession in accessions: 
        accessions_fixed.append(accession+'.1')
    accessions_data = metadata_database.loc[metadata_database['acc'].isin(accessions_fixed)]
    return accessions_data
    
def parse_args():
    parser = argparse.ArgumentParser(description='MetaData analyzer - output information for an assay')
    parser.add_argument('assay_path', 
        metavar='path_to_assay', 
        action='store', 
        type=str, 
        help = 'Path to assay.tab file'
    )
    parser.add_argument('archive_path', 
        metavar='path_to_archive', 
        action='store', 
        type=str, 
        help = 'Path to folder containing all three archive files (fw, pb, rev). Ensure that the files are suffixed _fw.txt, _probes.txt, _rev.txt '
    )
    parser.add_argument('database_path',
        metavar='path_to_database_metadata', 
        action='store', 
        type=str, 
        help = 'Path to database metadata file'
    )
    parser.add_argument('blastdb_regions',
        metavar='path_to_blastdb_regions', 
        action='store', 
        type=str, 
        help = 'Path to blast db region files'
    )
    parser.add_argument('genus_regions',
        metavar='path_to_genus_regions', 
        action='store', 
        type=str, 
        help = 'Path to genus region files'
    )
    parser.add_argument('genus_blastdb_map',
        metavar='path_to_genus_blastdb_map', 
        action='store', 
        type=str, 
        help = 'Path to the genus-blastdb map'
    )
    args = parser.parse_args()
    assay_path = pathlib.Path(args.assay_path)
    archive_path = pathlib.Path(args.archive_path)
    database_path = pathlib.Path(args.database_path)
    genus_regions_path = pathlib.Path(args.genus_regions)
    blastdb_regions_path = pathlib.Path(args.blastdb_regions)
    genus_blastdb_map_path = pathlib.Path(args.genus_blastdb_map)
    return(assay_path, archive_path, database_path, genus_regions_path, blastdb_regions_path, genus_blastdb_map_path)
def main(assay_path, archive_path, database_path, genus_regions_path, blastdb_regions_path, genus_blastdb_map_path):
    output_path = assay_path.parent
    assays = import_assays(assay_path)
    fw_dict, probe_dict, rev_dict = import_archive(archive_path)
    metadata_database = import_database(database_path)
    genus_regions, blastdb_regions, genus_blastdb_map = import_genus_blastdb_data(genus_regions_path, blastdb_regions_path, genus_blastdb_map_path)
    list_sens = []
    list_spec = []
    i = 0
    for assay in assays.iterrows(): 
        fw_seq = assay[1]['FW']
        rev_seq = assay[1]['REV']
        probe_seq = assay[1]['PB']
        fw_data = fw_dict[fw_seq]
        rev_data = rev_dict[rev_seq]
        probe_data = probe_dict[probe_seq]
        sens, targets, spec, non_targets = calc_assay_sens_spec(fw_data, probe_data, rev_data, genus_regions, blastdb_regions, genus_blastdb_map)
        list_sens.append(sens)
        list_spec.append(spec)
        target_metadata = get_metadata(targets, metadata_database)
        non_target_metadata = get_metadata(non_targets, metadata_database)
        target_metadata_file = open(output_path.joinpath(f'assay_{i}_target_metadata.csv'), 'w', newline='')
        non_target_metadata_file = open(output_path.joinpath(f'assay_{i}_non_target_metadata.csv'), 'w', newline='')
        target_metadata.to_csv(target_metadata_file)
        non_target_metadata.to_csv(non_target_metadata_file)
        i = i + 1
    
    assays = assays.assign(sens=pandas.Series(list_sens).values)
    assays = assays.assign(spec=pandas.Series(list_spec).values)
    output_assays = open(output_path.joinpath('assay_complete.csv'), 'w', newline='')
    assays.to_csv(output_assays)

if __name__ == '__main__': 
    assay_path, archive_path, database_path, genus_regions_path, blastdb_regions_path, genus_blastdb_map_path = parse_args()
    main(assay_path, archive_path, database_path, genus_regions_path, blastdb_regions_path, genus_blastdb_map_path)
# data.loc[data['acc'].isin(accs)] Import assays into list     Parameters:      assay_path - pathlib.Path - path to assay.xlsx file     Output:      list:      [0] - assay1     [1] - assay2     ...     [x] - assayx Import archive data     Parameters:     archive_path - path to archive files     Output:      fw_dict, probes_dict, rev_dict, key = sequence of oligo Identify which sequences contain the target region.         Note that the target region is considered continuous from          the specified start and end coordinates.         Parameters:          sequence_regions - dictionary returned from get_sequence_regions         [accession_id]:(seq_start, seq_end)         start - start of target regoin         end - end of target region         Output:          traversed_regions - list of accession IDs that contain         the target region #Calculate total possible targets #Identify accessions that are in all three #Calculate total possible non-targets #Remove accs that belong to the genus from possible_assay_targets #Identify accessions that are in all three #Quirk - metadata accession has .1 as a part of the accession id