Classifying pediatric IBD stool samples (work in progress)¶
This notebook is a recoding of the analysis used in the PLoSONE paper: Non-Invasive Mapping of the Gastrointestinal Microbiota Identifies Children with Inflammatory Bowel Disease using python, sklearn and pandas.
We thought that the SLiME package, as it was packaged for the publication of the paper became outdated and should probably not be available anymore. This notebook replaces it, replicating the analysis executed on the paper with more up-to-date tools and with a few extra figures. I hope this can be the starting point for others trying to follow the same approach and improve upon it.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
%matplotlib inline
from scipy import interp
from sklearn.ensemble import RandomForestClassifier
from sklearn.multiclass import OneVsRestClassifier
from sklearn.cross_validation import cross_val_score, StratifiedKFold, train_test_split, KFold
from sklearn.metrics import roc_curve, roc_auc_score, auc
from sklearn.preprocessing import LabelEncoder, label_binarize
print('pandas ' + pd.__version__)
print('mpl ' +mpl.__version__)
print('numpy' + np.__version__)
print('seaborn ' + sns.__version__)
Loading data¶
The data came from two rounds of 16S sequencing of previously collected stool samples. Here we will use the OTU tables directly, which were created by using the RDP classifier and were subsequently normalized (details in the paper's methods).
Sequencing was performed at the Broad Institute. The first round of sequencing was dubbed CHIMP (Children Hospital IBD Pediatric), while the second round of sequencing -- performed following the request of an anonymous peer reviewer -- was termed blind validation. Its purpose was to further validate the algorithm trained on the CHIMP dataset, as one of the reviewers did not think it was sufficient to use a "leave 20% out" approach on the CHIMP dataset to demonstrate robust prediction.
While these rounds were used as training and test set in the last figure of the paper respectively, it is more useful at this stage to join the two data sets and split them in different random combinations later.
#get the CHIMP training data
X_chimp = pd.read_csv('data/chimp/chimp.Qsorted.rdpout.xtab.norm', delimiter="\t", index_col=0)
y_chimp = pd.read_csv('data/chimp/sampledata.training.chimp.csv', index_col=0)
#just make sure the labels are the same
X_chimp.sort_index(inplace=True)
y_chimp.sort_index(inplace=True)
assert (X_chimp.index == y_chimp.index).all()
## do the same for the blind validation test data
X_blind = pd.read_csv('data/chimp/blind.sorted.rdpout.xtab.norm',
delimiter="\t", index_col=0)
y_blind = pd.read_csv('data/chimp/sampledata.validation.blind.csv',
index_col=0)
X_blind.sort_index(inplace=True)
y_blind.sort_index(inplace=True)
assert (X_blind.index == y_blind.index).all()
#concatenate using pandas
X = pd.concat([X_chimp, X_blind], keys=['chimp','blind'])
X.head()
X.fillna(value=0,inplace=True) #replace NAs with zeroes
y_dx = pd.concat([y_chimp.dx, y_blind.dx], keys=['chimp','blind'])
print(y_dx.head())
print(y_dx.value_counts())
#convert the training and testing labels to numerical values
le = LabelEncoder()
le.fit(y_dx)
y = le.transform(y_dx)
# just for reference, the columns of the binarized label read respectively:
le.inverse_transform([0,1,2])
Exploratory data analysis¶
# map CD,NM and UC to healthy and IBD classes
y_ibd = y_dx.map({'CD':'ibd','UC':'ibd','NM':'healthy'})
y_ibd.value_counts()
# split OTU by phylogenetic level
expl = X.copy()
expl.columns = expl.columns.str.partition('_')
expl.head()
# how many healthy people have a given OTU vs how many IBD patients have the same OTU
# we use the mean, dividing by the total number of patients in each group to get percentages
otupresence = expl.astype(bool).groupby(y_ibd).mean()
otup = otupresence.stack(level=[0,2]).reset_index()
plt.figure(figsize=(20,4))
f = sns.stripplot(data=otup,x='level_2',y='_',hue='dx')
f.axes.xaxis.grid(True)
f.axes.yaxis.grid(False)
ibdotu = otupresence.stack(level=[0,2]).unstack('dx').reset_index()
ibdotu.columns.values
ibdotu.columns=['level','name','healthy','ibd']
ibdotu_abd = ibdotu[(ibdotu.healthy < 0.1) & (ibdotu.ibd > 0.05) & (ibdotu.level != 'domain')]
print('We found %s OTUs (out of %s ~= %1.1f %%) that are \n \
present in IBD (> 5%% of patients) and underrepresented (< 10 %%) in healthy controls\n' % \
(ibdotu_abd.shape[0],ibdotu.shape[0],(ibdotu_abd.shape[0]/ibdotu.shape[0])*100))
print(ibdotu_abd)
# let's make it easier to see.
#multiply ibd percentages by -1
otup.loc[otup['dx'] == 'ibd',"_"] *= -1
plt.figure(figsize=(20,4))
f = sns.stripplot(data=otup,x='level_2',y='_',hue='dx')
f.axes.xaxis.grid(True)
f.axes.yaxis.grid(False)
diff = otup.groupby(['level_1','level_2']).sum().reset_index()
diff.head()
g = sns.FacetGrid(diff[diff.level_1 != 'domain'], col="level_1",sharey=False,size=20, aspect=.3)
g.map(sns.barplot,"_","level_2")
g.set_xlabels('ibd <== enriched in ==> healthy')
sns.set_style("whitegrid")
g = sns.FacetGrid(diff[diff.level_1 != 'domain'], col="level_1",sharey=False,size=20, aspect=.3)
g.map(sns.stripplot,"_","level_2")
g.set_xlabels('ibd <== enriched in ==> healthy')
for ax in g.axes.flat:
# Make the grid horizontal instead of vertical
# ax.xaxis.grid(False)
ax.yaxis.grid(True)
f = sns.factorplot(x='_',y='level_2',col='dx',data=otup[otup.level_1 == 'genus'],kind='bar',size=20,aspect=.25)
f.set_axis_labels('percentage of patients which has OTU',"")
plt.savefig('OTUs.png')