Chemotherapy is one of the most used treatment against cancer. It uses chemical substances (chemotherapeutic agents) which kill cells that divide too quickly. These chemical substances are often diluted in a particular solution and packaged in bags, diffusers, or syringes, before being administered. Wrong medication (wrong chemotherapeutic agent or wrong concentration) can have major impacts for patients. To prevent wrong medication, some recent French regulations impose the verification of anti-cancer drugs before their administration. The goal is to check that they contain the good chemotherapeutic agent with the good dosage.
Raman spectroscopy could be used to make this check, since, theoretically, i) each molecule has a specific spectral fingerprint by which the molecule can be identified; and ii) the Raman intensity increases with the concentration of the molecule. The main advantage of spectroscopy above other methods (for example, liquid chromatography) is that it is non-destructive and non-invasive (measures are made without opening the drug containers). However, this method is rarely used in hospital environment because of the complexity of the spectral signals to analyze. Automating the analysis of these spectral signals could significantly help. Eventually, a complete analytical system (from measuring Raman spectra to identifying the chemotherapeutic agent and its concentration) could be designed, which would be easy to use and would prevent wrong medication.
In this context, the goal of this project is to develop prediction models able to identify and quantify chemotherapeutic agents from their Raman spectra.
The Lip(Sys)² laboratory measured Raman spectra of 4 types of chemotherapeutic agents (called molecule) in 3 different packages (called vial), diluted in 9 different solutions (called solute gammes), and having different concentrations. A total of 360 spectra were measured for each agent, except for one (348 spectra).
Part of these data are saved in the file train.csv as follows (n_samples being the number of samples):
n_samples,)1, n_samples)1, n_samples)n_samples, 1)n_samples, 1866)To sum up, there are too objectives:
We can read spectra_0.csv, spectra_1.csv, spectra_2.csv, spectra_3.csv which will be used to train a model of machine learning in order to predict either a concentration of molecule or a molecule type.
In [1]:
import pandas
from numpy import *
s0=pandas.read_csv('data/spectra_0.csv')
c0=s0['concentration']
m0=s0['molecule']
s0=s0['spectra']
x=[]
for spec in s0:
x.append(fromstring(spec[1:-1], sep=','))
s0=pandas.DataFrame(x)
s1=pandas.read_csv('data/spectra_1.csv')
c1=s1['concentration']
m1=s1['molecule']
s1=s1['spectra']
x=[]
for spec in s1:
x.append(fromstring(spec[1:-1], sep=','))
s1=pandas.DataFrame(x)
s2=pandas.read_csv('data/spectra_2.csv')
c2=s2['concentration']
m2=s2['molecule']
s2=s2['spectra']
x=[]
for spec in s2:
x.append(fromstring(spec[1:-1], sep=','))
s2=pandas.DataFrame(x)
s3=pandas.read_csv('data/spectra_3.csv')
c3=s3['concentration']
m3=s3['molecule']
s3=s3['spectra']
x=[]
for spec in s3:
x.append(fromstring(spec[1:-1], sep=','))
s3=pandas.DataFrame(x)
f=pandas.read_csv('data/freq.csv')
stot=pandas.concat([s0, s1, s2, s3])
c=pandas.concat([c0, c1, c2, c3])
m=pandas.concat([m0, m1, m2, m3])
In [2]:
f.head()
Out[2]:
In [3]:
c.head()
Out[3]:
In [4]:
m.head()
Out[4]:
In [5]:
stot.head()
Out[5]:
In [6]:
from matplotlib import pyplot
%matplotlib inline
Fig, ax = pyplot.subplots()
ax.plot(f, stot.T)
ax.set_xlabel('Frequency')
ax.set_ylabel('Intensity')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
ax.spines['left'].set_position(('outward',10))
ax.spines['bottom'].set_position(('outward',10))
In [7]:
import numpy
numpy.unique(m)
Out[7]:
In [8]:
fig, ax = pyplot.subplots()
for m_type in numpy.unique(m):
i = numpy.nonzero(m==m_type)[0]
ax.plot(f,numpy.mean(stot.iloc[i],axis=0),label=m_type)
ax.fill_between(numpy.ravel(f),numpy.mean(stot.iloc[i],axis=0)+numpy.std(stot.iloc[i],axis=0),numpy.mean(stot.iloc[i],axis=0)-numpy.std(stot.iloc[i],axis=0),alpha=0.2)
ax.set_xlabel('Frequency')
ax.set_ylabel('Intensity')
ax.legend()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
ax.spines['left'].set_position(('outward',10))
ax.spines['bottom'].set_position(('outward',10))
In [9]:
numpy.unique(c)
Out[9]:
In [10]:
fig, ax = pyplot.subplots()
for c_type in numpy.unique(c):
i = numpy.nonzero(c == c_type)[0]
ax.plot(f, numpy.mean(stot.iloc[i], axis=0), label=c_type)
ax.fill_between(numpy.ravel(f), numpy.mean(stot.iloc[i], axis=0) + numpy.std(stot.iloc[i], axis=0), numpy.mean(stot.iloc[i], axis=0) - numpy.std(stot.iloc[i], axis=0), alpha=0.2)
ax.set_xlabel('Frequency')
ax.set_ylabel('Intensity')
ax.legend(loc='center left', bbox_to_anchor = (1,0.5))
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
ax.spines['left'].set_position(('outward', 10))
ax.spines['bottom'].set_position(('outward', 10))
We are going to use spectra_4.csv for testing a model learned on the previous data
In [11]:
s4 = pandas.read_csv('data/spectra_4.csv')
f = pandas.read_csv('data/freq.csv')
c4 = s4['concentration']
m4 = s4['molecule']
s4 = s4['spectra']
In [12]:
x = []
for spec in s4:
x.append(numpy.fromstring(spec[1:-1], sep=','))
s4 = pandas.DataFrame(x)
s4.head()
Out[12]:
In [13]:
fig, ax = pyplot.subplots()
ax.plot(f, s4.T)
ax.set_xlabel('Frequency')
ax.set_ylabel('Intensity')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
ax.spines['left'].set_position(('outward', 10))
ax.spines['bottom'].set_position(('outward', 10))
In [14]:
fig, ax = pyplot.subplots()
for m_type in numpy.unique(m4):
i = numpy.nonzero(m4 == m_type)[0]
ax.plot(f, numpy.mean(s4.iloc[i], axis=0), label=m_type)
ax.fill_between(numpy.ravel(f), numpy.mean(s4.iloc[i], axis=0) + numpy.std(s4.iloc[i], axis=0), numpy.mean(s4.iloc[i], axis=0) - numpy.std(s4.iloc[i], axis=0), alpha=0.2)
ax.set_xlabel('Frequency')
ax.set_ylabel('Intensity')
ax.legend()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
ax.spines['left'].set_position(('outward', 10))
ax.spines['bottom'].set_position(('outward', 10))
In [15]:
fig, ax = pyplot.subplots()
for c_type in numpy.unique(c4):
i = numpy.nonzero(c4 == c_type)[0]
ax.plot(f, numpy.mean(s4.iloc[i], axis=0), label=c_type)
ax.fill_between(numpy.ravel(f), numpy.mean(s4.iloc[i], axis=0) + numpy.std(s4.iloc[i], axis=0), numpy.mean(s4.iloc[i], axis=0) - numpy.std(s4.iloc[i], axis=0), alpha=0.2)
ax.set_xlabel('Frequency')
ax.set_ylabel('Intensity')
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
ax.spines['left'].set_position(('outward', 10))
ax.spines['bottom'].set_position(('outward', 10))
In [16]:
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.ensemble import RandomForestClassifier
from sklearn.pipeline import make_pipeline
p = make_pipeline(StandardScaler(), PCA(n_components=100, random_state=0), RandomForestClassifier(random_state=0))
p.fit(stot, m)
print('Accuracy score: {0:.2f}'.format(p.score(s4, m4)))
In [17]:
from sklearn.metrics import confusion_matrix
import itertools
pr = p.predict(s4)
cm = confusion_matrix(m4, pr)
print(cm)
pyplot.imshow(cm, interpolation='nearest', cmap='bwr')
pyplot.title('Confusion matrix')
pyplot.colorbar()
tick_marks = numpy.arange(len(numpy.unique(m4)))
pyplot.xticks(tick_marks, numpy.unique(m4), rotation=45)
pyplot.yticks(tick_marks, numpy.unique(m4))
fmt = 'd'
thresh = cm.max() / 2.
for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
pyplot.text(j, i, format(cm[i, j], fmt),horizontalalignment="center",color="white" if cm[i, j] > thresh else "black")
pyplot.tight_layout()
pyplot.ylabel('True label')
pyplot.xlabel('Predicted label')
Out[17]:
In [18]:
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.svm import LinearSVC
from sklearn.pipeline import make_pipeline
p = make_pipeline(StandardScaler(), PCA(n_components=100, random_state=0), LinearSVC(random_state=0))
p.fit(stot, m)
print('Accuracy score: {0:.2f}'.format(p.score(s4, m4)))
In [19]:
from sklearn.metrics import confusion_matrix
import itertools
pr = p.predict(s4)
cm = confusion_matrix(m4, pr)
print(cm)
pyplot.imshow(cm, interpolation='nearest', cmap='bwr')
pyplot.title('Confusion matrix')
pyplot.colorbar()
tick_marks = numpy.arange(len(numpy.unique(m4)))
pyplot.xticks(tick_marks, numpy.unique(m4), rotation=45)
pyplot.yticks(tick_marks, numpy.unique(m4))
fmt = 'd'
thresh = cm.max() / 2.
for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
pyplot.text(j, i, format(cm[i, j], fmt),horizontalalignment="center",color="white" if cm[i, j] > thresh else "black")
pyplot.tight_layout()
pyplot.ylabel('True label')
pyplot.xlabel('Predicted label')
Out[19]:
In [20]:
med_stot = numpy.median(stot, axis=0)
p_25_stot = numpy.percentile(stot,25,axis=0)
p_75_stot = numpy.percentile(stot,75,axis=0)
stot = (stot-med_stot)/(p_75_stot-p_25_stot)
s4 = (s4-med_stot)/(p_75_stot-p_25_stot)
In [21]:
from sklearn.decomposition import PCA
from sklearn.linear_model import RidgeCV
from sklearn.pipeline import make_pipeline
p = make_pipeline(PCA(n_components=100), RidgeCV())
p.fit(stot, c)
print('R2 score: {0:.2f}'.format(p.score(s4, c4)))
In [22]:
from sklearn.metrics import r2_score, median_absolute_error
y_pred = p.predict(s4)
fig, ax = pyplot.subplots()
ax.scatter(c4, y_pred)
ax.plot([0, 25000], [0, 25000], '--k')
ax.set_ylabel('Target predicted')
ax.set_xlabel('True Target')
ax.text(1000, 20000, r'$R^2$=%.2f, MAE=%.2f' % (r2_score(c4, y_pred), median_absolute_error(c4, y_pred)))
ax.set_xlim([0, 25000])
ax.set_ylim([0, 25000])
Out[22]:
In [23]:
from sklearn.decomposition import PCA
from sklearn.ensemble import RandomForestRegressor
from sklearn.pipeline import make_pipeline
p = make_pipeline(PCA(n_components=100, random_state=0), RandomForestRegressor(random_state=0))
p.fit(stot, c)
print('R2 score: {0:.2f}'.format(p.score(s4, c4)))
In [24]:
from sklearn.metrics import r2_score, median_absolute_error
y_pred = p.predict(s4)
fig, ax = pyplot.subplots()
ax.scatter(c4, y_pred)
ax.plot([0, 25000], [0, 25000], '--k')
ax.set_ylabel('Target predicted')
ax.set_xlabel('True Target')
ax.text(1000, 20000, r'$R^2$=%.2f, MAE=%.2f' % (r2_score(c4, y_pred), median_absolute_error(c4, y_pred)))
ax.set_xlim([0, 25000])
ax.set_ylim([0, 25000])
Out[24]: