compara-deep-learning/process_data.py

import numpy as np
import progressbar
from save_data import write_dict_json


def create_map_list(l):  # this function maps the indexes to values
    t = {}
    for i in range(len(l)):
        t[l[i]] = i
    return t


def create_chromosome_maps(a, n):
    cmap = []
    cimap = []
    for df in progressbar.progressbar(a):
        chmap = {}
        chindmap = {}
        for index, row in df.iterrows():
            g = row.gene_id
            try:
                _ = chmap[g]
            except BaseException:
                chmap[g] = str(row.Chr)
                if str(row.Chr) in chindmap:
                    chindmap[str(row.Chr)].append(index)
                else:
                    chindmap[str(row.Chr)] = []
                    chindmap[str(row.Chr)].append(index)
        cmap.append(chmap)
        cimap.append(chindmap)
    return cmap, cimap


def list_dict_genomes(a, n):
    lst = []
    ldt = []
    for x in a:
        ldgt = {}
        uc = list(x["gene_id"])
        for i, r in x.iterrows():
            ldgt[r.gene_id] = i
        lst.append(uc)
        ldt.append(ldgt)
    return lst, ldt


def get_nearest_neighbors(g, gs, n, a, d, ld, ldg, cmap, cimap):
    # print("Finding Neighbor Genes")
    ne = []  # list to store the backward genes
    nr = []  # list to store the forward genes
    # get the address of the corresponding species to which the gene belongs
    # whose neighbor has to be found
    gi = d[gs.capitalize()]
    sldf = a[gi]  # select the dataframe
    scmap = cmap[gi]  # select the correct chromosome map
    scimap = cimap[gi]  # select the correct index maps
    try:
        _ = ld[gi]  # see if the corresponding gene map exists
    except BaseException:
        # print("Length of Dataframes:{} \t Length of Loaded Genes:{} \t Length of Loaded Genomes Dictionaries:{}".format(len(a),len(ld),len(ldg)))
        return ne, nr
    sldg = ldg[gi]  # select the corresponding map
    if g not in sldg:  # if the gene is not present in the dataframe return empty lists
        # print(g,"\t",gs)
        return ne, nr
    i = sldg[g]  # find the index of the gene
    chromosome_id = scmap[g]  # get the chromosome no from the database.
    scimap = scimap[chromosome_id]  # select the correct chromosomes indexes
    sldf = sldf.loc[scimap]  # select only the same chromosome genes.
    # get the -n neighbors
    start = int(sldf.loc[i]['start'])  # get the start location of the gene
    flag = 0
    for j in range(n):
        if flag == 1:
            ne.append("NULL_GENE")
            continue
        itemp = 0
        # select the column
        end = list(sldf.end)
        end = np.array(end)
        assert(len(end) == len(sldf))
        end = end - start  # subtract start from it so as to get relative position
        end_s = np.argsort(end)  # sort them by the order of distance
        if end[end_s[0]] >= 0:  # if all the genes end ahead of the one in consideration
            flag = 1  # increment the pointer
            ne.append("NULL_GENE")  # append the NULL_GENE value
            continue
        for k in end_s:  # iterate through the sorted array
            # find the first value that is negative and the next one is
            # positive to get the nearest gene
            if end[k] < 0 and end[k + 1] >= 0:
                itemp = k
                break
        itemp = scimap[itemp]
        ne.append(sldf.loc[itemp].gene_id)
        # make "start" the start location of the current gene
        start = int(sldf.loc[itemp].start)
        # print(start)
    # get the +n neighbors
    flag = 0
    end = int(sldf.loc[i].end)
    for j in range(n):
        if flag == 1:
            nr.append("NULL_GENE")
            continue
        itemp = 0
        start = list(sldf.start)
        start = np.array(start)
        start = start - end
        start_s = np.argsort(start)
        if start[start_s[-1]] < 0:
            flag = 1
            nr.append("NULL_GENE")
            continue
        for k in start_s:
            if start[k] > 0:
                itemp = k
                break
        itemp = scimap[itemp]
        nr.append(sldf.loc[itemp].gene_id)
        end = int(sldf.loc[itemp].end)
    return ne, nr


def create_data_homology_ls(a_h, d_h, n, a, d, ld, ldg, cmap, cimap, to_write):
    # dictionary which stores +/- n genes of the given gene by id. Each key is
    # a gene id which corresponds to the one in center.
    lsy = {}
    lsytemp = {}
    name = "neighbor_genes"
    for df in a_h:
        for _, row in progressbar.progressbar(df.iterrows()):
            x = row["gene_stable_id"]
            y = row["homology_gene_stable_id"]
            xs = row["species"]
            ys = row["homology_species"]
            try:
                _ = lsy[x]
            except BaseException:
                try:
                    # see if the species exist in genomic maps
                    _ = d[xs.capitalize()]
                    xl, xr = get_nearest_neighbors(
                        x, xs, n, a, d, ld, ldg, cmap, cimap)
                    if len(
                            xl) != 0:  # check if neighboring genes were successfully found
                        lsy[x] = dict(b=xl, f=xr)
                        lsytemp[x] = dict(b=xl, f=xr)
                except BaseException:
                    continue
            try:
                _ = lsy[y]
            except BaseException:
                try:
                    _ = d[ys.capitalize()]
                    yl, yr = get_nearest_neighbors(
                        y, ys, n, a, d, ld, ldg, cmap, cimap)
                    if len(yl) != 0:
                        lsy[y] = dict(b=yl, f=yr)
                        lsytemp[y] = dict(b=yl, f=yr)
                except BaseException:
                    continue
    if to_write == 1:
        write_dict_json(name, "processed", lsy)
    return lsy