gliomas.py

# -*- coding: utf-8 -*-
"""Gliomas.ipynb

Automatically generated by Colab.

Original file is located at
    https://colab.research.google.com/drive/1urqNeoAGS4Gz-EGULpcFlTOyNUXFjgzA
"""

import pandas as pd
from sklearn.preprocessing import LabelEncoder, StandardScaler
from sklearn.model_selection import train_test_split,  cross_val_score
from sklearn.svm import SVC
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import classification_report, accuracy_score
from sklearn.model_selection import learning_curve
from sklearn.model_selection import GridSearchCV
from xgboost import XGBClassifier
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

# Load the data
file_path = '/content/TCGA_GBM_LGG_Mutations_all.csv'
df = pd.read_csv(file_path)

# Columns to be encoded
columns_to_encode = ['Grade', 'Project', 'Case_ID', 'Gender',
       'Primary_Diagnosis', 'Race', 'IDH1', 'TP53', 'ATRX', 'PTEN', 'EGFR',
       'CIC', 'MUC16', 'PIK3CA', 'NF1', 'PIK3R1', 'FUBP1', 'RB1', 'NOTCH1',
       'BCOR', 'CSMD3', 'SMARCA4', 'GRIN2A', 'IDH2', 'FAT4', 'PDGFRA']

# Initialize LabelEncoder
label_encoder = LabelEncoder()

# Create a DataFrame to store encoded data
encoded_df = pd.DataFrame()

# Encode categorical columns and store them in encoded_df
for col in columns_to_encode:
    encoded_df[col + '_encoded'] = label_encoder.fit_transform(df[col])

# Handle 'Age_at_diagnosis' separately
df['Age_at_diagnosis'].fillna('-1 years', inplace=True)
df['Age_at_diagnosis'] = df['Age_at_diagnosis'].astype(str).str.extract(r'(\d+)').fillna(-1).astype(int)
# Add the processed 'Age_at_diagnosis' to encoded_df
encoded_df['Age_at_diagnosis'] = df['Age_at_diagnosis']

# Drop rows with NaN or missing values in encoded_df
encoded_df.dropna(inplace=True)

# Display the DataFrame with the extracted age values
print("DataFrame with Extracted Age Values:")
print(encoded_df)

print(encoded_df.columns)

# List class names for Gender
gender_classes = label_encoder.fit(df['Gender']).classes_
print("Class names for Gender:", gender_classes)

# Drop rows where Gender is '--' in encoded_df
gender_encoded_class_index = list(gender_classes).index('--')
encoded_df = encoded_df[encoded_df['Gender_encoded'] != gender_encoded_class_index]

# List class names for Gender
gender_classes = label_encoder.fit(encoded_df['Gender_encoded']).classes_
print("Class names for Gender:", gender_classes)

# Features (all columns except 'Grade_encoded')
features = encoded_df.drop(columns=['Grade_encoded'])

# Label ('Grade_encoded')
label = encoded_df['Grade_encoded']

# Plot bar plots for features
for col in features.columns:
    plt.figure(figsize=(8, 6))
    sns.countplot(x=col, data=features)
    plt.title(f'Bar Plot for {col}')
    plt.xlabel(col)
    plt.ylabel('Count')
    plt.tight_layout()
    plt.show()

# Plot bar plot for the label
plt.figure(figsize=(8, 6))
sns.countplot(x=label)
plt.title('Bar Plot for Label (Grade_encoded)')
plt.xlabel('Grade_encoded')
plt.ylabel('Count')
plt.tight_layout()
plt.show()

# Compute the correlation matrix
correlation_matrix = encoded_df.corr()

# Print the correlation matrix
print("\nCorrelation Matrix:")
print(correlation_matrix)

# Print the numerical values in the correlation matrix greater than 0.80
print("\nCorrelation Values Greater Than 0.80 (excluding diagonal elements):")
for i in range(len(correlation_matrix)):
    for j in range(len(correlation_matrix.columns)):
        if i != j:  # Exclude diagonal elements
            correlation_value = correlation_matrix.iloc[i, j]
            if correlation_value > 0.80:
                print(f"Correlation between {correlation_matrix.index[i]} and {correlation_matrix.columns[j]}: {correlation_value}")

# Plot the heatmap
plt.figure(figsize=(12, 10))
sns.heatmap(correlation_matrix, cmap='coolwarm', annot=False, fmt=".2f")
plt.title('Correlation Heatmap of Encoded DataFrame')
plt.show()

print(encoded_df.columns)

# Drop 'Case_ID_encoded' and 'Project_encoded' from encoded_df
encoded_df.drop(columns=['Case_ID_encoded', 'Project_encoded', 'Primary_Diagnosis_encoded'], inplace=True)

print(encoded_df.columns)

# Select features (all columns except 'Grade_encoded')
features = encoded_df.drop(columns=['Grade_encoded'])

# Select label ('Grade_encoded')
label = encoded_df['Grade_encoded']

# Split the data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(features, label, test_size=0.3, random_state=42)

# Normalize the features
scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)

# Train an SVM model
svm_model = SVC(kernel='linear')
svm_model.fit(X_train, y_train)

# Predict on the test set
y_pred = svm_model.predict(X_test)

# Evaluate the model
accuracy = accuracy_score(y_test, y_pred)
print(f"SVM Accuracy: {accuracy:.2f}")

print("\nSVM Classification Report:")
print(classification_report(y_test, y_pred))

# Train a Random Forest model
rf_model = RandomForestClassifier(n_estimators=50, random_state=42)
rf_model.fit(X_train, y_train)

# Predict on the test set
y_pred_rf = rf_model.predict(X_test)

# Evaluate the Random Forest model
accuracy_rf = accuracy_score(y_test, y_pred_rf)
print(f"Random Forest Accuracy: {accuracy_rf:.2f}")
print("\nRandom Forest Classification Report:")
print(classification_report(y_test, y_pred_rf))

# Plot feature importance
plt.figure(figsize=(10, 8))
sns.barplot(x=rf_model.feature_importances_, y=features.columns)
plt.title('Feature Importance - Random Forest')
plt.xlabel('Importance')
plt.ylabel('Feature')
plt.show()

def plot_learning_curves(estimator, X, y, title='Learning Curves'):
    train_sizes, train_scores, test_scores = learning_curve(estimator, X, y, cv=5, n_jobs=-1,
                                                            train_sizes=np.linspace(0.1, 1.0, 10), scoring='accuracy')

    train_scores_mean = np.mean(train_scores, axis=1)
    train_scores_std = np.std(train_scores, axis=1)
    test_scores_mean = np.mean(test_scores, axis=1)
    test_scores_std = np.std(test_scores, axis=1)

    plt.figure(figsize=(12, 6))
    plt.title(title)
    plt.xlabel("Training examples")
    plt.ylabel("Score")
    plt.grid()

    plt.fill_between(train_sizes, train_scores_mean - train_scores_std, train_scores_mean + train_scores_std, alpha=0.1,
                     color="r")
    plt.fill_between(train_sizes, test_scores_mean - test_scores_std, test_scores_mean + test_scores_std, alpha=0.1,
                     color="g")
    plt.plot(train_sizes, train_scores_mean, 'o-', color="r", label="Training score")
    plt.plot(train_sizes, test_scores_mean, 'o-', color="g", label="Cross-validation score")

    plt.legend(loc="best")
    return plt


# Plot learning curves for the Random Forest model
plot_learning_curves(rf_model, X_train, y_train, title='Random Forest Learning Curves')
plt.show()

# Define the parameter grid
param_grid = {
    'n_estimators': [50, 100, 150],
    'max_depth': [None, 10, 20, 30],
    'min_samples_split': [2, 5, 10],
    'min_samples_leaf': [1, 2, 4],
    'max_features': ['auto', 'sqrt', 'log2']
}

# Initialize a RandomForestClassifier
rf = RandomForestClassifier(random_state=42)

# Initialize the GridSearchCV object
grid_search = GridSearchCV(estimator=rf, param_grid=param_grid, cv=5, scoring='accuracy', verbose=2, n_jobs=-1)

# Fit the grid search to the data
grid_search.fit(X_train, y_train)

# Best parameters and best score
print("Best parameters found: ", grid_search.best_params_)
print("Best accuracy found: ", grid_search.best_score_)

# Retrieve the best model
best_rf_model = grid_search.best_estimator_

# Predict on the test set using the best model
y_pred_rf = best_rf_model.predict(X_test)

# Evaluate the Random Forest model
accuracy_rf = accuracy_score(y_test, y_pred_rf)
print(f"Optimized Random Forest Accuracy: {accuracy_rf:.2f}")
print("\nOptimized Random Forest Classification Report:")
print(classification_report(y_test, y_pred_rf))

# Train an XGBoost model
xgb_model = XGBClassifier(use_label_encoder=False, eval_metric='logloss')
xgb_model.fit(X_train, y_train)

# Predict on the test set
y_pred = xgb_model.predict(X_test)

# Evaluate the model
accuracy = accuracy_score(y_test, y_pred)
print(f"XGBOOST Accuracy: {accuracy:.2f}")

print("\nXGBOOST Classification Report:")
print(classification_report(y_test, y_pred))

# Cross-validation
cross_val_scores = cross_val_score(xgb_model, X_train, y_train, cv=5)
print(f"Cross-validation scores: {cross_val_scores}")
print(f"Mean cross-validation score: {cross_val_scores.mean():.2f}")

# Plot feature importance
plt.figure(figsize=(10, 8))
sns.barplot(x=xgb_model.feature_importances_, y=features.columns)
plt.title('Feature Importance')
plt.xlabel('Importance')
plt.ylabel('Feature')
plt.show()

# Adjusting model parameters
xgb_model = XGBClassifier(
    use_label_encoder=False,
    eval_metric='logloss',
    max_depth=6,
    colsample_bytree=0.8,
    gamma=0.1,
    min_child_weight=1
)
xgb_model.fit(X_train, y_train)

# Evaluating the model
y_pred = xgb_model.predict(X_test)
accuracy = accuracy_score(y_test, y_pred)
print(f"Adjusted XGBOOST Accuracy: {accuracy:.2f}")
print("\nAdjusted XGBOOST Classification Report:")
print(classification_report(y_test, y_pred))

# Cross-validation to check stability
cross_val_scores = cross_val_score(xgb_model, X_train, y_train, cv=5)
print(f"Adjusted Cross-validation scores: {cross_val_scores}")
print(f"Mean Adjusted cross-validation score: {cross_val_scores.mean():.2f}")

# Feature importance
plt.figure(figsize=(10, 8))
sns.barplot(x=xgb_model.feature_importances_, y=features.columns)
plt.title('Adjusted Feature Importance')
plt.xlabel('Importance')
plt.ylabel('Feature')
plt.show()