Visualising metagenomic data collected over time
I’m going to play with some metagenomics data and rplotly. The animated rplotly, created in RMarkdown, is run on this static site.
I’m drawing on a wide undergraduate and coursework masters statistics background, a number of years of data analysis experience, and my Graduate Cert in Molecular Biology. My “go to” for how to analyse ’omics data is Gordon Smyth of the Walter and Eliza Hall Institute. Gordon is author of several R Bioconductor packages.
I looked for some shotgun metagenomics data to use for an example, and spent a couple of hours looking at some data that is supplied with an R package called “metamicrobiomeR” (not one of Professor Smyth’s). I will try to succinctly describe the data and what I see as some of the problems with it, before I do some exploratory analysis with r plotly and use the Bioconductor function limma for analysis.
The data are relative abundances in the metagenomes of fifty babies faecal samples. Most babies had their faecal metagenomes analysed 6 times, roughly once every month from birth, though one baby was only analysed 3 times and one baby had 10 microbiome analyses. There are a total of 322 observations. All 50 babies came from one location, Bangladesh, but the data has two levels for location; the other is Myanmar. I’m curious where the data for the Myanmar babies is, and knowing some data has been excluded means that understanding the nature of the dataset, as presented, is more complicated. Some demographic variables are provided, including the variable I have decided is of interest, the breastfed status of the baby. I’ve seen breastfeeding data before and really, the statistician needs to be able to ask the people who collected the data what was actually meant by certain responses. I haven’t seen any background information about this data, but I haven’t looked for any yet as this is a quick analysis to try out r plotly and limma.
There is suggestion in the data that archaea were found in some babies, as there are columns for archaea. However, all values are 0. The only data is for bacteria. This is unfortunate as I would have been interested to see some data that included archaea, viruses and microbial eukaryotes as well as bacteria. The data is in the form of relative abundance values, like shotgun metagenomic data would be.
The variable bf, breastfed status, has 3 levels: ExclusiveBF, Non-exclusiveBF and No_BF. I decided, since I am not fully informed about the nature of the data, and what I am aiming to do on this web page is show my thinking processes as to how I would analyse shotgun metagenomic data, to just use data from the 8 babies that were exclusively breastfed on all occasions their faecal sample was assessed, and the 6 babies who were always non-exclusively breastfed. Five of the 6 non-exclusively breastfed babies were twins or triplets; knowing this I would not normally go ahead with the simple sort of analysis I am planning, as it seems birth multiplicity is associated with breastfed status, but as I have said, I just want to have a look at the nature of metagenomic abundance data analysis. If there is a difference in the relative abundances within microbiomes, we will have to qualify it by saying it could be due to singleton versus multiple birth status rather than breastfed status.
The idea is to put together this analysis and web page over a couple of days, but if I planned to do anything more comprehensive, I would consider how to fit a repeated measures model to this data, as the babies were measured repeatedly over time. As it is with only a couple of days to work on this, as part of the exploratory analysis, I have decided to analyse the data from 14 babies in their sixth month, giving the microbiome time to develop. At each measurement, whether the baby received antibiotics was noted, as was whether they had received oral rehydrating fluid. I created variables to indicate if the babies had received antibiotics or oral rehydrating fluid at any time before age 6 months.
On to the actual relative abundances data. The github page where this package was developed suggests that the relative abundance data was trimmed. That is, where the relative abundance was small, it was set to 0. Further, if less than 5% of observations had the microbe, it was set to 0. There are no viruses or archaea in the data and I’m concerned they may have been trimmed out inappropriately.
Adding over columns of data shows that the relative abundances add to 1. If trimming was done, the relative abundances have been recalculated.
The data is presented firstly as a breakdown into abundance of phylum within kingdom (22 columns, the only kingdom being bacteria), then further into classes (42 columns), orders (77 columns), families (173 columns), and genus (427 columns). 285 of the 740 columns (38%) have only 0s in them. I’m not sure if this is from data being zeroed, whether if the Myanmar data was included there would be some non-zero values, or if perhaps some sort of phylogenetic tree is used to decide what might be in the microbiome.
All 8 of the exclusively breastfed babies received antibiotics at least once before they were 6 months old but only 2 of the 6 babies that were not exclusively breastfed had antibiotics before 6 months. It is possible that relative abundances in the microbiota would be different in babies given antibiotics.
Five of the 8 exclusively breastfed babies were given oral rehydrating solution before they were 6 months old. This means these babies were not actually exclusively breastfed. None of the “not exclusively breastfed” babies had oral rehydrating fluid.
Despite feeling miserable about the dodginess of the data, I will press on! First a look at the abundances. I have created this plot with r plotly. Hovering with your cursor over a bar will give you a “tooltip”, with the baby’s demographics. Clicking on the legend will toggle highlighting those sections of the stacked bar chart. Note that only four phyla were present in these babies’ faecal samples to any degree but some babies have small relative abundances of the other four phyla listed in the legend.
I ran a few analyses to see if there is a difference in relative abundance within microbiomes between exclusively breastfed babies and babies who were not exclusively breastfed. In the plot above, the babies who were exclusively breastfed were the first baby starting from the left, and the 3rd to the 9th babies. (Normally I’d rearrange the order of the babies so that all the breastfed ones were in a group, but with the data being the mess it is, I wondered if perhaps I should be asking if relative abundance is different in babies who had antibiotics, or in babies who were part of multiple births, or in babies who had received oral hydration. In the end I left the babies in the order of their index numbers.)
Once the dataset was reduced to 14 babies, I removed all rows of the expression set that were 0 for all 14 babies. Note that although no data can be seen for the phyla thermi, fusobacteria, tenericutes and unassigned.other in the plot, some babies had very small relative abundances of these. If you switch off all the other legend entries you may be able to see that baby 3 has a very small amount of unassigned.other, for example.
After removing rows with all zeros, I added 0.01 and took logs base 2. This is as per the suggestion of Gordon Smyth in a Bioconductor forum discussion from 3 months ago. I then did a limma-trend analysis, also as per Gordon’s suggestion.
Breastfed status did not make a difference to relative abundance at any level (from genus to phylum).
Neither did gender, antibiotics, rehydration, or multiple birth status. It wasn’t the intention to look at any variables except breastfed status, so my data dredging is unscientific. Normally my analysis would be informed by experts who had decided prior to the study what the research questions and hypotheses are.
Since I am playing around here, I decided to make a dummy variable that coincides with the babies that have high relative abundance of the proteobacteria phylum. This is the 4th, 5th, 7th and 13th babies. I just wanted to see that this dummy variable would come up as significant.
dummy <- factor(c(0,0,0,1,1,0,1,0,0,0,0,0,1,0))
design <- model.matrix(~dummy)
fit <- lmFit(exampleSet,design)
fit <- eBayes(fit,trend=TRUE)
results <- decideTests(fit)
summary(results) (Intercept) dummy1
Down 8 0
NotSig 0 7
Up 0 1topTable(fit, coef="dummy1", n=3) logFC AveExpr t
k__bacteria.p__proteobacteria 3.7125784838 -4.2489941 8.601509
k__bacteria.p__.thermi. 0.0008306631 -6.6436189 1.479136
k__bacteria.p__actinobacteria -0.3927525584 -0.6459723 -1.128654
P.Value adj.P.Val B
k__bacteria.p__proteobacteria 1.237392e-06 9.899134e-06 5.752242
k__bacteria.p__.thermi. 1.636173e-01 6.544692e-01 -6.030665
k__bacteria.p__actinobacteria 2.800341e-01 7.467576e-01 -6.461327The limma-trend analysis is performed by the two “fit” lines in the code snippet.
The “1” in the cell that is the intersection of the dummy1 column and the Up row in the table that is the output of the summary(results) command, indicates that the relative abundance of one phyla is increased in the babies in the “dummy” group. The topTable command then shows that it is indeed the proteobacteria phylum, with small adjusted p-value.
Normally I’d make a conclusion addressing the research questions. A proper analysis would have used much more of the data. I have only considered 14 of the 322 observations.
If I was actually analysing this data, I’d find out why the thermi phylum is written as .thermi., and check this was consistent in the data. I’d get a nice colour palette, although it’s funny that the main colour in the plot has defaulted to breastfed baby poo mustard!
Moving on now, to animating plotly.
Animating ggplotly, and r plotly generally, only works on scatter points, including geom_segment. That’s why in this animation, rather than a stacked bar chart, I have thick line segments connecting the cumulative relative abundances. This chart with thick line segments look the same as a stacked bar chart on my laptop, but on some devices the thickness of the lines will look wrong.
Hovering tooltips don’t work with segments - I didn’t want tooltips in this animation anyway.
I made the animation option “transition = 0” because with any sort of positive transition time, some of the line segments come flying in from the right! I would have liked to see all the abundances smoothly change as the animation plays, but I’ll have to think of how to get that. The problem arises when the abundances of some phyla are 0 at one time and not zero at the next. Also when “transition \(\gt\) 0”, segments that were very small but not zero at one time point, that become bigger at the next time point, begin the transition as thin, black lines and evolve to the thick, coloured state, which looks weird. I have in mind how I might do this by fixing the start of every line segment at y=0 and drawing the cumulative abundance y values over the top of one another in decreasing order.
Animating is unofficial in ggplotly; I got the error message “Ignoring unknown aesthetics: frame, ids” but it worked anyway.
g <- ggplot(pFPhyla,
aes(x=Baby, xend=Baby, y=abundstart, yend=abundend, colour=Bacteria)) +
geom_segment(stat="identity", size=4, aes(frame= Month, ids=Bacteria)) +
theme(axis.text.x = element_text(angle=30,hjust=1,vjust=1)) +
ylab("Abundance") +
scale_color_futurama()
ggplotly(g,tooltip="") %>% animation_opts(800, mode="afterall",transition=0)I quite like this animation. If nothing else, you can see that the babies had a large relative abundance of proteobacteria to begin with but later there is a predominance of actinobacteria.