Microbiome experiment step-by-step analysis¶

Import the calour module¶

In [1]:

import calour as ca

/Users/amnon/miniconda3/envs/calour/lib/python3.5/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters

(optional) Set the level of feedback messages from calour¶

can use:

• 1 for debug (lots of feedback on each command)
• 11 for info (useful information from some commands)
• 21 for warning (just warning messages)

The Calour default is warning (21)

In [2]:

ca.set_log_level(11)

Also enable interactive plots inside the jupyter notebook¶

In [3]:

%matplotlib notebook

Get rid of the features (bacteria) with small amount of reads¶

We throw away all features with total reads (over all samples) < 10 (after each sample was normalized to 10k reads/sample). So a bacteria present (with 1 read) in 10 samples will be kept, as well as a bacteria present in only one sample, with 10 reads in this sample. Note alternatively we could filter based on mean reads/sample or fraction of samples where the feature is present. Each method filters away slightly different bacteria. See filtering notebook for details on the filtering functions.

In [6]:

dat=dat.filter_abundance(10)

2018-03-04 12:40:30 INFO After filtering, 1100 remaining

Cluster (reorder) the features so similarly behaving bacteria are close to each other¶

Features are clustered (hierarchical clustering) based on euaclidian distance between features (over all samples) following normalizing each feature to mean 0 std 1. For more details and examples, see sorting notebook or cluster_features documentation

• Note that if we have a lot of features, clustering is slow, so it is recommended to first filter away the non-interesting features.

In [7]:

datc=dat.cluster_features()

2018-03-04 12:40:30 INFO After filtering, 1100 remaining

Sort the samples according to physical functioning and Disease state¶

Note that order within each group of similar value is maintained. We first sort by physical functioning, then sort by the disease state. So within each disease state, samples will still be sorted by physical functioning.

In [8]:

datc=datc.sort_samples('Physical_functioning')
datc=datc.sort_samples('Subject')

Adding a field to the top bar¶

Now let’s add the values of the “Sex” field into the xbar on top First we’ll also sort by sex, so values will be continuous (note we then sort by the disease state to get the two groups separated).

In [10]:

datc=datc.sort_samples('Sex')
datc=datc.sort_samples('Subject')

In [11]:

datc.plot(sample_field='Subject', gui='jupyter',barx_fields=['Sex'])

Out[11]:

<calour.heatmap.plotgui_jupyter.PlotGUI_Jupyter at 0x1a19dacc50>

Differential abundance testing¶

Let’s look for bacteria separating sick from healthy We ask it to find all bacteria significantly different between samples with ‘Control’ and ‘Patient’ in the ‘Subject’ field.

By default calour uses the mean of the ranks of each feature (over all samples), with dsFDR multiple hypothesis correction.

For more information, see notebook and function doc

In [12]:

dd=datc.diff_abundance(field='Subject',val1='Control',val2='Patient', random_seed=2018)

2018-03-04 12:40:36 INFO 87 samples with both values
2018-03-04 12:40:36 INFO After filtering, 1100 remaining
2018-03-04 12:40:36 INFO 39 samples with value 1 (['Control'])
2018-03-04 12:40:37 INFO method meandiff. number of higher in ['Control'] : 44. number of higher in ['Patient'] : 19. total 63

Plotting the differentially abundant features¶

Let’s plot to see the behavior of these bacteria. The output of diff_abundance is an Experiment with only the significant bacteria, which are sorted by the effect size. On the bottom is the bacteria with the largest effect size (higher in Control compared to Patient).

In [13]:

dd.plot(sample_field='Subject', gui='jupyter')

Out[13]:

<calour.heatmap.plotgui_jupyter.PlotGUI_Jupyter at 0x1a0ffeca20>

dbBact term enrichment¶

We can ask what is special in the bacteria significanly higher in the Control vs. the Patient group and vice versa.

We supply the parameter ignore_exp=[12] to ignore annotations regarding this experiment (expid=12) since it is already in the dbBact database.

• Note since we need to get the per-feature annotations from dbBact, we need a live internet connection to run this command.

In [14]:

ax, enriched=dd.plot_diff_abundance_enrichment(term_type='combined',ignore_exp=[12])

2018-03-04 12:40:41 INFO Getting dbBact annotations for 63 sequences, please wait...
2018-03-04 12:40:45 INFO Got 2383 annotations
2018-03-04 12:40:45 INFO Added annotation data to experiment. Total 631 annotations, 63 terms
2018-03-04 12:40:45 INFO removed 171 terms

The enriched terms are in a calour experiment class (terms are features, bacteria are samples), so we can see the list of enriched terms with the p-value (pval) and effect size (odif)

In [15]:

enriched.feature_metadata

Out[15]:

In [ ]: