compara-deep-learning/pfam_matrix.py

import pandas as pd
import numpy as np
import progressbar
import os
import json
import sys
from prepare_synteny_matrix import read_data_homology, load_neighbor_genes
from process_data import create_map_list
from process_negative import read_database_txt


def get_score_overlap(x, y, pfam_db, pfam_map):
    df_1 = pfam_db.loc[pfam_map[x]]
    df_2 = pfam_db.loc[pfam_map[y]]
    list_of_domains = list(df_2.domain)
    c = 0
    c_1 = 0
    for _, row in df_1.iterrows():
        # check if the domain exists in the list
        if row.domain in list_of_domains:
            c_1 += 1
            # get the start
            st = int(df_2[df_2["domain"] == row.domain].hmm_from)
            # get the end
            end = int(df_2[df_2["domain"] == row.domain].hmm_to)
            # check if the domain is a ovelapping domain
            if (int(
                row.hmm_from) > st and int(row.hmm_from) < end) or (int(
                    row.hmm_to) > st and int(row.hmm_from) < end):
                c += 1
    return c/max(len(df_1), len(df_2))


def pfam_matrix(g1, g2, n, pfam_db, gmap, pfam_map):
    pm = np.zeros((n, n))
    for i in range(n):
        if g1[i] == "NULL_GENE":
            continue
        try:
            _ = gmap[g1[i]]
        except Exception:
            continue
        for j in range(n):
            if g2[j] == "NULL_GENE":
                continue
            try:
                _ = gmap[g2[j]]
            except Exception:
                continue
            pm[i][j] = get_score_overlap(g1[i], g2[j], pfam_db, pfam_map)
    return pm


def create_pfam_matrix(df, lsy, pfam_db, pfam_map):
    n = 3
    glist = list(pfam_db.gene_stable_id)
    gmap = create_map_list(glist)
    pg = []
    indexes = []
    for index, row in progressbar.progressbar(df.iterrows()):
        g1 = str(row["gene_stable_id"])
        g2 = str(row["homology_gene_stable_id"])
        x = []
        y = []
        try:
            _ = lsy[g1]
            _ = lsy[g2]
        except Exception:
            continue
        for i in range(len(lsy[g1]['b'])-1, -1, -1):
            x.append(lsy[g1]['b'][i])
        x.append(g1)
        for k in lsy[g1]['f']:
            x.append(k)

        for i in range(len(lsy[g2]['b'])-1, -1, -1):
            y.append(lsy[g2]['b'][i])
        y.append(g2)
        for k in lsy[g2]['f']:
            y.append(k)

        assert(len(x) == len(y))
        assert(len(x) == (2*n+1))
        pmtemp = pfam_matrix(x, y, 2*n+1, pfam_db, gmap, pfam_map)
        pg.append(pmtemp)
        indexes.append(index)
    return np.array(pg), np.array(indexes)


def create_pfam_map(pfam_db):
    pfam_map = {}
    for index, row in progressbar.progressbar(pfam_db.iterrows()):
        try:
            _ = pfam_map[row.gene_stable_id]
        except Exception:
            pfam_map[row.gene_stable_id] = []
        pfam_map[row.gene_stable_id].append(index)
    return pfam_map


def main_positive():
    if not os.path.isdir("processed/pfam_matrices"):
        os.mkdir("processed/pfam_matrices")
    a_h, d_h = read_data_homology("data_homology")
    lsy = load_neighbor_genes()
    pfam_db = pd.read_hdf("pfam_db_positive.h5")
    pfam_map = create_pfam_map(pfam_db)
    ndir = "processed/pfam_matrices/"
    nf1 = "pfam_matrices"
    nf3 = "pfam_indexes"
    for i in range(len(a_h)):
        df = a_h[i]
        print(len(df))
        pfam_matrices, indexes = create_pfam_matrix(df, lsy, pfam_db, pfam_map)
        np.save(ndir+str(d_h[i])+"_"+nf1, pfam_matrices)
        np.save(ndir+str(d_h[i])+"_"+nf3, indexes)
        print(len(indexes))


def read_data_negative(arg):
    df = read_database_txt(arg[-1])
    name = arg[-1].split(".")[0]
    ind = np.load("processed/synteny_matrices/"+name+"_indexes.npy")
    df = df.loc[ind]
    pfam_db = pd.read_hdf("pfam_db_negative.h5")
    pfam_map = create_pfam_map(pfam_db)
    with open("processed/neighbor_genes_negative.json", "r") as file:
        lsy = dict(json.load(file))
    return df, pfam_db, pfam_map, lsy, name


def main_negative(arg):
    df, pfam_db, pfam_map, lsy, name = read_data_negative(arg)
    ndir = "processed/pfam_matrices/"
    nf1 = "pfam_matrices"
    nf3 = "pfam_indexes"
    print(len(df))
    pfam_matrices, indexes = create_pfam_matrix(df, lsy, pfam_db, pfam_map)
    np.save(ndir+name+"_"+nf1, pfam_matrices)
    np.save(ndir+name+"_"+nf3, indexes)
    print(len(indexes))


def main():
    arg = sys.argv
    main_positive()
    main_negative(arg)


if __name__ == "__main__":
    main()