Communityc2st.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from sklearn.neighbors import KNeighborsClassifier
import os
import keras
import numpy as np
import random
import pandas as pd
from sklearn.cluster import AgglomerativeClustering
from sklearn.cluster import KMeans, SpectralClustering
from sklearn.preprocessing import MinMaxScaler
from sklearn.decomposition import PCA
from sklearn.model_selection import ShuffleSplit
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LogisticRegression
from keras import optimizers
from keras.layers.core import Dense, Dropout, Activation, Flatten
from keras.utils import np_utils
from keras.models import Model
from keras.models import load_model, model_from_json
from keras import backend as K
import shutil
from scipy import stats
from sklearn import svm
from scipy.stats import norm
import csv
from sklearn.metrics import calinski_harabasz_score, davies_bouldin_score, silhouette_score
import time
from sklearn.model_selection import permutation_test_score, cross_val_score, cross_validate, KFold, StratifiedKFold
from Utils import plot_roc_curve, save_accuracy_data, plot_confusion_matrix
import community
import networkx as nx
import sys
import multiprocessing
import time
import gc
import heapq

SNP_NUMBER = 404078
CONTROL_NUMBER = 283


def freeMemory(Vari):
    del (Vari)
    gc.collect()


def getCommunitySnp(Communityfile):
    Communitydict = {}
    with open(Communityfile, 'r') as datafile:
        lines = datafile.read().strip().split("\n")
        for positionid in range(len(lines)):
            lst = lines[positionid].strip().split('\t')
            list_int = [eval(i) for i in lst]
            Communitydict[positionid] = list_int
    return Communitydict


def getPersonCommunity(Communitydict, PersonSnpfile):
    with open(PersonSnpfile) as inputfile:
        ComPerdict = {}
        for line in inputfile:
            parts = line.strip().split('\t')
            for Communityid in range(len(Communitydict)):
                if Communityid in ComPerdict:
                    list_singlesnp = []
                    for col in Communitydict[Communityid]:
                        if (col > SNP_NUMBER):
                            str2 = parts[col - SNP_NUMBER].strip().split('/')[1]
                            list_singlesnp.append(str2)
                        else:
                            str1 = parts[col].strip().split('/')[0]
                            list_singlesnp.append(str1)
                    ComPerdict[Communityid] = np.concatenate((ComPerdict[Communityid], list_singlesnp), axis=0)
                else:
                    list_singlesnp = []
                    for col in Communitydict[Communityid]:
                        if (col > SNP_NUMBER):
                            str2 = parts[col - SNP_NUMBER].strip().split('/')[1]
                            list_singlesnp.append(str2)
                        else:
                            str1 = parts[col].strip().split('/')[0]
                            list_singlesnp.append(str1)
                    ComPerdict[Communityid] = list_singlesnp
    return ComPerdict


def getOneHotSNP(dict_):
    one_hot = {
        'A': [0, 0, 0, 1],
        'T': [0, 0, 1, 0],
        '0': [0, 0, 0, 0],
        'C': [0, 1, 0, 0],
        'G': [1, 0, 0, 0]
    }
    onehotdic = {}
    for communityid in range(len(dict_)):
        if communityid in onehotdic:
            tmp = []
            for i in range(len(dict_[communityid])):
                tmp.extend(one_hot[dict_[communityid][i]])
            onehotdic[communityid] = np.concatenate((onehotdic[communityid], tmp), axis=0)
        else:
            tmp = []
            for i in range(len(dict_[communityid])):
                tmp.extend(one_hot[dict_[communityid][i]])
            onehotdic[communityid] = tmp
    return onehotdic


def getone_hot_community_control(one_hotdict):
    for i in range(len(one_hotdict)):
        one_hotdict[i] = np.array(one_hotdict[i]).reshape(CONTROL_NUMBER, -1)
    return one_hotdict


def getone_hot_community_case(one_hotdict, key_num):
    for i in range(len(one_hotdict)):
        one_hotdict[i] = np.array(one_hotdict[i]).reshape(key_num, -1)
    return one_hotdict



def _iter_edge_weights(edge_file):
    with open(edge_file, 'r') as f:
        for line in f:
            parts = line.strip().split('\t')
            if len(parts) < 3:
                continue
            try:
                yield float(parts[2])
            except Exception:
                continue


def compute_threshold_by_target_edge_count(edge_file, target_edges=SNP_NUMBER):

    heap = []
    total = 0

    for w in _iter_edge_weights(edge_file):
        total += 1
        if len(heap) < target_edges:
            heapq.heappush(heap, w)
        else:
            if w > heap[0]:
                heapq.heapreplace(heap, w)

    if total == 0:
        return None, 0

    if total < target_edges:
        return min(heap), total

    return heap[0], total


def get_or_compute_global_threshold(cache_path, reference_edge_file, target_edges=SNP_NUMBER):

    if cache_path and os.path.exists(cache_path):
        try:
            with open(cache_path, 'r') as f:
                val = float(f.read().strip())
            return val
        except Exception:
            pass 

    thr, total = compute_threshold_by_target_edge_count(reference_edge_file, target_edges=target_edges)
    if thr is None:
        return None

    if cache_path:
        os.makedirs(os.path.dirname(cache_path), exist_ok=True)
        with open(cache_path, 'w') as f:
            f.write(str(thr))

    return thr


def filter_by_weight_threshold(originsubnetfile, subnetfile, threshold):
    """weight >= threshold"""
    with open(originsubnetfile, 'r') as FromFile:
        with open(subnetfile, 'w') as OutFile:
            for line in FromFile:
                parts = line.strip().split('\t')
                if len(parts) < 3:
                    continue
                try:
                    weight = float(parts[2])
                except Exception:
                    continue
                if weight >= threshold:
                    OutFile.writelines(line)


#########################################################
# Get Community, one_hot
#########################################################

def getCommunity(netfile, communityfile):
    initnode = [-1 for i in range(808156)]
    G = nx.read_weighted_edgelist(netfile, delimiter='\t', nodetype=str, encoding='utf-8')
    partition = community.best_partition(G)
    communitylist = set(partition.values())
    cluster_num = len(communitylist)

    for nodeid, communityid in partition.items():
        initnode[int(nodeid) - 1] = int(communityid)
    line = " ".join(str(i) for i in initnode)
    with open(communityfile, 'w') as f:
        f.writelines("808156 nodes " + str(cluster_num) + " clusters 4444444 edges\n")
        f.writelines(line)


#########################################################
# classifier
#########################################################

def knn_test(p, q, count, n_splits, test_size):
    print("====knn====")
    n_neighbors = int((p.shape[0] * (1 - test_size)) ** 0.5)
    sklearn_knn_clf = KNeighborsClassifier(n_neighbors=n_neighbors)
    cv = ShuffleSplit(n_splits=n_splits, test_size=test_size, random_state=0)
    scores = cross_val_score(sklearn_knn_clf, p, q, cv=cv, scoring='accuracy')
    scores = sorted(scores)[1:-2]
    score = np.mean(scores)
    std = (0.25 / p.shape[0] * test_size) ** 0.5
    mean = 0.5
    normalDistribution = norm(mean, std)
    p_value = 1 - normalDistribution.cdf(score)
    return score, p_value


def rf_test(p, q, count, n_splits, test_size):
    print("====rf====")
    from sklearn.ensemble import RandomForestClassifier
    n_estimators = 200
    max_depth = 15
    min_samples_leaf = 2
    rfc = RandomForestClassifier(
        n_estimators=n_estimators,
        max_depth=max_depth,
        min_samples_leaf=min_samples_leaf,
        random_state=90
    )
    cv = ShuffleSplit(n_splits=n_splits, test_size=test_size, random_state=0)
    scores = cross_val_score(rfc, p, q, cv=cv, scoring='accuracy')
    scores = sorted(scores)[1:-2]
    score = np.mean(scores)
    std = (0.25 / p.shape[0] * test_size) ** 0.5
    mean = 0.5
    normalDistribution = norm(mean, std)
    p_value = 1 - normalDistribution.cdf(score)
    return score, p_value


def svm_test(datanumpy, labellist):
    label = np.array(labellist)
    data = np.mat(datanumpy)

    kernel = 'rbf'
    if (data.shape[0] < data.shape[1]):
        kernel = 'linear'
    clf = svm.SVC(
        C=100.,
        kernel=kernel,
        decision_function_shape='ovr',
        max_iter=1000,
        gamma='auto',
        class_weight='balanced',
        probability=True
    )

    KF = StratifiedKFold(n_splits=5, random_state=0, shuffle=True)

    scores = []
    yPredList = []
    yTrueList = []

    for train_index, test_index in KF.split(data, label):
        X_train, X_test = data[train_index], data[test_index]
        Y_train, Y_test = label[train_index], label[test_index]

        clf.fit(X_train, Y_train)

        y_pred = clf.predict_proba(X_test)
        meanScore = clf.score(X_test, Y_test)
        scores.append(meanScore)
        yPredList.append(y_pred)
        yTrueList.append(Y_test)

    return yPredList, yTrueList, np.mean(scores)


def lr_test(p, q, count, n_splits, test_size):
    print("====lr====")
    logreg = LogisticRegression(solver='liblinear')
    cv = ShuffleSplit(n_splits=n_splits, test_size=test_size, random_state=0)
    scores = cross_val_score(logreg, p, q, cv=cv, scoring='accuracy')
    scores = sorted(scores)[1:-2]
    score = np.mean(scores)
    std = (0.25 / p.shape[0] * test_size) ** 0.5
    mean = 0.5
    normalDistribution = norm(mean, std)
    p_value = 1 - normalDistribution.cdf(score)
    print(scores)
    return score, p_value


def c2st(cluter_num, data_first, data_second, label_first, label_second,
         community_num, first_cluster_num, second_cluster_num=-1,
         num_train=50, base_fileName='/res/', epochs=50, batch_size=16,
         test_size=0.1, classifier_type='svm'):

    length = min(len(label_first), len(label_second))
    print("__________________________start_________________________")
    acc_dict = {}
    p_value_dict = {}
    index = 0
    community_list = range(0, community_num)
    some_accuracy_data = pd.DataFrame(columns=['accuracy', 'F1', 'precession', 'recall'], index=community_list)

    for key, value in data_first.items():
        health_npy = value
        patient_npy = data_second.get(key)
        health_npy = health_npy.reshape((len(label_first), np.size(health_npy) // len(label_first)))
        patient_npy = patient_npy.reshape((len(label_second), np.size(patient_npy) // len(label_second)))

        np.random.seed(7)
        row_rand_array = np.arange(health_npy.shape[0])
        np.random.shuffle(row_rand_array)
        health_value = health_npy[row_rand_array[0:length]]
        row_rand_array = np.arange(patient_npy.shape[0])
        np.random.shuffle(row_rand_array)
        patient_value = patient_npy[row_rand_array[0:length]]

        x_data = np.concatenate((health_value, patient_value), axis=0)
        y_data = random.sample(label_first, length) + random.sample(label_second, length)

        index += 1
        acc = 0
        p_value = 0

        if (classifier_type == "knn"):
            acc, p_value = knn_test(x_data, y_data, index, num_train, test_size=test_size)
        elif (classifier_type == "svm"):
            yPredList, yTrueList, acc = svm_test(x_data, y_data)
            save_accuracy_data(yTrueList, yPredList, key, some_accuracy_data)
            if acc > 0.58:
                plot_roc_curve(
                    yTrueList, yPredList, 2,
                    save_path=base_fileName + 'roc/',
                    file_name="{}_{}_roc.jpg".format(str(first_cluster_num), str(key)),
                    image_flag=True
                )
                plot_confusion_matrix(
                    yTrueList, yPredList, 2,
                    save_path=base_fileName + 'confusion_matrix/',
                    file_name="{}_{}_matrix.jpg".format(str(first_cluster_num), str(key))
                )
        elif (classifier_type == "lr"):
            acc, p_value = lr_test(x_data, y_data, index, num_train, test_size=test_size)
        elif (classifier_type == "rf"):
            acc, p_value = rf_test(x_data, y_data, index, num_train, test_size=test_size)

        acc_dict[key] = acc
        p_value_dict[key] = 1 - p_value

    acc_fileName = base_fileName + 'c2st_accuracy/'
    dict_to_csv(cluter_num, acc_fileName, acc_dict, first_cluster_num, second_cluster_num)

    some_accuracy_datapath = base_fileName + 'some_accuracy_data/'
    if not os.path.exists(some_accuracy_datapath):
        os.makedirs(some_accuracy_datapath)
    fileName = some_accuracy_datapath + '{}.xlsx'.format(first_cluster_num)
    writer = pd.ExcelWriter(fileName)
    some_accuracy_data.to_excel(writer, index=True)
    writer.save()

    return acc_dict, p_value_dict


def dict_to_csv(cluter_num, path, datas, first_cluster_num, second_cluster_num):
    if not os.path.exists(path):
        os.makedirs(path)
    if (second_cluster_num == -1):
        filename = path + 'cluster_{}_class{}.csv'.format(str(cluter_num), str(first_cluster_num))
    else:
        filename = path + 'cluster_{}_class{}_{}.csv'.format(str(cluter_num), str(first_cluster_num), str(second_cluster_num))
    with open(filename, 'a', newline='') as f:
        writer = csv.writer(f)
        for key, value in datas.items():
            writer.writerow([key, value])


def getclusterdata(key, column, subname):
    if os.path.exists('bloc/%s/%d/%s_cluster_genotype_%d.txt' % (subname, key, subname, key)):
        return
    with open('bloc/genotype_sz_case_QC_exchange.txt', 'r') as openfile:
        with open('bloc/%s/%d/%s_cluster_genotype_%d.txt' % (subname, key, subname, key), 'w') as inputfile:
            index = 0
            for line in openfile:
                if index in column:
                    index = index + 1
                    inputfile.writelines(line)
                else:
                    index = index + 1
                    continue
    freeMemory(inputfile)
    freeMemory(openfile)


def add(num, subname, key, key_num, threshold):
    startnum = 15000 * num + 1
    mid = 15000 * (num + 1)
    end = startnum + 1
    print(num, startnum, mid, end)
    inputfile = 'bloc/%s/%d/%s_cluster_genotype_%d.txt' % (subname, key, subname, key)

    subnetfile = 'bloc/%s/%d/net/%d.txt' % (subname, key, num)
    cmd = "./ccc %s %s %f %d SNP_NUMBER 0 1 %d %d %d SNP_NUMBER" % (inputfile, subnetfile, threshold, key_num, startnum, mid, end)
    os.system(cmd)
    os.system("exit")
    freeMemory(inputfile)
    freeMemory(subnetfile)


def mergefile(path, filename):
    filedir = path
    filenames = os.listdir(filedir)
    f = open(filename, 'w')
    for fn in filenames:
        filepath = filedir + '/' + fn
        for line in open(filepath):
            f.writelines(line)
    f.close()
    shutil.rmtree(path)


def accuracy(resultor_communityInfo, inputfile, key_num, savepath, key):
    community_data = getCommunitySnp(resultor_communityInfo)

    community_num = len(community_data.keys())
    dict_case = getPersonCommunity(community_data, inputfile)
    dict_control = getPersonCommunity(community_data, 'singlesnp/genotype_sz_control_QC_exchange.txt')

    one_hotdict_case = getOneHotSNP(dict_case)
    one_hotdict_control = getOneHotSNP(dict_control)
    data_healthy = getone_hot_community_control(one_hotdict_control)
    data_patient = getone_hot_community_case(one_hotdict_case, key_num)

    label_healthy = [0 for _ in range(0, CONTROL_NUMBER)]
    label_patient = [1 for _ in range(0, key_num)]

    basefile = savepath + '/res/'
    if not os.path.exists(basefile):
        os.makedirs(basefile)
    c2st(1, data_healthy, data_patient, label_healthy, label_patient,
         community_num, key, second_cluster_num=-1, num_train=5,
         base_fileName=basefile, epochs=50)


def communityBycarriers(subname, subnetfile, inputfile, key_num, key, resultor, resultor_communityInfo, threshold):
    subnetcommunity = 'bloc/%s/%d/%d_cluster_community_%f.bfs' % (subname, key, key, threshold)
    getCommunity(subnetfile, subnetcommunity)
    cmd_carries = "./carriers %s %s bloc/genotype_sz_control_QC_exchange.txt 0 1 %d CONTROL_NUMBER SNP_NUMBER %s %s" % (
        subnetcommunity, inputfile, key_num, resultor, resultor_communityInfo
    )
    os.system(cmd_carries)


def ccc(subname, key, key_num, inputfile, output_netfile, threshold):
    cores = multiprocessing.cpu_count()
    pool = multiprocessing.Pool(processes=cores)
    results = []
    path = 'bloc/%s/%d/net/' % (subname, key)
    if not os.path.exists(path):
        os.makedirs(path)

    print(os.getcwd())
    for i in range(0, 26):
        r = pool.apply_async(add, args=(i, subname, key, key_num, threshold,))
        results.append(r)
    _ = [xx.get() for xx in results]
    pool.close()
    pool.join()

    # CCC
    finalnetfile = 'bloc/%s/%d/net/27.txt' % (subname, key)
    cmd = "./ccc %s %s %f %d SNP_NUMBER 0 1 390000 SNP_NUMBER 390001 SNP_NUMBER" % (inputfile, finalnetfile, threshold, key_num)
    os.system(cmd)

    # merge
    mergefile("bloc/%s/%d/net" % (subname, key), output_netfile)
    freeMemory(inputfile)
    freeMemory(finalnetfile)


def getAccuracyByCommunity(cluster, threshold, savepath, subname,
                           reference_edge_file=None,
                           SNP_NUMBER=SNP_NUMBER):
                           
    # 1) global threshold*
    thr_cache_file = os.path.join('bloc', subname, f'global_threshold_top{SNP_NUMBER}.txt')
    global_threshold = None

    if threshold is not None:
        global_threshold = float(threshold)
    else:
        if reference_edge_file is None:
            guess = os.path.join('bloc', subname, 'origin_network.txt')
            if os.path.exists(guess):
                reference_edge_file = guess

        if reference_edge_file is not None and os.path.exists(reference_edge_file):
            global_threshold = get_or_compute_global_threshold(
                thr_cache_file, reference_edge_file, target_edges=SNP_NUMBER
            )
            print(f"threshold={global_threshold}")
        else:
            global_threshold = None
            print("reference_edge_file not found!")

    y_list = cluster.reshape(-1).tolist()
    dict_idx = {}
    for index in range(len(y_list)):
        if y_list[index] in dict_idx:
            dict_idx[y_list[index]].append(index)
        else:
            dict_idx[y_list[index]] = [index]

    for key in sorted(dict_idx.keys()):
        subtypenetworkdir = 'bloc/%s/%d/' % (subname, key)
        if not os.path.exists(subtypenetworkdir):
            os.makedirs(subtypenetworkdir)

        key_num = len(dict_idx[key])

        # genotype SNP
        inputfile = 'bloc/%s/%d/%s_cluster_genotype_%d.txt' % (subname, key, subname, key)
        if not os.path.exists(inputfile):
            getclusterdata(key, dict_idx[key], subname)

        originsubnetfile = 'bloc/%s/%d/%d_cluster_network.txt' % (subname, key, key)

        if not os.path.exists(originsubnetfile):
            raw_threshold = 7.0
            print(f"{originsubnetfile} not found -> run ccc with threshold={raw_threshold}")
            ccc(subname, key, key_num, inputfile, originsubnetfile, raw_threshold)

        if global_threshold is None:
            global_threshold = get_or_compute_global_threshold(
                thr_cache_file, originsubnetfile, target_edges=SNP_NUMBER
            )

        subnetfile = 'bloc/%s/%d/%d_cluster_network_%f.txt' % (subname, key, key, global_threshold)
        if not os.path.exists(subnetfile):
            filter_by_weight_threshold(originsubnetfile, subnetfile, global_threshold)

        resultor = 'bloc/%s/%d/%d_cluster_result_%f.xlsx' % (subname, key, key, global_threshold)
        resultor_communityInfo = 'bloc/%s/%d/%d_cluster_communityid_%f.txt' % (subname, key, key, global_threshold)

        if not os.path.exists(resultor_communityInfo):
            communityBycarriers(subname, subnetfile, inputfile, key_num, key, resultor, resultor_communityInfo, global_threshold)

        if not os.path.exists(savepath + '/res/some_accuracy_data/%s.xlsx' % (key)):
            accuracy(resultor_communityInfo, inputfile, key_num, savepath, key)

        freeMemory(inputfile)
        freeMemory(subnetfile)
        freeMemory(originsubnetfile)