6_discriminatory_taxa.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###################################################\n",
    "# Taxa significantly contributing to classification \n",
    "###################################################\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol>\n",
       "\t<li>TRUE</li>\n",
       "\t<li>TRUE</li>\n",
       "\t<li>TRUE</li>\n",
       "\t<li>TRUE</li>\n",
       "\t<li>TRUE</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate}\n",
       "\\item TRUE\n",
       "\\item TRUE\n",
       "\\item TRUE\n",
       "\\item TRUE\n",
       "\\item TRUE\n",
       "\\end{enumerate}\n"
      ],
      "text/markdown": [
       "1. TRUE\n",
       "2. TRUE\n",
       "3. TRUE\n",
       "4. TRUE\n",
       "5. TRUE\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[[1]]\n",
       "[1] TRUE\n",
       "\n",
       "[[2]]\n",
       "[1] TRUE\n",
       "\n",
       "[[3]]\n",
       "[1] TRUE\n",
       "\n",
       "[[4]]\n",
       "[1] TRUE\n",
       "\n",
       "[[5]]\n",
       "[1] TRUE\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message in dir.create(\"output/6_discriminatory_taxa\"):\n",
      "“'output/6_discriminatory_taxa' already exists”"
     ]
    }
   ],
   "source": [
    "# Set seed \n",
    "set.seed(9356947)\n",
    "\n",
    "# Load required libraries \n",
    "source(\"scr/functions/general/load_abs_install_pkg.R\")\n",
    "\n",
    "load_abs_install_pkg(c(\"phyloseq\", \"tidyverse\", \"rfPermute\", \"ComplexHeatmap\", \"gmodels\"))\n",
    "\n",
    "# Create output directory \n",
    "dir.create(\"output/6_discriminatory_taxa\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1. Build random forest model with permutations\n",
    "#    rfPremute package\n",
    "################################################\n",
    "\n",
    "# Read in filtered and normalized phyloseq object \n",
    "ps.tf.css.01 <- readRDS(\"output/3_filt_norm_phyl/ps_tf2_css.RDS\")\n",
    "\n",
    "# Load custom function to format data from phyloseq for RF analysis \n",
    "source(\"scr/functions/data_for_rf2.R\")\n",
    "\n",
    "# Format data for RF regression analysis\n",
    "rf.data.css <- data_for_rf2(phyloseq = ps.tf.css.01, \n",
    "                              class.column = 'WeightedScoreII', \n",
    "                              remove.taxa.prev.less.than = 1)\n",
    "\n",
    "# Change name of the \"WeightedScoreII\" to \"Shedder\"\n",
    "colnames(rf.data.css)[colnames(rf.data.css) %in% \"WeightedScoreII\"] <- \"Shedder\"\n",
    "\n",
    "# Define all samples with Shedding score less than 0.51 as \"Low\" Shedders \n",
    "rf.data.css$Shedder <- as.factor(ifelse(rf.data.css$Shedder < 0.51, \"Low\", \"High\"))\n",
    "\n",
    "# Copy \"rf.data.css\" object into a new object \n",
    "rf.data.perm <- rf.data.css\n",
    "\n",
    "# Change row names \n",
    "rownames(rf.data.perm) <- paste0(\"Sample_\", 1:nrow(rf.data.css))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# 2. Build the random forest model with optimal parameters\n",
    "#    and find taxa significantly contributing to classification \n",
    "############################################################\n",
    "\n",
    "set.seed(43957)\n",
    "rf.perm.obj <- rfPermute(Shedder ~ ., \n",
    "                         data = rf.data.perm, \n",
    "                         importance=TRUE,\n",
    "                         proximity=TRUE, \n",
    "                         mtry = 565,\n",
    "                         ntree= 151, \n",
    "                         nrep = 9,\n",
    "                         num.cores = 12)\n",
    "\n",
    "# Save the resulted object \n",
    "#save(rf.perm.obj, file = \"output/6_discriminatory_taxa/srf_perm_obj.RData\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# 3. Extract features that are significantly affecting mean accuracy and the Gini index\n",
    "#######################################################################################\n",
    "\n",
    "# Extract P values from RF Premute object in as a dataframe \n",
    "pval.rf.perm <- data.frame(rf.perm.obj$pval)\n",
    "\n",
    "# Subset taxa significantly contributing to MeanDecreaseAccuracy (p < 0.05)\n",
    "pval.rf.perm.s0 <- pval.rf.perm[pval.rf.perm$MeanDecreaseAccuracy.scaled < 0.05, ]\n",
    "\n",
    "# Subset taxa significantly contributing to MeanDecreaseGini (p < 0.05)\n",
    "pval.rf.perm.s1 <- pval.rf.perm.s0[pval.rf.perm.s0$MeanDecreaseGini.scaled < 0.05, ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>png:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{png:} 2"
      ],
      "text/markdown": [
       "**png:** 2"
      ],
      "text/plain": [
       "png \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "# 4. Visualize abundance of significantly contributing taxa as a heatmap \n",
    "########################################################################\n",
    "\n",
    "# Extract data about normalized count of significantly contributing taxa \n",
    "abund.heat <- rf.data.perm[, colnames(rf.data.perm) %in% rownames(pval.rf.perm.s1)]\n",
    "\n",
    "# Order rows by shedding status of corresponding cows  \n",
    "abund.heat <- abund.heat[order(rf.data.perm$Shedder), ]\n",
    "\n",
    "# Order columns by total abundance\n",
    "abund.heat <- abund.heat[,order(colSums(abund.heat)) ]\n",
    "\n",
    "# Prepare data for colored bar \n",
    "col.stat <- c(\"steelblue\", \"gold3\")\n",
    "\n",
    "# Adjust names in bar data \n",
    "names(col.stat) <- c(\"High\", \"Low\")\n",
    "\n",
    "# Extract data about shedding status \n",
    "heat.ant.0 <- data.frame(as.character(rf.data.perm$Shedder))\n",
    "\n",
    "# Add Cow IDs as a column for shedding status \n",
    "heat.ant.0$rnam <- rownames(rf.data.perm)\n",
    "\n",
    "# Add Cow IDs as row names  \n",
    "rownames(heat.ant.0) <- rownames(rf.data.perm)\n",
    "\n",
    "heat.ant.0 <- heat.ant.0[rownames(abund.heat), ]\n",
    "\n",
    "heat.ant <- data.frame(heat.ant.0[,1])\n",
    "\n",
    "# Adjust column name of heat annotation df \n",
    "colnames(heat.ant) <- \"Shedding\"\n",
    "\n",
    "# Add Cow IDs as row names for heat annotation df \n",
    "rownames(heat.ant) <- heat.ant.0$rnam\n",
    "\n",
    "# Combine annotation \n",
    "ha =  HeatmapAnnotation(df = heat.ant, col = list(Shedding = col.stat))\n",
    "\n",
    "# Plot the heatmap \n",
    "set.seed(69734)\n",
    "g.heat.abound <- ComplexHeatmap::Heatmap(t(abund.heat), \n",
    "                        top_annotation = ha, \n",
    "                        show_column_names = FALSE, \n",
    "                        #km = 2,\n",
    "                        show_row_names = FALSE,\n",
    "                        #col = my_palette, \n",
    "                        cluster_columns = FALSE,\n",
    "                        cluster_rows = TRUE,\n",
    "                        name = \"Abundance\" )\n",
    "\n",
    "# Save the heatmap plot \n",
    "pdf(\"output/6_discriminatory_taxa/Figure_5A.pdf\", width = 6, height = 3.5, paper='special')\n",
    "g.heat.abound\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# 5. Create a local boxplot function \n",
    "####################################\n",
    "\n",
    "box_sig <- function(data.matrix, heat.ant, prev) {\n",
    "    \n",
    "    # Prepare data \n",
    "    # If ploting prevalence plot (prevalence data supplyed)\n",
    "    if (prev == TRUE) {\n",
    "        \n",
    "        # Calculate prevalence in the \"High\" shedding group by finding sum for each taxa in all samples\n",
    "        #           and dividing by the number of samples\n",
    "        High <- (as.vector(colSums(data.matrix[heat.ant$Shedding %in% \"High\", ]))/\n",
    "                  length(heat.ant$Shedding[heat.ant$Shedding %in% \"High\"]))*100\n",
    "    \n",
    "        # Calculate prevalence in the \"High\" shedding group by finding sum for each taxa in all samples\n",
    "        #           and dividing by the number of samples \n",
    "        Low <- (as.vector(colSums(data.matrix[heat.ant$Shedding %in% \"Low\", ]))/\n",
    "                  length(heat.ant$Shedding[heat.ant$Shedding %in% \"Low\"]))*100\n",
    "     \n",
    "    # If plotting abundance (abundance data supplied)\n",
    "    } else {\n",
    "        \n",
    "        # Find taxa mean abundance in the \"High\" shedding group \n",
    "        High <- as.vector(colMeans(data.matrix[heat.ant$Shedding %in% \"High\", ]))\n",
    "        \n",
    "        # Find taxa mean abundance in the \"Low\" shedding group \n",
    "        Low <- as.vector(colMeans(data.matrix[heat.ant$Shedding %in% \"Low\", ]))   \n",
    "        \n",
    "    }\n",
    "\n",
    "    # Create a datatable \n",
    "    d <- data.frame(High = High, Low = Low)\n",
    "\n",
    "    # Add samples names as column names \n",
    "    d$Taxa <- colnames(data.matrix)[1:ncol(data.matrix)]\n",
    "\n",
    "    # Add column with an indication if taxa decrease or increase in the groups  \n",
    "    d$colr <- ifelse(d$High <= d$Low, \"Increase\", \"Dicrease\")\n",
    "\n",
    "    # Melt the dataframe \n",
    "    d.m <- reshape2::melt(d) \n",
    "\n",
    "    # Convert column \"Taxa\" into a factor \n",
    "    d.m$Taxa <- factor(d.m$Taxa)\n",
    "\n",
    "    # Add value for a proper log transformation if abundance is plotted\n",
    "    if (prev == FALSE) {\n",
    "        d.m$value <- d.m$value + 0.01\n",
    "    }\n",
    "\n",
    "\n",
    "    # Plot the results for prevalence mode \n",
    "    if (prev == TRUE) {\n",
    "        \n",
    "        p.box <- ggplot(d.m, aes(x = variable, y = value)) +   \n",
    "                  geom_point() + \n",
    "                  geom_line(aes(group = Taxa, color=colr), size = 0.8, alpha=0.75) +   \n",
    "                  geom_boxplot(alpha=0.01) + theme_bw() + theme(legend.position = \"none\") + \n",
    "                  ylab(\"Prevalence (%)\") + \n",
    "                  xlab(\"Status\")\n",
    "    \n",
    "    # Plot the results for abundance mode (log transformed)     \n",
    "    } else {\n",
    "        \n",
    "        p.box <- ggplot(d.m, aes(x = variable, y = log(value))) +   \n",
    "                  geom_point() + \n",
    "                  geom_line(aes(group = Taxa, color=colr), size = 0.8, alpha=0.75) +   \n",
    "                  geom_boxplot(alpha=0.01) + theme_bw() + theme(legend.position = \"none\") + \n",
    "                  ylab(\"Abundance (Log)\") + \n",
    "                  xlab(\"Status\")\n",
    "        \n",
    "    }\n",
    "  \n",
    "# Return the boxplot    \n",
    "return(p.box)\n",
    "    \n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using Taxa, colr as id variables\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# 6. Plot abundance difference of taxa significantly contributes to classification \n",
    "#         between \"High\" and \"Low\" shedding groups \n",
    "#####################################################################################\n",
    "\n",
    "# Generate abundance boxplot using function above \n",
    "sig.abund.box <- box_sig(data.matrix = abund.heat, heat.ant = heat.ant, prev = FALSE)\n",
    "\n",
    "# Save the plot \n",
    "ggsave(\"output/6_discriminatory_taxa/Figure_5B.pdf\", sig.abund.box, width = 3, height = 5)\n",
    "ggsave(\"output/6_discriminatory_taxa/Figure_5B.png\", sig.abund.box, width = 3, height = 5, dpi = 300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using Taxa, colr as id variables\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# 7. Plot prevalence difference of taxa significantly contributes to classification \n",
    "#         between \"High\" and \"Low\" shedding groups \n",
    "####################################################################################\n",
    "\n",
    "# Copy heatmap data in an object \n",
    "prev.heat <- abund.heat\n",
    "\n",
    "# Convert into prevalence \n",
    "prev.heat[prev.heat > 0] <- 1 \n",
    "\n",
    "# Generate prevalence boxplot using function above \n",
    "sig.prev.box <- box_sig(data.matrix = prev.heat, heat.ant = heat.ant, prev = TRUE)\n",
    "\n",
    "# Save the plot \n",
    "ggsave(\"output/6_discriminatory_taxa/Figure_5C.pdf\", sig.prev.box, width = 3, height = 5)\n",
    "ggsave(\"output/6_discriminatory_taxa/Figure_5C.png\", sig.prev.box,  width = 3, height = 5, dpi = 400)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# 8. Format taxonomy of significantly contributing taxa \n",
    "#######################################################\n",
    "\n",
    "# Extract a taxonomic table from the phyloseq object  \n",
    "tax.table <- tax_table(ps.tf.css.01)\n",
    "\n",
    "# Subset significantly contributing taxa \n",
    "sig.tax <- data.frame(tax.table[rownames(tax.table) %in% colnames(abund.heat), ])\n",
    "\n",
    "# Convert every column into character class \n",
    "sig.tax.m <- mutate_all(sig.tax, as.character)\n",
    "\n",
    "# Create vector with prefixes for taxonomic level\n",
    "phy.ind <- c(\"k_\", \"p_\", \"c_\", \"o_\", \"f_\", \"g_\")\n",
    "\n",
    "# Add prefixes to taxonomic levels \n",
    "for (i in 1:ncol(sig.tax.m)) {\n",
    "    \n",
    "    sig.tax.m[, i] <- paste0(phy.ind[i], sig.tax.m[, i])\n",
    "    \n",
    "}\n",
    "\n",
    "# Add row names \n",
    "rownames(sig.tax.m) <- rownames(sig.tax)\n",
    "\n",
    "# Add column GenusID \n",
    "sig.tax.m$GenusID <- sig.tax.m$Genus\n",
    "\n",
    "# Format GenusID column \n",
    "# \"i\" is a row number \n",
    "for (i in 1:nrow(sig.tax.m)) {\n",
    "        \n",
    "        # Subset a row \n",
    "        t.na <- sig.tax.m[i,]\n",
    "        \n",
    "        # Grep objects that are not \"NA\"\n",
    "        t.na <- t.na[!grepl(\"NA\", t.na)]\n",
    "        \n",
    "        # Grep objects that are not \"-\"\n",
    "        t.na <- t.na[!grepl(\"-\", t.na)]\n",
    "        \n",
    "        # Select last object (last identified taxonomic level) and add \n",
    "        #         ASV \n",
    "        sig.tax.m$GenusID[i] <- paste0(t.na[length(t.na)], \"-\", \"ASV\") \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# 9. Summary of information about taxa significant for classification.\n",
    "######################################################################\n",
    "\n",
    "# Subset data about abundance of significantly contributing taxa \n",
    "#        from the data used for RF model with permutations\n",
    "sig.tax.abund <- rf.data.perm[, rownames(sig.tax.m)]\n",
    "\n",
    "# Copy to make a prevalence object \n",
    "sig.tax.prev <- sig.tax.abund\n",
    "\n",
    "# Convert abundance into prevalence in the prevalence object\n",
    "sig.tax.prev[sig.tax.prev > 0] <- 1 \n",
    "\n",
    "# Calculate mean abundance in the \"Low\" shedding group \n",
    "sig.tax.m$MeanAdunanceLow <- colMeans(sig.tax.abund[rf.data.perm$Shedder %in% \"Low\", ])\n",
    "\n",
    "# Calculate mean abundance in the \"High\" shedding group \n",
    "sig.tax.m$MeanAdunanceHigh <- colMeans(sig.tax.abund[rf.data.perm$Shedder %in% \"High\", ])\n",
    "\n",
    "# Calculate prevalence expressed in % in the \"Low\" shedding group\n",
    "sig.tax.m$PrevalenceLow <- colMeans(sig.tax.prev[rf.data.perm$Shedder %in% \"Low\", ]) * 100\n",
    "\n",
    "# Calculate prevalence expressed in % in the \"High\" shedding group \n",
    "sig.tax.m$PrevalenceHigh <- colMeans(sig.tax.prev[rf.data.perm$Shedder %in% \"High\", ]) * 100\n",
    "\n",
    "# Make an empty vector for significance testing of differences in \n",
    "#      abundance between \"Low\" and \"High\" shedding group\n",
    "wt.prev <- c()\n",
    "\n",
    "# Test difference significance (loop)\n",
    "# \"i\" number of a column \n",
    "for (i in 1:ncol(sig.tax.abund)) {\n",
    "   \n",
    "   # Perform Wilcoxon test  \n",
    "   pw <- wilcox.test(sig.tax.abund[rf.data.perm$Shedder %in% \"Low\", i], \n",
    "                     sig.tax.abund[rf.data.perm$Shedder %in% \"High\", i])\n",
    "    \n",
    "   # Combine obtained P values into a vector   \n",
    "   wt.prev <- c(wt.prev, pw$p.value)  \n",
    "}\n",
    "\n",
    "# Add P values to the summary table \n",
    "sig.tax.m$WilcoxAbundPval <- wt.prev \n",
    "\n",
    "# Add adjusted P values to the summary table \n",
    "sig.tax.m$WilcoxAbundQvalFDR <- p.adjust(wt.prev, method = \"fdr\")\n",
    "\n",
    "# Add information about contribution to the classification model\n",
    "rf.imp.tax <- rf.perm.obj$importance\n",
    "\n",
    "# Subset only significantly contributing taxa \n",
    "rf.imp.tax <- rf.imp.tax[rownames(rf.imp.tax) %in% rownames(sig.tax.m), ]\n",
    "\n",
    "# Combine data from both tables \n",
    "sig.tax.f <- cbind(sig.tax.m, rf.imp.tax[, c(3:4)])\n",
    "\n",
    "# Make GenusID unique\n",
    "sig.tax.f$GenusID <- make.unique(sig.tax.f$GenusID)\n",
    "\n",
    "# Save the final table \n",
    "write.table(sig.tax.f, \"output/6_discriminatory_taxa/sig_contr_taxa.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 10. Plot Mean Decrease Accuracy and Gini indexes\n",
    "###################################################\n",
    "\n",
    "# Prepare data \n",
    "# Subset data for plotting from table with information about significant contributing taxa\n",
    "sig.p.d <- sig.tax.f[, c(\"Phylum\", \"GenusID\", \"MeanDecreaseAccuracy\", \"MeanDecreaseGini\")]\n",
    "\n",
    "# Order taxa by mean decrease accuracy \n",
    "sig.p.d <- sig.p.d[order(sig.p.d$MeanDecreaseAccuracy), ]\n",
    "\n",
    "# Subset only top 20 contributing taxa \n",
    "sig.p.d <- sig.p.d[(nrow(sig.p.d) - 20):nrow(sig.p.d), ]\n",
    "\n",
    "# Melt (gather) data into long format \n",
    "sig.p.dm <- gather(sig.p.d, Index, Value, -Phylum, -GenusID)\n",
    "\n",
    "# Arrange appropriate order of level in GenusID\n",
    "sig.p.dm$GenusID <- factor(sig.p.dm$GenusID, levels=sig.p.d$GenusID)\n",
    "\n",
    "# Plot and save \n",
    "mda.plot <- ggplot(sig.p.dm, aes(x = GenusID, y = Value, fill = Phylum)) + \n",
    "            geom_bar(stat = \"identity\", width = 0.6, color = \"black\") +\n",
    "            facet_grid(.~ Index, scales = \"free\") + \n",
    "            theme_bw() + \n",
    "            theme(axis.text.x = element_text(angle = 45, hjust = 1)) + \n",
    "            ylab(\"\") + \n",
    "            xlab(\"AVSs with the corresponding genus names\") +\n",
    "            scale_fill_brewer(palette = \"Dark2\") + \n",
    "            coord_flip()\n",
    "\n",
    "ggsave(filename = \"output/6_discriminatory_taxa/Figure_4.pdf\", plot = mda.plot, width = 8, height = 5.5)\n",
    "ggsave(filename = \"output/6_discriminatory_taxa/Figure_4.jpg\", plot = mda.plot, width = 8, height = 5.5, dpi = 400)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message in ci.numeric(sig.tax.m$MeanAdunanceLow):\n",
      "“No class or unkown class.  Using default calcuation.”Warning message in ci.numeric(sig.tax.f$MeanAdunanceHigh):\n",
      "“No class or unkown class.  Using default calcuation.”"
     ]
    }
   ],
   "source": [
    "# 11. Differences in abundance and prevalence between \"Low\" & \"High\" shedding groups\n",
    "####################################################################################\n",
    "# Will gather summary and statistics used in the manuscript text into a single object (list)\n",
    "\n",
    "all.summ <- list()\n",
    "\n",
    "all.summ[[\"abund.dif\"]] <- summary(abs(sig.tax.m$MeanAdunanceLow - sig.tax.f$MeanAdunanceHigh))\n",
    "\n",
    "all.summ[[\"prev.dif\"]] <- summary(abs(sig.tax.m$PrevalenceLow - sig.tax.f$PrevalenceHigh))\n",
    "\n",
    "all.summ[[\"abund.dif.test\"]] <- wilcox.test(sig.tax.m$MeanAdunanceLow, sig.tax.f$MeanAdunanceHigh)\n",
    "\n",
    "all.summ[[\"prev.dif.test\"]] <-wilcox.test(sig.tax.m$PrevalenceLow, sig.tax.f$PrevalenceHigh)\n",
    "\n",
    "all.summ[[\"abund.low.ci\"]] <- ci(sig.tax.m$MeanAdunanceLow)\n",
    "\n",
    "all.summ[[\"abund.high.ci\"]] <- ci(sig.tax.f$MeanAdunanceHigh)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}