MAREDAT_consistency_check/.ipynb_checkpoints/consistency_with_tara_oceans-checkpoint.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load dependencies\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Consistency check between the MAREDAT data and *Tara* Oceans data\n",
    "Our estimates of the global biomass of several marine taxa are based on data from the MAREDAT database. As stated in the specific sections relying on data from the database, there are many sources of uncertainty associated with the estimates stemming from the MAREDAT data. Many of those sources of uncertainty are hard to quantify, and no uncertainty estimate is provided in the literature on estimates based on the MAREDAT database.\n",
    "\n",
    "Thus, we perform consistency checks for the MAREDAT data against independent sources of data, to increases our confidence in our estimates and to provide a measure of the uncertainty associated with our estimates.\n",
    "\n",
    "In this document we conduct a comparison between the estimates of biomass based on the MAREDAT database and data from [de Vargas et al.](http://dx.doi.org/10.1126/science.1261605). the data in de Vargas et al. is based on 18S rDNA sequencing of different populations of plankton collected by the *Tara* Oceans campaign. \n",
    "\n",
    "The dataset in de Vargas et al. divides the plankton community based on size ranges (pico-nano-, nano-, micro- and meso-plankton). de Vargas et al. provides only number of reads for each taxon. The fraction of reads that a taxon has out of the total number of reads can be used as a proxy for the biomass fraction of the taxon, but not as a proxy of its absolute biomass. Relying on 18S rDNA sequence abundance as a proxy for biomass is not a well established practice, and has its own biases, but we chose to use it for the sake of comparing it to independent approaches such as the MAREDAT database. Each plankton size fraction sampled in the study was sequenced to approximately the same sequencing depth (≈120 million reads). This means that the 18S read data can provide a possible proxy for the biomass fraction of a certain taxon within a size fraction, but not across size fractions.\n",
    "\n",
    "We focus on comparing the MAREDAT and de Vargas et al. data in two case studies: the biomass of diatoms and the total biomass of nanoplakton and microplankton.\n",
    "\n",
    "## Diatoms\n",
    "We begin by describing how to compare the biomass estimates of diatoms based on de Vargas et al. and the MAREDAT dataset. Our aim is to calculate the relative fraction of diatoms out of the total biomass of organisms in the same size range as diatoms. As the data in those two datasets is structured differently, we first need to make corrections to the data so a valid comparison will be available.\n",
    "\n",
    "### MAREDAT\n",
    "In de Vargas et al., diatoms appear mainly in the nanoplankton (5-20 µm in diameter) and microplankton (20-180 µm) size fractions. In order to make a comparison to the MAREDAT database we need to find the corresponding groups in the MAREDAT database. The corresponding groups in the MAREDAT database are the microzooplankton and the diatom groups (zooplankton between 5 and 200 µm in diameter). As calculated in the marine protists sections, our estimates for the respective biomass of microzooplankton and diatoms are ≈0.6 Gt C and ≈0.3 Gt C. Thus, according to the MAREDAT data diatoms account for about 30% of the total biomass of plankton in the 5-200 µm size fraction.\n",
    "\n",
    "### de Vargas et al.\n",
    "We use data on the total number of reads of different taxa in each size fraction. The data originates from de Vargas et al. from Database W6 in the companion website, as well as from Figure 3 in the main text. Here is a sample of the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Group</th>\n",
       "      <th>Total # of reads</th>\n",
       "      <th>Piconano fraction</th>\n",
       "      <th>Nano fraction</th>\n",
       "      <th>Micro fraction</th>\n",
       "      <th>Meso fraction</th>\n",
       "      <th>Rhizaria</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Metazoa</td>\n",
       "      <td>246962646</td>\n",
       "      <td>0.050202</td>\n",
       "      <td>0.242105</td>\n",
       "      <td>0.370850</td>\n",
       "      <td>0.336032</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Bacillariophyta</td>\n",
       "      <td>14586250</td>\n",
       "      <td>0.123622</td>\n",
       "      <td>0.540945</td>\n",
       "      <td>0.292913</td>\n",
       "      <td>0.039370</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Collodaria</td>\n",
       "      <td>97085064</td>\n",
       "      <td>0.052061</td>\n",
       "      <td>0.238612</td>\n",
       "      <td>0.160521</td>\n",
       "      <td>0.548807</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Nassellaria &amp; Eucyrtidium</td>\n",
       "      <td>656961</td>\n",
       "      <td>0.429501</td>\n",
       "      <td>0.119306</td>\n",
       "      <td>0.234273</td>\n",
       "      <td>0.216920</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Spumellaria</td>\n",
       "      <td>5199377</td>\n",
       "      <td>0.791757</td>\n",
       "      <td>0.054230</td>\n",
       "      <td>0.071584</td>\n",
       "      <td>0.080260</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                       Group  Total # of reads  Piconano fraction  \\\n",
       "0                    Metazoa         246962646           0.050202   \n",
       "1            Bacillariophyta          14586250           0.123622   \n",
       "2                 Collodaria          97085064           0.052061   \n",
       "3  Nassellaria & Eucyrtidium            656961           0.429501   \n",
       "4                Spumellaria           5199377           0.791757   \n",
       "\n",
       "   Nano fraction  Micro fraction  Meso fraction  Rhizaria  \n",
       "0       0.242105        0.370850       0.336032     False  \n",
       "1       0.540945        0.292913       0.039370     False  \n",
       "2       0.238612        0.160521       0.548807      True  \n",
       "3       0.119306        0.234273       0.216920      True  \n",
       "4       0.054230        0.071584       0.080260      True  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load data on the total number of reads of each taxon from de Vargas et al.\n",
    "data = pd.read_excel('tara_oceans_data.xlsx','de Vargas W6',skiprows=1)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also use data on the total number of reads from each size fraction from Figure 2 in de Vargas et al.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Piconano reads</th>\n",
       "      <th>Nano reads</th>\n",
       "      <th>Micro reads</th>\n",
       "      <th>Meso reads</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>114000000</td>\n",
       "      <td>135000000</td>\n",
       "      <td>121000000</td>\n",
       "      <td>135000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Piconano reads  Nano reads  Micro reads  Meso reads\n",
       "0       114000000   135000000    121000000   135000000"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Load data on the total number of reads in each size fraction\n",
    "tot_reads = pd.read_excel('tara_oceans_data.xlsx','Total number of reads', skiprows=1)\n",
    "tot_reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In de Vargas et al., diatoms appear mainly in the nanoplankton (5-20 µm in diameter) and microplankton (20-180 µm) size fractions. However, the microzooplankton biomass estimates in the MAREDAT database do not include copepods, which were moved to the mesozooplankton group. Fragile protists such as Rhizaria, are probably also undersampled in the MAREDAT database. Therefore, to correct for these effects such that we could compare the MAREDAT and Tara Oceans datasets, we remove metazoa (dominated by arthropods) and Rhizaria reads from the relevant size fractions in the Tara Oceans dataset (nano and microplakton):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The fraction of diatoms out of the total number of reads in nanoplankton and microplankton\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Nano reads</th>\n",
       "      <th>Micro reads</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.156385</td>\n",
       "      <td>0.337967</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Nano reads  Micro reads\n",
       "0    0.156385     0.337967"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate the total number of reads for the Nano and Micro fractions\n",
    "read_data = pd.DataFrame()\n",
    "read_data['Nano reads'] = data['Total # of reads']*data['Nano fraction']\n",
    "read_data['Micro reads'] = data['Total # of reads']*data['Micro fraction']\n",
    "\n",
    "# Subtract the total sum of rhizaria and metazoa from the total number of reads\n",
    "corrected_total_reads = tot_reads[['Nano reads','Micro reads']] - read_data[data['Group'] == 'Metazoa'] - read_data[data['Rhizaria'] == True].sum()\n",
    "\n",
    "print('The fraction of diatoms out of the total number of reads in nanoplankton and microplankton')\n",
    "read_data.loc[1]/corrected_total_reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After correcting for those biases, the biomass fraction of diatoms in microplankton in the Tara Oceans dataset is between 16%-33%, which fits well with the estimate from the MAREDAT database of ≈30%.\n",
    "\n",
    "## Nanoplankton and Microplankton biomass\n",
    "In this section we generate an independent estimate of the total biomass of nanoplankton and microplankton, based on several data sources. We begin with the independently measured biomass of Rhizaria. The independent measurement using microscopy by  [Biard et al.](http://dx.doi.org/10.1038/nature17652) has estimated ≈0.2 Gt C of Rhizaria above 600 µm in diameter.\n",
    "\n",
    "We assume that this biomass represents the biomass of Rhizaria in mesozooplankton. As we calculated in the marine arthropod section, Rhizaria represent ≈40% of the total mesoplankton biomass:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The average fraction of Rhizaria in 18S rDNA sequencing data in surface waters is 40%\n",
      "The average fraction of Rhizaria in 18S rDNA sequencing data in the deep chlorophyll maximum is 35%\n"
     ]
    }
   ],
   "source": [
    "# Load 18S sequecing data of mesozooplankton\n",
    "seq_data = pd.read_excel('../animals/arthropods/marine_arthropods/marine_arthropods_data.xlsx',sheet_name='de Vargas',skiprows=1)\n",
    "\n",
    "print('The average fraction of Rhizaria in 18S rDNA sequencing data in surface waters is ' + '{:,.0f}%'.format(seq_data['Rhizaria surface'].mean()*100))\n",
    "print('The average fraction of Rhizaria in 18S rDNA sequencing data in the deep chlorophyll maximum is ' + '{:,.0f}%'.format(seq_data['Rhizaria DCM'].mean()*100))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The remaining 60% are made up mainly of arthropods. This would put the total mesozooplankton arthropods biomass at ≈0.3 Gt C. Our estimate for the total biomass of arthropods in the nano, micro and mesozooplankton size fraction is ≈0.56 Gt C (see the marine arthropod section for details). Subtracting the fraction of As which leaves ≈0.2 Gt C of nano and microzooplankton arthropod biomass."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# The estimate of the biomass of rhizaria based on Biard et al.\n",
    "rhizaria_biomass = 0.2e15\n",
    "\n",
    "# Calculate the average fraction of rhizaria out of the biomass\n",
    "# of mesozooplankton\n",
    "rhizaria_fraction = np.mean([seq_data['Rhizaria surface'].mean(),seq_data['Rhizaria DCM'].mean()])\n",
    "\n",
    "# Calculate the biomass of mesozooplankton arthropods\n",
    "meso_arth = rhizaria_biomass/(1-rhizaria_fraction)\n",
    "\n",
    "# Our estimate for the total biomass of arthropods in nano\n",
    "# micro and mezozooplankton size fractions\n",
    "nano_micro_mezo_arthropod = 0.56e15\n",
    "\n",
    "# Subtract the mesozooplankton arthropod biomass to estimate\n",
    "# The nanozooplankton and microzooplankton arthropod biomass\n",
    "nano_micro_arth = nano_micro_mezo_arthropod - meso_arth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the Tara Oceans data, the nano and microzooplankton arthropod biomass accounts for ≈40-75% of the total nano and microplankton biomass:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The fraction of arthropods out of the total number of reads in nanoplankton and microplankton\n",
      "The mean fraction of arthropods out of the total number of reads in nanoplankton and microplankton is ≈60%\n"
     ]
    }
   ],
   "source": [
    "print('The fraction of arthropods out of the total number of reads in nanoplankton and microplankton')\n",
    "metazoa_frac = read_data[data['Group'] == 'Metazoa']/tot_reads[['Nano reads','Micro reads']]\n",
    "\n",
    "print('The mean fraction of arthropods out of the total number of reads in nanoplankton and microplankton is ≈' + '{:,.0f}%'.format(float(metazoa_frac.mean(axis=1)*100)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the estimate we just calculated of ≈0.2 Gt C of arthropods in the nano and microplankton size fractions, and combine it with the estimate of the biomass fraction of arthropods in the nano and microplankton size fractions from the Tara Oceans dataset. This yields an estimate for the total nano and microplankton biomass which is about ≈0.4 Gt C:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The toal biomass of nano and microplankton we estimate is ≈0.4 Gt C\n"
     ]
    }
   ],
   "source": [
    "tot_nano_micro_biomass = nano_micro_arth/metazoa_frac.mean(axis=1)\n",
    "\n",
    "print('The toal biomass of nano and microplankton we estimate is ≈%.1f Gt C' %(tot_nano_micro_biomass/1e15))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we stated in the section regarding the biomass of diatoms, the biomass of nano and microplankton is estimated at ≈1 Gt C based on the MAREDAT database, which is about 2-fold larger than the estimate we got based on combination of information from Biard et al. the Tara Oceans dataset and the MAREDAT database."
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}