{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nPractical Deep Learning for Genomic Prediction\n==============================================\n\nCredit: A Grigis\n\nBased on:\n\n- https://github.com/miguelperezenciso/DLpipeline\n\nLoad the data\n-------------\n\nLoad some data.\nYou may need to change the 'datasetdir' parameter.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import os\nfrom pynet.datasets import DataManager, fetch_genomic_pred\nfrom pynet.utils import setup_logging\n\nsetup_logging(level=\"info\")\n\ndata = fetch_genomic_pred(\n datasetdir=\"/tmp/genomic_pred\")\nmanager = DataManager(\n input_path=data.input_path,\n labels=[\"env0\"],\n metadata_path=data.metadata_path,\n number_of_folds=2,\n batch_size=5,\n test_size=0.2,\n continuous_labels=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Basic inspection\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nfrom sklearn.decomposition import PCA\nimport matplotlib.pyplot as plt\n\ntrain_dataset = manager[\"train\"][0]\nX_train = train_dataset.inputs[train_dataset.indices]\ny_train = train_dataset.labels[train_dataset.indices]\ntest_dataset = manager[\"test\"]\nX_test = test_dataset.inputs[test_dataset.indices]\ny_test = test_dataset.labels[test_dataset.indices]\nprint(X_train.shape, y_train.shape)\nprint(X_test.shape, y_test.shape)\nprint(\" min max mean sd\")\nprint(\"Train:\", y_train.min(), y_train.max(), y_train.mean(),\n np.sqrt(y_train.var()))\nprint(\"Test:\", y_test.min(), y_test.max(), y_test.mean(),\n np.sqrt(y_test.var()))\nplt.figure()\nplt.title(\"Train / test data\")\nplt.hist(y_train, label=\"Train\")\nplt.hist(y_test, label=\"Test\")\nplt.legend(loc=\"best\")\nX = np.concatenate((X_train, X_test))\npca = PCA(n_components=2)\np = pca.fit(X).fit_transform(X)\nNtrain = X_train.shape[0]\nplt.figure()\nplt.title(\"PCA decomposition\")\nplt.scatter(p[0:Ntrain, 0], p[0:Ntrain, 1], label=\"Train\")\nplt.scatter(p[Ntrain:, 0], p[Ntrain:, 1], label=\"Test\", color=\"orange\")\nplt.legend(loc=\"best\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "SNP preselection according to a simple GWAS: select N_best most\nassociated SNPs or select by min_P_value.\nOptional: not used after.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from scipy import stats\n\npvals = []\nfor idx in range(X_train.shape[1]):\n b, intercept, r_value, p_value, std_err = stats.linregress(\n X_train[:, idx], y_train)\n pvals.append(-np.log10(p_value))\npvals = np.array(pvals)\nplt.figure()\nplt.ylabel(\"-log10 P-value\")\nplt.xlabel(\"SNP\")\nplt.plot(pvals, marker=\"o\")\nN_best = 100\nsnp_list = pvals.argsort()[-N_best:].squeeze().tolist()\nmin_P_value = 2 # P = 0.01\nprint(np.nonzero(pvals > min_P_value))\nsnp_list = np.nonzero(pvals > min_P_value)[0].squeeze().tolist()\nX_train_filter = X_train[:, snp_list]\nX_test_filter = X_test[:, snp_list]\nprint(X_train_filter.shape, y_train.shape)\nprint(X_test_filter.shape, y_test.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Apply standard penalized methods (lasso using scikit-learn).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from sklearn import linear_model\nfrom sklearn.metrics import mean_squared_error\n\nlasso = linear_model.Lasso(alpha=0.01)\nlasso.fit(X_train, y_train)\ny_hat = lasso.predict(X_test)\nmse = mean_squared_error(y_test, y_hat)\nprint(\"MSE in prediction =\", mse)\ncorr = np.corrcoef(y_test, y_hat)[0, 1]\nprint(\"Corr obs vs pred =\", corr)\nplt.figure()\nplt.title(\"Lasso: Observed vs Predicted Y\")\nplt.ylabel(\"Predicted\")\nplt.xlabel(\"Observed\")\nplt.scatter(y_test, y_hat, marker=\"o\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Implements a standard fully connected network (MLP) for a quantitative\ntarget.\nUse Mean Squared Error as loss, ie, quantitative variable, regression.\nWe apply a kernel regularization on the first linear layer to punish the\nweights which are very large causing the network to overfit, after applying\nthis regularization the weights will become smaller.\nWe also apply an activity regularization on the first layer that tries to\nmake the output smaller so as to remove overfitting.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import collections\nimport torch\nimport torch.nn as nn\nfrom pynet.utils import get_named_layers\nfrom pynet.interfaces import DeepLearningInterface\n\n\nclass TwoLayersMLP(nn.Module):\n \"\"\" Simple two hidden layers percetron.\n \"\"\"\n def __init__(self, data_size, nb_neurons, nb_classes, drop_rate=0.2):\n \"\"\" Initialize the instance.\n\n Parameters\n ----------\n data_size: int\n the number of elements in the data.\n nb_neurons: 2-uplet with int\n the number of neurons of the hidden layers.\n nb_classes: int\n the number of classes.\n drop_rate: float, default 0.2\n the dropout rate.\n \"\"\"\n super(TwoLayersMLP, self).__init__()\n self.nb_classes = nb_classes\n self.layers = nn.Sequential(collections.OrderedDict([\n (\"linear1\", nn.Linear(data_size, nb_neurons[0])),\n (\"activation1\", nn.ReLU()),\n (\"linear2\", nn.Linear(nb_neurons[0], nb_neurons[1])),\n (\"activation2\", nn.Softplus()),\n (\"drop1\", nn.Dropout(drop_rate)),\n (\"linear3\", nn.Linear(nb_neurons[1], nb_classes))\n ]))\n\n def forward(self, x):\n layer1_out = self.layers[0](x)\n x = self.layers[1:](layer1_out)\n if self.nb_classes == 1:\n x = x.view(x.size(0))\n return x, {\"layer1\": layer1_out}\n\n\ndef linear1_l2_kernel_regularizer(signal):\n lambda2 = 0.01\n model = signal.object.model\n all_linear2_params = torch.cat([\n x.view(-1) for x in model.layers[0].parameters()])\n l2_regularization = lambda2 * torch.norm(all_linear2_params, 2)\n return l2_regularization\n\n\ndef linear1_l1_activity_regularizer(signal):\n lambda1 = 0.01\n layer1_out = model = signal.layer_outputs[\"layer1\"]\n l1_regularization = lambda1 * torch.norm(layer1_out, 1)\n return l1_regularization\n\n\nnb_snps = X_train.shape[1]\nmodel = TwoLayersMLP(nb_snps, nb_neurons=[64, 32], nb_classes=1)\nprint(model)\ncl = DeepLearningInterface(\n optimizer_name=\"SGD\",\n learning_rate=5e-4,\n loss_name=\"MSELoss\",\n metrics=[\"pearson_correlation\"],\n model=model)\ncl.add_observer(\"regularizer\", linear1_l2_kernel_regularizer)\ncl.add_observer(\"regularizer\", linear1_l1_activity_regularizer)\ntest_history, train_history = cl.training(\n manager=manager,\n nb_epochs=(100 if \"CI_MODE\" not in os.environ else 10),\n checkpointdir=\"/tmp/genomic_pred\",\n fold_index=0,\n with_validation=True)\ny_hat, X, y_true, loss, values = cl.testing(\n manager=manager,\n with_logit=False,\n predict=False)\nprint(y_hat.shape, y_true.shape)\nprint(y_hat)\nprint(y_true)\nprint(\"MSE in prediction =\", loss)\ncorr = np.corrcoef(y_true, y_hat)[0, 1]\nprint(\"Corr obs vs pred =\", corr)\nplt.figure()\nplt.title(\"MLP: Observed vs Predicted Y\")\nplt.ylabel(\"Predicted\")\nplt.xlabel(\"Observed\")\nplt.scatter(y_test, y_hat, marker=\"o\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Implements the same probblem but with a Convolutional Neural Network (CNN)\nfor a quantitative target.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class MyNet(torch.nn.Module):\n def __init__(self):\n super(MyNet, self).__init__()\n self.conv1 = torch.nn.Conv1d(1, 32, kernel_size=3, stride=3, padding=1)\n self.maxpool = torch.nn.MaxPool1d(kernel_size=2)\n self.linear = nn.Sequential(collections.OrderedDict([\n (\"linear1\", nn.Linear(32 * 213, 64)),\n (\"activation1\", nn.ReLU()),\n (\"linear2\", nn.Linear(64, 32)),\n (\"activation2\", nn.Softplus()),\n (\"linear3\", nn.Linear(32, 1))\n ]))\n\n def forward(self, x):\n x = x.view(x.shape[0], 1, x.shape[1])\n x = self.maxpool(self.conv1(x))\n x = x.view(-1, 32 * 213)\n x = self.linear(x)\n x = x.view(x.size(0))\n return x\n\n\nmodel = MyNet()\nprint(model)\ncl = DeepLearningInterface(\n optimizer_name=\"SGD\",\n learning_rate=5e-4,\n loss_name=\"MSELoss\",\n metrics=[\"pearson_correlation\"],\n model=model)\ntest_history, train_history = cl.training(\n manager=manager,\n nb_epochs=(50 if \"CI_MODE\" not in os.environ else 10),\n checkpointdir=\"/tmp/genomic_pred\",\n fold_index=0,\n with_validation=True)\ny_hat, X, y_true, loss, values = cl.testing(\n manager=manager,\n with_logit=False,\n predict=False)\nprint(y_hat.shape, y_true.shape)\nprint(y_hat)\nprint(y_true)\nprint(\"MSE in prediction =\", loss)\ncorr = np.corrcoef(y_true, y_hat)[0, 1]\nprint(\"Corr obs vs pred =\", corr)\nplt.figure()\nplt.title(\"MLP: Observed vs Predicted Y\")\nplt.ylabel(\"Predicted\")\nplt.xlabel(\"Observed\")\nplt.scatter(y_test, y_hat, marker=\"o\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Implements the same fully connected network (MLP) for a quantitative\ntarget but in the context of multiclass target.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data = fetch_genomic_pred(\n datasetdir=\"/tmp/genomic_pred\",\n to_categorical=True)\nmanager = DataManager(\n input_path=data.input_path,\n labels=[\"env0_cat0\", \"env0_cat1\", \"env0_cat2\"],\n stratify_label=\"env0\",\n projection_labels={\"env0\": [0, 1, 2]},\n metadata_path=data.metadata_path,\n number_of_folds=2,\n batch_size=5,\n test_size=0.2)\ntrain_dataset = manager[\"train\"][0]\nX_train = train_dataset.inputs[train_dataset.indices]\ny_train = train_dataset.labels[train_dataset.indices]\ntest_dataset = manager[\"test\"]\nX_test = test_dataset.inputs[test_dataset.indices]\ny_test = test_dataset.labels[test_dataset.indices]\nprint(X_train.shape, y_train.shape)\nprint(X_test.shape, y_test.shape)\nnb_snps = X_train.shape[1]\ny_train = manager[\"train\"][0].labels[train_dataset.indices]\nprint(y_train.shape)\nmodel = TwoLayersMLP(nb_snps, nb_neurons=[64, 32], nb_classes=3)\nprint(model)\n\n\ndef my_loss(x, y):\n \"\"\" nn.CrossEntropyLoss expects a torch.LongTensor containing the class\n indices without the channel dimension.\n \"\"\"\n device = y.get_device()\n y = torch.argmax(y, dim=1).type(torch.LongTensor)\n if device != -1:\n y = y.to(device)\n criterion = nn.CrossEntropyLoss()\n return criterion(x, y)\n\n\ncl = DeepLearningInterface(\n optimizer_name=\"Adam\",\n learning_rate=5e-4,\n loss=my_loss,\n model=model)\ntest_history, train_history = cl.training(\n manager=manager,\n nb_epochs=(100 if \"CI_MODE\" not in os.environ else 10),\n checkpointdir=\"/tmp/genomic_pred\",\n fold_index=0,\n with_validation=True)\ny_hat, X, y_true, loss, values = cl.testing(\n manager=manager,\n with_logit=True,\n predict=False)\nprint(y_hat.shape, y_true.shape)\nprint(y_hat)\nprint(y_true)\nprint(\"MSE in prediction =\", loss)\nheat = np.zeros([3, 3])\nfor i in range(3):\n klass = np.nonzero(y_true[:, i] > 0)\n for j in range(3):\n heat[i, j] = np.mean(y_hat[klass, j])\nprint(\"Probabilities matrix\", heat)\nplt.figure()\nplot = plt.imshow(heat, cmap=\"Blues\")\nplt.ylabel(\"Predicted class\")\nplt.xlabel(\"Observed class\")\n\nif \"CI_MODE\" not in os.environ:\n plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.12" } }, "nbformat": 4, "nbformat_minor": 0 }