<jupyter_start><jupyter_text>Data source:<NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, Lau LMS, <NAME>, <NAME>, <NAME>, Grimmond SM, <NAME>, <NAME>, <NAME>, Scolyer RA, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>. **Telomere sequence content can be used to determine ALT activity in tumours.**_Nucleic Acids Res._ 2018 Jun 1;46(10):4903-4918.  STEP 1: Preprocessing  Load data<jupyter_code>import pandas as pd

# Load data
data = pd.read_csv("data/telomere_ALT/telomere.csv", sep='\t')
data.head(10)    # Show first ten samples<jupyter_output><empty_output><jupyter_text>Split data into training and test sets<jupyter_code>from sklearn import model_selection

num_row = data.shape[0]     # number of samples in the dataset
num_col = data.shape[1]     # number of features in the dataset (plus the label column)

X = data.iloc[:, 0: num_col-1]    # feature columns
y = data['TMM']                   # label column

X_train, X_test, y_train, y_test = \
    model_selection.train_test_split(X, y, 
                                     test_size=0.2,      # reserve 20 percent data for testing
                                     stratify=y,         # stratified sampling
                                     random_state=0)
X_train.shape
X_test.shape
y_train.value_counts()
y_test.value_counts()<jupyter_output><empty_output><jupyter_text>STEP 2: Learning  Train a Random Forest classifier<jupyter_code>import numpy as np
from sklearn import ensemble

rf = ensemble.RandomForestClassifier(
    n_estimators = 10,           # 10 random trees in the forest
    criterion = 'entropy',       # use entropy as the measure of uncertainty
    max_depth = 3,               # maximum depth of each tree is 3
    min_samples_split = 5,       # generate a split only when there are at least 5 samples at current node
    class_weight = 'balanced',   # class weight is inversely proportional to class frequencies
    random_state = 0)

k = 3   # number of folds

# split data into k folds
kfold = model_selection.StratifiedKFold(
    n_splits = k,
    shuffle = True,
    random_state = 0)

cv = list(kfold.split(X_train, y_train))    # indices of samples in each fold
for j, (train_idx, val_idx) in enumerate(cv):
    print('validation fold %d:' % (j+1))
    print(val_idx)<jupyter_output>validation fold 1:
[  3   6   7  10  12  13  20  24  31  34  36  37  39  42  45  56  57  59
  62  63  69  70  71  75  77  81  83  88  90  91  95  97 100 104 105 106
 107 108 111 113 116 117 127]
validation fold 2:
[  0   2   4   8   9  11  14  16  17  18  22  25  26  29  33  38  40  43
  46  48  51  53  54  55  65  66  68  72  73  76  78  80  82  84  93  98
 103 112 114 119 120 121 123]
validation fold 3:
[  1   5  15  19  21  23  27  28  30  32  35  41  44  47  49  50  52  58
  60  61  64  67  74  79  85  86  87  89  92  94  96  99 101 102 109 110
 115 118 122 124 125 126]<jupyter_text>Compute classifier accuracy<jupyter_code>accuracy = []

# Compute classifier's accuracy on each fold
for j, (train_idx, val_idx) in enumerate(cv):
    rf.fit(X_train.iloc[train_idx], y_train.iloc[train_idx])
    accuracy_per_fold = rf.score(X_train.iloc[val_idx], y_train.iloc[val_idx])
    accuracy.append(accuracy_per_fold)
    print('Fold %d accuracy: %.3f' % (j+1, accuracy_per_fold))

print('Average accuracy: %.3f' % np.mean(accuracy))<jupyter_output>Fold 1 accuracy: 0.860
Fold 2 accuracy: 0.860
Fold 3 accuracy: 0.857
Average accuracy: 0.859<jupyter_text>Plot receiver operating characteristic (ROC) curve<jupyter_code>from matplotlib import pyplot as plt
from sklearn import metrics

fpr = dict()       # false positive rate
tpr = dict()       # true positive rate
auroc = dict()     # area under the ROC curve (AUROC)

# Compute an ROC curve (with respect to the positive class) for each fold
# Compute AUROC for each fold
for j, (train_idx, val_idx) in enumerate(cv):
    rf.fit(X_train.iloc[train_idx], y_train.iloc[train_idx])
    y_prob = rf.predict_proba(X_train.iloc[val_idx])
    fpr[j], tpr[j], _ = metrics.roc_curve(y_train.iloc[val_idx], y_prob[:, 0], pos_label='+')
    auroc[j] = metrics.auc(fpr[j], tpr[j])

# Compute an average ROC curve for all folds
fpr['avg'] = np.unique(np.concatenate([fpr[j] for j in range(k)]))
tpr['avg'] = np.zeros_like(fpr['avg'])

for j in range(k):
    tpr['avg'] += np.interp(fpr['avg'], fpr[j], tpr[j])

tpr['avg'] /= k

# Compute AUROC of the average ROC curve
auroc['avg'] = metrics.auc(fpr['avg'], tpr['avg'])
# Plot the ROC curves
for j in range(k):
    plt.plot(fpr[j], tpr[j], label='fold %d (area = %.3f)' % (j+1, auroc[j]))

plt.plot(fpr['avg'], tpr['avg'], 'k--', label='mean ROC (area = %.3f)' % auroc['avg'], lw=2)
plt.plot([0,1], [0,1], linestyle='--', color=(0.6, 0.6, 0.6))     # ROC curve for a random classifier
plt.plot([0,0,1], [0,1,1], linestyle='--', color=(0.6, 0.6, 0.6)) # ROC curve for a perfect classifier

plt.xlim([-0.05, 1.05])
plt.ylim([-0.05, 1.05])
plt.xlabel('FPR')
plt.ylabel('TPR')
plt.legend(loc='lower right')

plt.show()<jupyter_output><empty_output><jupyter_text>Plot precision-recall (PR) curve<jupyter_code>precision_plus = dict()
recall_plus = dict()
auprc_plus = dict()

precision_minus = dict()
recall_minus = dict()
auprc_minus = dict()

# Compute PR curves (with respect to the positive and negative classes) for each fold
# Compute AUPRC for each fold
for j, (train_idx, val_idx) in enumerate(cv):
    rf.fit(X_train.iloc[train_idx], y_train.iloc[train_idx])
    y_prob = rf.predict_proba(X_train.iloc[val_idx])
    precision_plus[j], recall_plus[j], _ = metrics.precision_recall_curve(y_train.iloc[val_idx], y_prob[:, 0], pos_label='+')
    precision_minus[j], recall_minus[j], _ = metrics.precision_recall_curve(y_train.iloc[val_idx], y_prob[:, 1], pos_label='-')
    auprc_plus[j] = metrics.auc(recall_plus[j], precision_plus[j])
    auprc_minus[j] = metrics.auc(recall_minus[j], precision_minus[j])

# Compute average PR curves for all folds
recall_plus['avg'] = np.unique(np.concatenate([recall_plus[j] for j in range(k)]))
recall_minus['avg'] = np.unique(np.concatenate([recall_minus[j] for j in range(k)]))
precision_plus['avg'] = np.zeros_like(recall_plus['avg'])
precision_minus['avg'] = np.zeros_like(recall_minus['avg'])

for j in range(k):
    precision_plus['avg'] += np.interp(recall_plus['avg'], np.flip(recall_plus[j], 0), np.flip(precision_plus[j], 0))
    precision_minus['avg'] += np.interp(recall_minus['avg'], np.flip(recall_minus[j], 0), np.flip(precision_minus[j], 0))

precision_plus['avg'] /= k
precision_minus['avg'] /= k

# Compute AUPRCs for the average curves
auprc_plus['avg'] = metrics.auc(recall_plus['avg'], precision_plus['avg'])
auprc_minus['avg'] = metrics.auc(recall_minus['avg'], precision_minus['avg'])<jupyter_output><empty_output><jupyter_text>PR curve w.r.t the positive class<jupyter_code># Plot the PR curves for the positive class
for j in range(k):
    plt.plot(recall_plus[j], precision_plus[j], label='fold %d (area = %.3f)' % (j+1, auprc_plus[j]))

plt.plot(recall_plus['avg'], precision_plus['avg'], 'k--', label='mean PR (area = %.3f)' % auprc_plus['avg'], lw=2)
plt.plot([0,1,1], [1,1,0], linestyle='--', color=(0.6, 0.6, 0.6))   # PR curve for a perfect classifier

plt.xlim([-0.05, 1.05])
plt.ylim([-0.05, 1.05])
plt.xlabel('recall')
plt.ylabel('precision')
plt.legend(loc='lower left')

plt.show()<jupyter_output><empty_output><jupyter_text>PR curve w.r.t the negative class<jupyter_code># Plot the PR curves for the negative class
for j in range(k):
    plt.plot(recall_minus[j], precision_minus[j], label='fold %d (area = %.3f)' % (j+1, auprc_minus[j]))

plt.plot(recall_minus['avg'], precision_minus['avg'], 'k--', label='mean PR (area = %.3f)' % auprc_minus['avg'], lw=2)
plt.plot([0,1,1], [1,1,0], linestyle='--', color=(0.6, 0.6, 0.6))   # PR curve for a perfect classifier

plt.xlim([-0.05, 1.05])
plt.ylim([-0.05, 1.05])
plt.xlabel('recall')
plt.ylabel('precision')
plt.legend(loc='lower left')

plt.show()<jupyter_output><empty_output><jupyter_text>Plot confusion matrix<jupyter_code># Compute a 2x2 confusion matrix
confusion_matrix = np.zeros([2, 2])

for j, (train_idx, val_idx) in enumerate(cv):
    rf.fit(X_train.iloc[train_idx], y_train.iloc[train_idx])
    y = rf.predict(X_train.iloc[val_idx])
    confusion_matrix += metrics.confusion_matrix(y_train.iloc[val_idx], y)

# Average over all folds
confusion_matrix /= k
confusion_matrix = np.around(confusion_matrix, decimals=3)
# Plot the confusion matrix
plt.matshow(confusion_matrix, interpolation='nearest', cmap=plt.cm.Blues, alpha=0.5)
plt.xticks(np.arange(2), ('+', '-'))
plt.yticks(np.arange(2), ('+', '-'))

for i in range(2):
    for j in range(2):
        plt.text(j, i, confusion_matrix[i, j], ha='center', va='center')

plt.xlabel('Predicted label')
plt.ylabel('True label')

plt.show()<jupyter_output><empty_output><jupyter_text>STEP 3: Evaluation<jupyter_code># Refit the classifier using the full training set
rf = rf.fit(X_train, y_train)<jupyter_output><empty_output><jupyter_text>Compute classifier accuracy on the test set<jupyter_code># Compute classifier's accuracy on the test set
test_accuracy = rf.score(X_test, y_test)
print('Test accuracy: %.3f' % test_accuracy)<jupyter_output>Test accuracy: 0.848<jupyter_text>Generate plots for the test set<jupyter_code># Predict labels for the test set
y = rf.predict(X_test)             # Solid prediction
y_prob = rf.predict_proba(X_test)  # Soft prediction

# Compute the ROC and PR curves for the test set
fpr, tpr, _ = metrics.roc_curve(y_test, y_prob[:, 0], pos_label='+')
precision_plus, recall_plus, _ = metrics.precision_recall_curve(y_test, y_prob[:, 0], pos_label='+')
precision_minus, recall_minus, _ = metrics.precision_recall_curve(y_test, y_prob[:, 1], pos_label='-')

# Compute the AUROC and AUPRCs for the test set
auroc = metrics.auc(fpr, tpr)
auprc_plus = metrics.auc(recall_plus, precision_plus)
auprc_minus = metrics.auc(recall_minus, precision_minus)

# Compute the confusion matrix for the test set
cm = metrics.confusion_matrix(y_test, y)
# Plot the ROC curve for the test set
plt.plot(fpr, tpr, label='test (area = %.3f)' % auroc)
plt.plot([0,1], [0,1], linestyle='--', color=(0.6, 0.6, 0.6))
plt.plot([0,0,1], [0,1,1], linestyle='--', color=(0.6, 0.6, 0.6))

plt.xlim([-0.05, 1.05])
plt.ylim([-0.05, 1.05])
plt.xlabel('FPR')
plt.ylabel('TPR')
plt.legend(loc='lower right')
plt.title('ROC curve')

plt.show()
# Plot the PR curve for the test set
plt.plot(recall_plus, precision_plus, label='positive class (area = %.3f)' % auprc_plus)
plt.plot(recall_minus, precision_minus, label='negative class (area = %.3f)' % auprc_minus)

plt.plot([0,1,1], [1,1,0], linestyle='--', color=(0.6, 0.6, 0.6))

plt.xlim([-0.05, 1.05])
plt.ylim([-0.05, 1.05])
plt.xlabel('recall')
plt.ylabel('precision')
plt.legend(loc='lower left')
plt.title('Precision-recall curve')

plt.show()
# Plot the confusion matrix for the test set
plt.matshow(cm, interpolation='nearest', cmap=plt.cm.Blues, alpha=0.5)
plt.xticks(np.arange(2), ('+', '-'))
plt.yticks(np.arange(2), ('+', '-'))

for i in range(2):
    for j in range(2):
        plt.text(j, i, cm[i, j], ha='center', va='center')

plt.xlabel('Predicted label')
plt.ylabel('True label')

plt.show()<jupyter_output><empty_output><jupyter_text>STEP 4: Prediction<jupyter_code># Enter unlabeled data for prediction
newdata = pd.read_csv("data/telomere_ALT/telomere_new.csv", sep='\t')
newdata.head(6)
# predict labels
y_new = rf.predict(newdata)

y_new<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15483758599669076, 'fraction_numerical': 0.046242271183488635, 'mean_word_length': 4.5, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 67, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '92722', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 4812, 'n_tokens_neox': 4364, 'n_words': 1093, 'path': 'illustration.ipynb', 'repo_name': 'carpentries-incubator/ml4bio-workshop', 'size': '11525'}