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}")