<jupyter_start><jupyter_text>Motif Kernel / SVM TutorialThis Tutorial shows how you can combine the motif kernel of the *strkernel* package with a Support Vector Machine (SVM) to predict the cell population based on the motif content of a read sequence. There are two FASTA files each filled with sequences from two different cell popultions. If you are not familar with the FASTA format, [here](http://genetics.bwh.harvard.edu/pph/FASTA.html) is a short explantion. *fibroblast.fa* are sequences obtained from fibroblast while *stemcells.fa* contains sequences from stemcells. The goal of this tutorial is to show that we can use prior knowledge to construct motifs and use those motifs to classify new sequences into the two cell populations. The accuracy of the prediction is in this case secondary as long as we get a result that shows that we can use the motif content of sequences to compare their similarity and therefore classify them.First, you will need some packages for preprocessing, classification and to plot the results.<jupyter_code># preprocessing
import numpy as np
from Bio.Seq import Seq
import Bio.SeqIO as sio

# SVM
from sklearn.svm import SVC
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report

# ROC and precision-recall curve
from sklearn.metrics import roc_curve, auc, precision_recall_curve, average_precision_score

# plotting
import matplotlib.pyplot as plt  # plotting

# motif kernel
from strkernel.motifkernel import motifKernel<jupyter_output><empty_output><jupyter_text>PreprocessingIn order to use the data with the SVM provided by *sklearn* we need to read in the FASTA files with *Biopython* and add some labels. In this case we will label the stemcells as positive (1) and the fibroblast as negative (0).<jupyter_code># load the data
# stemcells
pos_data = [seq.seq for seq in sio.parse('notebook_data/stemcells.fa', 'fasta')]
# fibroblasts
neg_data = [seq.seq for seq in sio.parse('notebook_data/fibroblast.fa', 'fasta')]

pos_labels = np.ones(len(pos_data), dtype=int)
neg_labels = np.zeros(len(neg_data), dtype=int)

y = np.concatenate((pos_labels, neg_labels), axis=0)<jupyter_output><empty_output><jupyter_text>MotifsNow we have to decide which collection of motifs we will be using to construct the motif kernel. There are specific transcription factors binding sites which tend to be enriched in one of the groups. Naturally not every sequence contains the binding site for these transcription factors but we should be able to correctly classify those that do. In stemcells, we expect the **oct4** binding site to be enriched because this protein has been shown to be heavily involved in the pluripotency of cells. In fibroblast lung tissue on the other hand, we expect the **mafk** binding site to be enriched.The [jasper](http://jaspar.genereg.net/) database allows to check for the binding motifs of the transcription factors. If we search for **oct4** and **mafk** we get the following binding sites: **oct4**![](notebook_data/oct4.png) **mafk**![](notebook_data/mafk.png)We can now add the most prominent motifs of these bindings sites to our motif collection. For this tutorial we chose the motif *TCAGCA* of the **mafk** binding site and the motifs *ATGCAA* and *TTGT* of the **oct4** binding site. This decision is based on which motifs are the most prominent. Since the reads originate from both strands we also have to include the reverse complement of the motifs. In normal use cases the number of motifs is usally a lot higher but for this tutorial the six motifs will be enough.<jupyter_code>motif_collection = ["TCAGCA","TGCTGA","ATGCAA","TTGCAT","TTGT","ACAA"]

#create the motif kernel
motif_kernel = motifKernel(motif_collection)

#use the motif kernel to compute the motif content matrix for all sequences
motif_matrix = motif_kernel.compute_matrix(pos_data + neg_data)<jupyter_output><empty_output><jupyter_text>ClassificationWe can now split the matrix into test and training data. The training data can then be used to train the SVM. We will keep 30% of the data as test data to evaluate the model.<jupyter_code>#split the data into test and training set
X_train, X_test, y_train, y_test = train_test_split(motif_matrix, y, test_size=0.3, random_state=42, stratify=y)
#train the classifier
clf = SVC()
clf.fit(X_train, y_train)<jupyter_output><empty_output><jupyter_text>ResultsThe only thing left to do is to analyze the trained model. *sklearn* provides a function which can be used to produce a classifcation report. This report shows us the precision, recall and f1-score when we apply the model to our test data. We will also plot the ROC and PRC but first we will have to define a couple of wrapper functions.<jupyter_code>def plot_roc_curve(y_test, y_score):
    '''Plots a roc curve including a baseline'''
    fpr, tpr, thresholds = roc_curve(y_test, y_score)
    roc_auc = auc(fpr, tpr)

    plt.figure()
    lw = 2
    plt.plot(fpr, tpr, color='darkorange',
             lw=lw, label='ROC curve (area = %0.2f)' % roc_auc)
    plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
    plt.xlim([0.0, 1.0])
    plt.ylim([0.0, 1.05])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    plt.title('Receiver operating curve')
    plt.legend(loc="lower right")
    plt.show()


def plot_prec_recall_curve(y_test, y_scores):
    '''Plots a precision-recall curve including a baseline'''
    precision, recall, thresholds = precision_recall_curve(y_test, y_scores)
    average_precision = average_precision_score(y_test, y_scores)
    baseline = np.bincount(y_test)[1] / sum(np.bincount(y_test))
    plt.figure()
    plt.step(recall, precision, color='b', alpha=0.2,
             where='post')
    plt.fill_between(recall, precision, step='post', alpha=0.2,
                     color='b')
    plt.axhline(y=baseline, linewidth=2, color='navy', linestyle='--')
    plt.xlabel('Recall')
    plt.ylabel('Precision')
    plt.ylim([0.0, 1.05])
    plt.xlim([0.0, 1.0])
    plt.title('2-class Precision-Recall curve: AP={0:0.2f}'.format(
              average_precision))
    plt.show()


def reportclassfication(clf, X_test, y_test):
    '''Reports classification results with the given model and testdata'''

    print("Detailed classification report:")
    print()
    y_true, y_pred = y_test, clf.predict(X_test)
    print(classification_report(y_true, y_pred))
    print()
    

# get classfication results with the tuned parameters
reportclassfication(clf, X_test, y_test)

y_scores = clf.decision_function(X_test)
# plot ROC
plot_roc_curve(y_test, y_scores)

# plot PRC
plot_prec_recall_curve(y_test, y_scores)<jupyter_output>Detailed classification report:

             precision    recall  f1-score   support

          0       0.77      0.65      0.70      1200
          1       0.70      0.80      0.75      1199

avg / total       0.73      0.73      0.72      2399
{'fraction_non_alphanumeric': 0.08540103866128101, 'fraction_numerical': 0.013271783035199077, 'mean_word_length': 4.374418604651162, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '19248', 'license': 'MIT', 'n_tokens_mistral': 2223, 'n_tokens_neox': 2062, 'n_words': 804, 'path': 'docs/source/.ipynb_checkpoints/motiftutorial-checkpoint.ipynb', 'repo_name': 'jakobhertzberg/string-kernel', 'size': '6932'}