#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""

Created on Sun Jun  7 14:33:03 2020

#Classifier to distinguish between patients with a pouch phenotypes (normal pouch vs. pouchitis) based on bacterial speceis, metabolic pathways or enzmes profiles.
#Author: Vadim (Dani) Dubinsky (<EMAIL>)

"""
import sys
import pandas as pd
import numpy as np
from numpy import interp #returns the one-dimensional piecewise linear interpolant to a function with given discrete data points
import matplotlib.pyplot as plt
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import confusion_matrix, roc_curve, auc

try:
    import xgboost as xgb ##XGBoost, extreme gradient boosting model
except:
    sys.exit("This program requires Python3 xgboost module, please install it and try again")


#Initialize lists to hold the results obtained from running the code below for repeated cross validation
mean_tpr_l, tprs_lower_l, tprs_upper_l, mean_auc_l, std_auc_l, mean_fpr_l = [],[],[],[],[],[]

#Repeat the 5-fold cross validation x100 times, to generate results based on average ROC curves (over folds and over repetitions) for plotting:
mean_fpr = np.linspace(0, 1, 100)
tprs, aucs = [],[]
i = 0
for _ in range(50):
    i += 1
    kfold = StratifiedKFold(n_splits=5, random_state=None, shuffle=True)
    #Split the training set into k-folds for cross validation, and calculate all the parameters for each fold repeatedly (x100)
    for train, test in kfold.split(taxa_data_sub.values, target_labels["labels_binary"].values):
        probs_cv = model.fit(taxa_data_sub.values[train], target_labels["labels_binary"].values[train]).predict_proba(taxa_data_sub.values[test])
        # Compute ROC curve and area the curve
        fpr, tpr, thresholds = roc_curve(target_labels["labels_binary"].values[test], probs_cv[:, 1])
        tprs.append(interp(mean_fpr, fpr, tpr))
        tprs[-1][0] = 0.0
        roc_auc = auc(fpr, tpr)
        aucs.append(roc_auc)

#Calculate all the needed variables to plot a mean ROC AUC with standard deviation area
mean_tpr, mean_auc, std_auc, tprs_upper, tprs_lower = calc_vars_AUC_std (tprs, aucs)    

#Plot multiple mean ROC AUC with standard deviation area
plot_mean_ROC_AUC_std_multiple_models(mean_fpr_l, mean_tpr_l, tprs_lower_l,
                                      tprs_upper_l, mean_auc_l, std_auc_l)


#----------------------Utility functions:-----------------------------------###

def calc_vars_AUC_std (tprs, aucs):
    '''Function to calculate all the needed variables to plot a ROC AUC curve with STD
       input arguments: tprs, auc (list) and mean_tpr (numpy array)'''
    mean_tpr = np.mean(tprs, axis=0)
    mean_tpr[-1] = 1.0
    mean_auc = np.mean(aucs)
    std_auc = np.std(aucs)
    std_tpr = np.std(tprs, axis=0)
    tprs_upper = np.minimum(mean_tpr + std_tpr, 1)
    tprs_lower = np.maximum(mean_tpr - std_tpr, 0)
    
    mean_tpr_l.append(mean_tpr)
    mean_auc_l.append(mean_auc)
    std_auc_l.append(std_auc)
    tprs_upper_l.append(tprs_upper)
    tprs_lower_l.append(tprs_lower)
    mean_fpr_l.append(mean_fpr)
    return (mean_tpr, mean_auc, std_auc, tprs_upper, tprs_lower)

###-------------------------------------------------------------------------###    

def plot_mean_ROC_AUC_std_multiple_models (mean_fpr_l, mean_tpr_l, tprs_lower_l, tprs_upper_l, mean_auc_l, std_auc_l):
    '''Function to plot ROC AUC curves with std intervals and with several models (1-3 AUCs). Input arguments are the output from calc_vars_AUC_std() function.
       You need to manually define the colors and the labels for each model ROC AUC. The number of curves that will be plotted depends on len(mean_auc_l) - number of supplied models. Thus the labels/colors needs to be changed accordingly
       With mean area under curve, we can see the variance of the curve when the training set is split into different subsets. This can show how the classifier output is affected by changes in the training data, and how different the splits generated by K-fold cross-validation are from one another.'''
    colors = ["darkviolet", "darkgreen", "darkred"] #"navy"
    labels = ["GBT species, full set", "GBT species, top k=40", "GBT species, top k=40 + calpro"]
    #colors = ["darkorange"] #for calprotectin as sole predictor
    #labels = ["calprotectin only"] #for calprotectin as sole predictor
    
    plt.rcParams["font.weight"] = "bold"
    plt.rcParams["axes.labelweight"] = "bold"
    plt.figure(figsize=(10, 10))   

    plt.plot([0, 1], [0, 1], linestyle="--", lw=2, color="darkblue", label="Random sample", alpha=.8)
    
    for i in range(len(mean_auc_l)):
        plt.plot(mean_fpr_l[i], mean_tpr_l[i], color = colors[i], label = labels[i] + r" (AUC = %0.3f $\pm$ %0.3f)" % (mean_auc_l[i], std_auc_l[i]), lw=3)
        plt.fill_between(mean_fpr_l[i], tprs_lower_l[i], tprs_upper_l[i], color = colors[i], alpha=.1) # label=r"$\pm$ 1 Standard Deviation"

    plt.xlim([-0.01, 1.01])
    plt.ylim([-0.01, 1.01])
    plt.xlabel("False Positive Rate", fontweight = "bold" , fontsize=18)
    plt.ylabel("True Positive Rate", fontweight = "bold" , fontsize=18)
    plt.tick_params(axis="both", which="major", labelsize=16)
    plt.legend( prop={"size":13} , loc = "lower right")
    plt.title("Classifier for Pouchitis vs. Normal pouch (species)", fontsize=20, pad=10)
    #plt.savefig("/home/.../xxx.png", format='png', dpi=350, bbox_inches='tight')

###-------------------------------------------------------------------------###

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""

Created on Sun Jun  7 14:33:03 2020

#Classifier to distinguish between patients with a pouch phenotypes (normal pouch vs. pouchitis) based on bacterial speceis, metabolic pathways or enzmes profiles.
#Author: Vadim (Dani) Dubinsky (<EMAIL>)

"""
import sys
import pandas as pd
import numpy as np
from numpy import interp #returns the one-dimensional piecewise linear interpolant to a function with given discrete data points
import matplotlib.pyplot as plt
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import confusion_matrix, roc_curve, auc

try:
    import xgboost as xgb ##XGBoost, extreme gradient boosting model
except:
    sys.exit("This program requires Python3 xgboost module, please install it and try again")


#Initialize lists to hold the results obtained from running the code below for repeated cross validation
mean_tpr_l, tprs_lower_l, tprs_upper_l, mean_auc_l, std_auc_l, mean_fpr_l = [],[],[],[],[],[]

#Repeat the 5-fold cross validation x100 times, to generate results based on average ROC curves (over folds and over repetitions) for plotting:
mean_fpr = np.linspace(0, 1, 100)
tprs, aucs = [],[]
i = 0
for _ in range(50):
    i += 1
    kfold = StratifiedKFold(n_splits=5, random_state=None, shuffle=True)
    #Split the training set into k-folds for cross validation, and calculate all the parameters for each fold repeatedly (x100)
    for train, test in kfold.split(taxa_data_sub.values, target_labels["labels_binary"].values):
        probs_cv = model.fit(taxa_data_sub.values[train], target_labels["labels_binary"].values[train]).predict_proba(taxa_data_sub.values[test])
        # Compute ROC curve and area the curve
        fpr, tpr, thresholds = roc_curve(target_labels["labels_binary"].values[test], probs_cv[:, 1])
        tprs.append(interp(mean_fpr, fpr, tpr))
        tprs[-1][0] = 0.0
        roc_auc = auc(fpr, tpr)
        aucs.append(roc_auc)

#Calculate all the needed variables to plot a mean ROC AUC with standard deviation area
mean_tpr, mean_auc, std_auc, tprs_upper, tprs_lower = calc_vars_AUC_std (tprs, aucs)    

#Plot multiple mean ROC AUC with standard deviation area
plot_mean_ROC_AUC_std_multiple_models(mean_fpr_l, mean_tpr_l, tprs_lower_l,
                                      tprs_upper_l, mean_auc_l, std_auc_l)


#----------------------Utility functions:-----------------------------------###

def calc_vars_AUC_std (tprs, aucs):
    '''Function to calculate all the needed variables to plot a ROC AUC curve with STD
       input arguments: tprs, auc (list) and mean_tpr (numpy array)'''
    mean_tpr = np.mean(tprs, axis=0)
    mean_tpr[-1] = 1.0
    mean_auc = np.mean(aucs)
    std_auc = np.std(aucs)
    std_tpr = np.std(tprs, axis=0)
    tprs_upper = np.minimum(mean_tpr + std_tpr, 1)
    tprs_lower = np.maximum(mean_tpr - std_tpr, 0)
    
    mean_tpr_l.append(mean_tpr)
    mean_auc_l.append(mean_auc)
    std_auc_l.append(std_auc)
    tprs_upper_l.append(tprs_upper)
    tprs_lower_l.append(tprs_lower)
    mean_fpr_l.append(mean_fpr)
    return (mean_tpr, mean_auc, std_auc, tprs_upper, tprs_lower)

###-------------------------------------------------------------------------###    

def plot_mean_ROC_AUC_std_multiple_models (mean_fpr_l, mean_tpr_l, tprs_lower_l, tprs_upper_l, mean_auc_l, std_auc_l):
    '''Function to plot ROC AUC curves with std intervals and with several models (1-3 AUCs). Input arguments are the output from calc_vars_AUC_std() function.
       You need to manually define the colors and the labels for each model ROC AUC. The number of curves that will be plotted depends on len(mean_auc_l) - number of supplied models. Thus the labels/colors needs to be changed accordingly
       With mean area under curve, we can see the variance of the curve when the training set is split into different subsets. This can show how the classifier output is affected by changes in the training data, and how different the splits generated by K-fold cross-validation are from one another.'''
    colors = ["darkviolet", "darkgreen", "darkred"] #"navy"
    labels = ["GBT species, full set", "GBT species, top k=40", "GBT species, top k=40 + calpro"]
    #colors = ["darkorange"] #for calprotectin as sole predictor
    #labels = ["calprotectin only"] #for calprotectin as sole predictor
    
    plt.rcParams["font.weight"] = "bold"
    plt.rcParams["axes.labelweight"] = "bold"
    plt.figure(figsize=(10, 10))   

    plt.plot([0, 1], [0, 1], linestyle="--", lw=2, color="darkblue", label="Random sample", alpha=.8)
    
    for i in range(len(mean_auc_l)):
        plt.plot(mean_fpr_l[i], mean_tpr_l[i], color = colors[i], label = labels[i] + r" (AUC = %0.3f $\pm$ %0.3f)" % (mean_auc_l[i], std_auc_l[i]), lw=3)
        plt.fill_between(mean_fpr_l[i], tprs_lower_l[i], tprs_upper_l[i], color = colors[i], alpha=.1) # label=r"$\pm$ 1 Standard Deviation"

    plt.xlim([-0.01, 1.01])
    plt.ylim([-0.01, 1.01])
    plt.xlabel("False Positive Rate", fontweight = "bold" , fontsize=18)
    plt.ylabel("True Positive Rate", fontweight = "bold" , fontsize=18)
    plt.tick_params(axis="both", which="major", labelsize=16)
    plt.legend( prop={"size":13} , loc = "lower right")
    plt.title("Classifier for Pouchitis vs. Normal pouch (species)", fontsize=20, pad=10)
    #plt.savefig("/home/.../xxx.png", format='png', dpi=350, bbox_inches='tight')

###-------------------------------------------------------------------------###

#!/usr/bin/env python3 # -*- coding: utf-8 -*- Created on Sun Jun  7 14:33:03 2020  #Classifier to distinguish between patients with a pouch phenotypes (normal pouch vs. pouchitis) based on bacterial speceis, metabolic pathways or enzmes profiles. #Author: Vadim (Dani) Dubinsky (<EMAIL>) #returns the one-dimensional piecewise linear interpolant to a function with given discrete data points ##XGBoost, extreme gradient boosting model #Initialize lists to hold the results obtained from running the code below for repeated cross validation #Repeat the 5-fold cross validation x100 times, to generate results based on average ROC curves (over folds and over repetitions) for plotting: #Split the training set into k-folds for cross validation, and calculate all the parameters for each fold repeatedly (x100) # Compute ROC curve and area the curve #Calculate all the needed variables to plot a mean ROC AUC with standard deviation area #Plot multiple mean ROC AUC with standard deviation area #----------------------Utility functions:-----------------------------------### Function to calculate all the needed variables to plot a ROC AUC curve with STD        input arguments: tprs, auc (list) and mean_tpr (numpy array) ###-------------------------------------------------------------------------### Function to plot ROC AUC curves with std intervals and with several models (1-3 AUCs). Input arguments are the output from calc_vars_AUC_std() function.        You need to manually define the colors and the labels for each model ROC AUC. The number of curves that will be plotted depends on len(mean_auc_l) - number of supplied models. Thus the labels/colors needs to be changed accordingly        With mean area under curve, we can see the variance of the curve when the training set is split into different subsets. This can show how the classifier output is affected by changes in the training data, and how different the splits generated by K-fold cross-validation are from one another. #"navy" #colors = ["darkorange"] #for calprotectin as sole predictor #labels = ["calprotectin only"] #for calprotectin as sole predictor # label=r"$\pm$ 1 Standard Deviation" #plt.savefig("/home/.../xxx.png", format='png', dpi=350, bbox_inches='tight') ###-------------------------------------------------------------------------###