import pandas as pd
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import RFE
from sklearn.model_selection import StratifiedKFold, cross_val_predict
from sklearn.metrics import roc_auc_score, roc_curve, accuracy_score, confusion_matrix
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib.lines import Line2D
# Load the 111 peptide dataset
file_path = '/Users/rachitbisht/Desktop/research/data/BC_Machine_Learning/datasets/Adjusted_062223_MH_Data_For_Training.xlsx'
data = pd.read_excel(file_path)
def get_feature_type(feature_name):
"""Determine if the feature is Acetyl, Glyco, Phospho, or Other."""
if 'Acetyl' in feature_name:
return 'Acetyl'
elif 'Glyco' in feature_name:
return 'Glyco'
elif 'Phospho' in feature_name:
return 'Phospho'
else:
return 'Other'
def perform_analysis(data, comparison_groups, target_group, n_features_to_select, exclude_features=[]):
# Adjust binary labeling: '1' for target group, '0' for comparison groups
data['Binary_Label'] = data['Type'].apply(lambda x: 1 if x == target_group else 0 if x in comparison_groups else None)
# Filter data to include only the relevant groups
data_filtered = data[data['Type'].isin(comparison_groups + [target_group])].copy()
# Prepare the data and exclude specified features
X = data_filtered.drop(['Type', 'Binary_Label'] + exclude_features, axis=1)
y = data_filtered['Binary_Label']
# Initialize the Random Forest classifier with class_weight='balanced'
rf_clf = RandomForestClassifier(max_depth=4, n_estimators=300, min_samples_split=2,min_samples_leaf=3, max_features='sqrt', random_state=42, class_weight='balanced')
#LAB 0.83 auc
# rf_clf = RandomForestClassifier(n_estimators=300, min_samples_split=2,min_samples_leaf=1, max_features='sqrt', random_state=42, class_weight='balanced')
# Use RFE to select the top features
rfe = RFE(estimator=rf_clf, n_features_to_select=n_features_to_select, step=1)
rfe.fit(X, y)
# Extract the selected feature names
selected_features = X.columns[rfe.support_]
print("Selected features:")
for feature in selected_features:
print(f"{feature}")
# Filter the features to the ones selected by RFE and refit RandomForest on these features
X_rfe = X.iloc[:, rfe.support_]
rf_clf.fit(X_rfe, y) # Refitting the model on the selected features
# Setup stratified cross-validation
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
y_prob = cross_val_predict(rf_clf, X_rfe, y, cv=cv, method='predict_proba')
# Compute AUC score
auc_score = roc_auc_score(y, y_prob[:, 1])
# Compute ROC curve
fpr, tpr, thresholds = roc_curve(y, y_prob[:, 1])
# Compute accuracy
y_pred = cross_val_predict(rf_clf, X_rfe, y, cv=cv)
accuracy = accuracy_score(y, y_pred)
# Compute confusion matrix
cm = confusion_matrix(y, y_pred)
# Plot the ROC curve
plt.figure()
plt.figure(figsize=(6, 6)) # Square plot
plt.plot(fpr, tpr, color='blue', lw=2, label=f'ROC curve (AUC = {auc_score:.2f})')
plt.plot([0, 1], [0, 1], color='gray', lw=2, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.0])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC Curve')
plt.legend(loc='lower right')
plt.show()
# Print and plot the confusion matrix
print("Confusion Matrix:")
print(cm)
# Plotting confusion matrix as a heatmap
plt.figure(figsize=(6, 6))
sns.heatmap(cm, annot=True, fmt="d", cmap="Blues", cbar=False, xticklabels=["Rest", "TN"], yticklabels=["Rest", "TN"])
plt.ylabel('True Label')
plt.xlabel('Predicted Label')
plt.title('Confusion Matrix')
plt.show()
# VIP (dot plot for feature importance)
feature_importances = rf_clf.feature_importances_
# Sort importances and features (descending order to have most important at the top)
sorted_idx = feature_importances.argsort()[::-1]
sorted_importances = feature_importances[sorted_idx]
sorted_features = selected_features[sorted_idx]
# Get the color for each feature based on its type
colors = ['red' if get_feature_type(feature) == 'Acetyl'
else 'green' if get_feature_type(feature) == 'Glyco'
else 'blue' if get_feature_type(feature) == 'Phospho'
else 'gray' for feature in sorted_features]
# Print feature importance scores
print("\nFeature importance scores:")
for feature, importance in zip(sorted_features, sorted_importances):
print(f"{feature}: {importance:.4f}")
# Plot the VIP as a dot plot with colors based on the feature type
plt.figure(figsize=(10, 6))
plt.scatter(sorted_importances, sorted_features, color=colors, s=100) # Scatter plot (dot plot)
plt.xlabel('Importance')
plt.title('Variable Importance Plot')
# Remove the background and turn off grid lines
plt.gca().patch.set_visible(False) # Removes the background
plt.grid(False) # Turns off grid lines
plt.gca().invert_yaxis() # Invert y-axis to have the most important feature at the top
# Add a legend in the bottom-right corner
legend_elements = [
Line2D([0], [0], marker='o', color='w', label='Acetyl', markerfacecolor='red', markersize=10),
Line2D([0], [0], marker='o', color='w', label='Glyco', markerfacecolor='green', markersize=10),
Line2D([0], [0], marker='o', color='w', label='Phospho', markerfacecolor='blue', markersize=10)
]
plt.legend(handles=legend_elements, loc='lower right', title='Feature Type')
plt.tight_layout() # Automatically adjusts the plot to make space for the labels
plt.show()
return auc_score, accuracy, cm
# features to exclude due to target leakage (not needed)
exclude_features = []
# Perform the analysis for LAB and Control vs TN excluding the specified features
# if running LAB experiment use ['TN', 'Control'], 'LAB' for parameters
# else for TN experiment use ['LAB', 'Control'], 'TN'
auc_score, accuracy, cm = perform_analysis(data, ['TN', 'Control'], 'LAB', 30, exclude_features)
print(f"AUC score: {auc_score:.2f}")
print(f"Accuracy: {accuracy:.2f}")