BNW-skin-microbiome-16S_merged.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conda activate qiime2-2021.4"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# activate qiime2-2021.4 environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir 2021_BNW_skin\n",
    "mkdir 16SV4_all\n",
    "cd 16SV4_all"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# create working folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime tools import \\\n",
    "  --type 'SampleData[PairedEndSequencesWithQuality]' \\ \n",
    "  --input-path /Users/christinanaesborg/Desktop/Raw_Samples/16SV4_BNW/Collected \\\n",
    "  --input-format CasavaOneEightSingleLanePerSampleDirFmt \\\n",
    "  --output-path demux.qza"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# import files into 16SV4 folder. Output file: demux-paired-end.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime demux summarize \\\n",
    "  --i-data demux-paired-end.qza \\\n",
    "  --o-visualization demux-paired-end.qzv\n",
    "\n",
    "qiime tools view demux-paired-end.qzv"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# create file for visualisation. Output file: demux.qzv. Visualise on https://view-qiime2.org or if on mac see code below"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Forward reads were of good quality, so we truncated it at position 240. Reverse reads were of lesser quality, therefore we truncated at position 160 were it crashed in quality."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "qiime dada2 denoise-paired \\\n",
    "--i-demultiplexed-seqs demux-paired-end.qza \\\n",
    "--p-trim-left-f 0 \\\n",
    "--p-trim-left-r 0 \\\n",
    "--p-trunc-len-f 240 \\\n",
    "--p-trunc-len-r 160 \\\n",
    "--o-representative-sequences BNW-skin-rep-seqs.qza \\\n",
    "--o-table BNW-skin-table.qza \\\n",
    "--o-denoising-stats BNW-skin-stats-dada2.qza \\\n",
    "--p-n-threads 2 "
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# dada2 is used to trim low quality bp, generated feature table and representative sequences. Output files:\n",
    "Saved FeatureTable[Frequency] to BNW-skin-table.qza\n",
    "Saved FeatureData[Sequence] to BNW-skin-rep-seqs.qza\n",
    "Saved DenoiseStats to BNW-skin-stats-dada2.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime metadata tabulate \\\n",
    "--m-input-file BNW-skin-stats-dada2.qza \\\n",
    "--o-visualization BNW-skin-stats-dada2.qzv \n",
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-table.qza \\\n",
    "--o-visualization BNW-skin-table.qzv\n",
    "qiime feature-table tabulate-seqs \\\n",
    "--i-data BNW-skin-rep-seqs.qza \\\n",
    "--o-visualization BNW-skin-rep-seqs.qzv"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# create visualisation files for all new .qza files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# since sample id were only numbers, I had to reindex. Otherwise I would have incountered problems later on\n",
    "qiime feature-table group \\\n",
    "--i-table BNW-skin-table.qza \\\n",
    "--p-axis sample \\\n",
    "--m-metadata-file BNW-skin-metadata.tsv \\\n",
    "--m-metadata-column SampleID-new \\\n",
    "--p-mode sum \\\n",
    "--o-grouped-table BNW-skin-reindexed-table.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-reindexed-table.qza \\\n",
    "--m-sample-metadata-file BNW-skin-metadata.tsv \\\n",
    "--o-visualization BNW-skin-reindexed-table.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Training SILVA database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import Qiime2 artifacts\n",
    "qiime tools import \\\n",
    "  --type 'FeatureData[Sequence]' \\\n",
    "  --input-path /Users/christinanaesborg/Desktop/SILVA_132_QIIME_release/rep_set/rep_set_16S_only/99/silva_132_99_16S.fna  \\\n",
    "  --output-path silva_99_otus.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime tools import \\\n",
    "  --type 'FeatureData[Taxonomy]' \\\n",
    "  --input-format HeaderlessTSVTaxonomyFormat \\\n",
    "  --input-path /Users/christinanaesborg/Desktop/SILVA_132_QIIME_release/taxonomy/16S_only/99/taxonomy_all_levels.txt  \\        \n",
    "  --output-path silva_99_ref-taxonomy.qza"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# output files: silva_99_otus.qza silva_99_ref-taxonomy.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract reference reads\n",
    "qiime feature-classifier extract-reads \\\n",
    "  --i-sequences silva_99_otus.qza \\\n",
    "  --p-f-primer GTGCCAGCMGCCGCGGTAA \\\n",
    "  --p-r-primer GGACTACHVGGGTWTCTAAT \\\n",
    "  --p-trunc-len 253  \\\n",
    "  --p-min-length 100 \\\n",
    "  --p-max-length 400 \\\n",
    "  --o-reads silva_ref-seqs.qza"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# output file: silva_ref-seqs.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# train naive bayes classifier\n",
    "qiime feature-classifier fit-classifier-naive-bayes \\\n",
    "  --i-reference-reads silva_ref-seqs.qza \\\n",
    "  --i-reference-taxonomy silva_99_ref-taxonomy.qza \\\n",
    "  --o-classifier silva_132_99_classifier.qza"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# output file: silva_132_99_classifier.qza"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Assign taxonomy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-classifier classify-sklearn \\\n",
    "  --i-classifier silva_132_99_classifier.qza \\\n",
    "  --i-reads BNW-skin-rep-seqs.qza \\\n",
    "  --o-classification BNW-skin-silva-132.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime metadata tabulate \\\n",
    "--m-input-file BNW-skin-silva-132.qza \\\n",
    "--o-visualization BNW-skin-silva-132.qzv"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Can take a few hours to run. Output files: BNW-skin-silva-132.qza BNW-skin-silva-132.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pre-trained SILVA database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-classifier classify-sklearn \\\n",
    "--i-classifier silva-138-99-515-806-nb-classifier.qza \\\n",
    "--i-reads BNW-skin-rep-seqs.qza \\\n",
    "--o-classification BNW-skin-silva-138-taxonomy.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime metadata tabulate \\\n",
    "--m-input-file BNW-skin-silva-138-taxonomy.qza \\\n",
    "--o-visualization BNW-skin-silva-138-taxonomy.qzv"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "There is not a huge difference (if at all) between self-trained and pre-trained database eventhough the pretrained is an updated version of the silva database. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create phylogenetic tree using sepp. Took a couple of days to run. \n",
    "#Output files: BNW-skin-sepp-tree.qza BNW-skin-sepp-placements.qza\n",
    "qiime fragment-insertion sepp \\\n",
    "--i-representative-sequences BNW-skin-rep-seqs.qza \\\n",
    "--i-reference-database sepp-refs-silva-128.qza \\\n",
    "--o-tree BNW-skin-sepp-tree.qza \\\n",
    "--o-placements BNW-skin-sepp-placements.qza \\\n",
    "--verbose \\\n",
    "--p-threads 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime fragment-insertion filter-features \\\n",
    "--i-table BNW-skin-table.qza \\\n",
    "--i-tree BNW-skin-sepp-tree.qza \\\n",
    "--o-filtered-table BNW-skin-filtered_table.qza \\\n",
    "--o-removed-table BNW-skin-removed_table.qza \\\n",
    "--verbose"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# filter out features that couldn't be inserted on the tree. Output files: BNW-skin-filtered_table.qza BNW-skin-removed_table.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-removed_table.qza \\\n",
    "--o-visualization BNW-skin-removed_table.qzv"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# visualise table. Output files: BNW-skin-removed_table.qzv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Looks like nothing was removed, so will just use the unfiltered table (BNW-skin-table.qza) for subsequent analyses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime diversity alpha-rarefaction \\\n",
    "--i-table BNW-skin-reindexed-table.qza \\\n",
    "--m-metadata-file /Users/christinanaesborg/Desktop/Microbiome/BNW-skin-metadata.tsv \\\n",
    "--o-visualization BNW-skin-reindexed-table-rarefaction.qzv \\\n",
    "--p-min-depth 500 \\\n",
    "--p-max-depth 50000\n",
    "\n",
    "# output: BNW-skin-reindexed-table-rarefaction.qzv"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Plateaus around 25K reads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table filter-samples \\\n",
    "--i-table BNW-skin-reindexed-table.qza \\\n",
    "--m-metadata-file BNW-skin-metadata.tsv \\\n",
    "--p-min-frequency 12201 \\\n",
    "--o-filtered-table BNW-table-atleast-12201-reads.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-table-atleast-12201-reads.qza \\\n",
    "--o-visualization BNW-table-atleast-12201-reads.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Remove toe samples and groups of Uncertain/low mange and Observational healthy from analysis as well as controls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table filter-samples \\\n",
    "--i-table BNW-skin-reindexed-table.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--p-where \"[bodysite] IN ('toe') OR [Severity] IN ('Uncertain/low mange', 'Observational healthy')\" \\\n",
    "--p-exclude-ids \\\n",
    "--o-filtered-table BNW-skin-table_Final-sev.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table filter-samples \\\n",
    "--i-table BNW-skin-table_Final-sev.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--p-where \"[sampletype]='Control'\" \\\n",
    "--p-exclude-ids \\\n",
    "--o-filtered-table BNW-skin-table_Final.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove samples with low reads\n",
    "qiime feature-table filter-samples \\\n",
    "--i-table BNW-skin-table_Final.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--p-min-frequency 12201 \\\n",
    "--o-filtered-table BNW-skin-table-at-least-12201-reads_Final.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-table-at-least-12201-reads_Final.qza \\\n",
    "--o-visualization BNW-skin-table-at-least-12201-reads_Final.qzv "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create taxa barplot per individual sample\n",
    "qiime taxa barplot \\\n",
    "--i-table BNW-skin-table-at-least-12201-reads_Final.qza \\\n",
    "--i-taxonomy BNW-skin-silva-138-taxonomy.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--o-visualization BNW-skin-table-at-least-12201-reads-Final_BARPLOT-PER-SAMPLE.qzv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# collapse samples by severity and make barplots\n",
    "qiime feature-table group \\\n",
    "--i-table BNW-skin-table-at-least-12201-reads_Final.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--m-metadata-column Severity \\\n",
    "--p-mode mean-ceiling \\\n",
    "--p-axis sample \\\n",
    "--o-grouped-table BNW-skin-table-SEVERITY_Final.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime taxa barplot \\\n",
    "--i-table BNW-skin-table-SEVERITY_Final.qza \\\n",
    "--i-taxonomy BNW-skin-silva-138-taxonomy.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-SEVERITY.tsv \\\n",
    "--o-visualization BNW-skin-table-Final-SEVERITY-BARPLOT.qzv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# core diversity on filtered table using 13084 rarefaction depth\n",
    "qiime diversity core-metrics-phylogenetic \\\n",
    "--i-table BNW-skin-table-at-least-12201-reads_Final.qza \\\n",
    "--i-phylogeny BNW-skin-sepp-tree.qza \\\n",
    "--p-sampling-depth 13084 \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--output-dir BNW-skin-Core-metrics-final-filtered-table-13084"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in BNW-skin-Core-metrics-final-filtered-table-13084/*vector.qza; do qiime diversity alpha-group-significance --i-alpha-diversity $i --m-metadata-file BNW-skin-metadata.tsv --o-visualization ${i/.qza/.qzv}; done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime diversity alpha \\\n",
    "  --i-table BNW-skin-table-at-least-12201-reads_Final.qza \\\n",
    "  --p-metric chao1 \\\n",
    "  --o-alpha-diversity BNW-skin-Core-metrics-final-filtered-table-13084/chao1_index.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime diversity alpha-group-significance \\\n",
    "  --i-alpha-diversity BNW-skin-Core-metrics-final-filtered-table-13084/chao1_index.qza \\\n",
    "  --m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "  --o-visualization BNW-skin-Core-metrics-final-filtered-table-13084/chao1_index.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Permanova test for microbial composition differences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime diversity beta-group-significance \\\n",
    "--i-distance-matrix BNW-skin-Core-metrics-final-filtered-table-13084/bray_curtis_distance_matrix.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--m-metadata-column Severity \\\n",
    "--o-visualization BNW-skin-Core-metrics-final-filtered-table-13084/bray_curtis_PERMANOVA.qzv \\\n",
    "--p-method permanova \\\n",
    "--p-pairwise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime diversity beta-group-significance \\\n",
    "--i-distance-matrix BNW-skin-Core-metrics-final-filtered-table-13084/jaccard_distance_matrix.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--m-metadata-column Severity \\\n",
    "--o-visualization BNW-skin-Core-metrics-final-filtered-table-13084/jaccard_PERMANOVA.qzv \\\n",
    "--p-method permanova \\\n",
    "--p-pairwise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ANCOM analysis to identify differentially abundant features/OTUS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#remove low abundance/prevalence features\n",
    "qiime feature-table filter-features \\\n",
    "--i-table BNW-skin-table-at-least-12201-reads_Final.qza \\\n",
    "--o-filtered-table BNW-skin-table-FILTERED-ANCOM-ms8-mf500_Final.qza \\\n",
    "--p-min-samples 8 \\\n",
    "--p-min-frequency 500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-table-FILTERED-ANCOM-ms8-mf500_Final.qza \\\n",
    "--o-visualization BNW-skin-table-FILTERED-ANCOM-ms8-mf500_Final.qzv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime composition add-pseudocount \\\n",
    "--i-table BNW-skin-table-FILTERED-ANCOM-ms8-mf500_Final.qza \\\n",
    "--o-composition-table BNW-skin-table-FILTERED-ANCOM-ms8-mf500-pseudo_Final.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime composition ancom \\\n",
    "--i-table BNW-skin-table-FILTERED-ANCOM-ms8-mf500-pseudo_Final.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--m-metadata-column Severity \\\n",
    "--o-visualization BNW-skin-table-FILTERED-ANCOM-ms8-mf500-pseudo-SEVERITY_Final.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ANCOM collapsed at family level (i.e. not at the feature level)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime taxa collapse \\\n",
    "--i-table BNW-skin-table-at-least-12201-reads_Final.qza  \\\n",
    "--i-taxonomy BNW-skin-silva-138-taxonomy.qza \\\n",
    "--p-level 5 \\\n",
    "--o-collapsed-table BNW-skin-table-FAMILY_Final.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table filter-features \\\n",
    "--i-table BNW-skin-table-FAMILY_Final.qza \\\n",
    "--o-filtered-table BNW-skin-table-FAMILY-FILTERED-ANCOM-ms8-mf500_Final.qza \\\n",
    "--p-min-samples 8 \\\n",
    "--p-min-frequency 500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime composition add-pseudocount \\\n",
    "--i-table BNW-skin-table-FAMILY-FILTERED-ANCOM-ms8-mf500_Final.qza  \\\n",
    "--o-composition-table BNW-skin-table-FAMILY-FILTERED-ANCOM-ms8-mf500-pseudo_Final.qza "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime composition ancom \\\n",
    "--i-table BNW-skin-table-FAMILY-FILTERED-ANCOM-ms8-mf500-pseudo_Final.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--m-metadata-column Severity \\\n",
    "--o-visualization BNW-skin-table-FAMILY-FILTERED-ANCOM-ms8-mf500-pseudo-SEVERITY_Final.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unique feature analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table filter-samples \\\n",
    "--i-table BNW-skin-table-at-least-12201-reads_Final.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--p-where \"[Severity]='Confidently healthy'\" \\\n",
    "--o-filtered-table BNW-skin-table-CONFHEALTHYonly.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-table-CONFHEALTHYonly.qza \\\n",
    "--o-visualization BNW-skin-table-CONFHEALTHYonly.qzv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table filter-samples \\\n",
    "--i-table BNW-skin-table-at-least-12201-reads_Final.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--p-where \"[Severity]='Mangy'\" \\\n",
    "--o-filtered-table BNW-skin-table-MANGYonly.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-table-MANGYonly.qza \\\n",
    "--o-visualization BNW-skin-table-MANGYonly.qzv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table filter-samples \\\n",
    "--i-table BNW-skin-table-at-least-12201-reads_Final.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--p-where \"[Severity]='Severe mange'\" \\\n",
    "--o-filtered-table BNW-skin-table-SEVEREonly.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-table-SEVEREonly.qza \\\n",
    "--o-visualization BNW-skin-table-SEVEREonly.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mitigating cross contamination/filter out common contaminats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pull out the NTCs to check if these abundant features spilled over into them as would be expected if cross-sample contamination was happening."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table filter-samples \\\n",
    "--i-table BNW-skin-reindexed-table.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--p-where \"SampleOrControl IN ('Negative control')\" \\\n",
    "--o-filtered-table BNW-skin-table-NTCs-ONLY.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-table-NTCs-ONLY.qza \\\n",
    "--o-visualization BNW-skin-table-NTCs-ONLY.qzv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime tools export \\\n",
    "--input-path BNW-skin-table-NTCs-ONLY.qza \\\n",
    "--output-path BNW-skin-table-NTCs-ONLY.biom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "biom convert -i BNW-skin-table-NTCs-ONLY.biom/feature-table.biom -o BNW-skin-table-NTCs-ONLY.biom/feature-table.tsv --to-ts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "biom convert -i bray_curtis.biom/feature-table.biom -o bray_curtis.biom/feature-table.tsv --to-tsv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# There seems to be some cross-contamination. \n",
    "# Will compare BNW-skin-reindexed-table.qza to BNW-skin-table-NTCs-ONLY.qza\n",
    "# to identify low abundance microbes as these are most likely contaminants arrising from extraction or sequencing.\n",
    "# High abundance microbes could be a transfer from true biological samples, and will be retained in further analysis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime taxa filter-table \\\n",
    "  --i-table BNW-skin-table_Final.qza \\\n",
    "  --i-taxonomy BNW-skin-silva-138-taxonomy.qza \\\n",
    "  --p-exclude mitochondria,8e09caaee61bd482522f087e96e093de,30ad487e1e86be71afce4e03fd7553b5,107273e73071282e4d0f4a0b46548da5,516a67156dacfa4341327a7e2a22f877,f2e69758d4f38addecad9194d39796e1,c489a5dd3e6b99cbbe986a7965f907ff,9ed691d9294238d730b23ab5b461835b,d829bee4984f82ffc2453212157caf96,882b5b49ef18215727a0900958cdec5d,f8485916bce858f455e79011d9a001f2,26ba1c3be9fb102889afcb9d61e730e7,06a7fe5af160203d234638f31f792cd0 \\\n",
    "  --o-filtered-table table-no-contaminants.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime metadata tabulate \\\n",
    "--m-input-file BNW-skin-silva-138-taxonomy.qza \\\n",
    "--o-visualization BNW-skin-silva-138-taxonomy.qzv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove samples with low reads\n",
    "qiime feature-table filter-samples \\\n",
    "--i-table table-no-contaminants.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--p-min-frequency 12201 \\\n",
    "--o-filtered-table BNW-skin-table-no-contaminants-at-least-12201-reads_Final.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-table-no-contaminants-at-least-12201-reads_Final.qza \\\n",
    "--o-visualization BNW-skin-table-no-contaminants-at-least-12201-reads_Final.qzv "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert to biom file for further analysis in R\n",
    "qiime tools export \\\n",
    "--input-path BNW-skin-table-no-contaminants-at-least-12201-reads_Final.qza \\\n",
    "--output-path BNW-skin-table-no-contaminants-at-least-12201-reads_Final.biom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create taxa barplot per individual sample\n",
    "qiime taxa barplot \\\n",
    "--i-table BNW-skin-table-no-contaminants-at-least-12201-reads_Final.qza \\\n",
    "--i-taxonomy BNW-skin-silva-138-taxonomy.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--o-visualization BNW-skin-table-no-cont-at-least-12201-reads-Final_BARPLOT-PER-SAMPLE.qzv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# collapse samples by severity and make barplots\n",
    "qiime feature-table group \\\n",
    "--i-table BNW-skin-table-no-contaminants-at-least-12201-reads_Final.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--m-metadata-column Severity \\\n",
    "--p-mode mean-ceiling \\\n",
    "--p-axis sample \\\n",
    "--o-grouped-table BNW-skin-table-no-cont-SEVERITY_Final.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime taxa barplot \\\n",
    "--i-table BNW-skin-table-no-cont-SEVERITY_Final.qza \\\n",
    "--i-taxonomy BNW-skin-silva-138-taxonomy.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-SEVERITY.tsv \\\n",
    "--o-visualization BNW-skin-table-Final-no-cont-SEVERITY-BARPLOT.qzv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# core diversity on filtered table using 13084 rarefaction depth\n",
    "qiime diversity core-metrics-phylogenetic \\\n",
    "--i-table BNW-skin-table-no-contaminants-at-least-12201-reads_Final.qza \\\n",
    "--i-phylogeny BNW-skin-sepp-tree.qza \\\n",
    "--p-sampling-depth 13084 \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--output-dir BNW-skin-Core-metrics-NoCont-filtered-table-13084"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in BNW-skin-Core-metrics-NoCont-filtered-table-13084/*vector.qza; do qiime diversity alpha-group-significance --i-alpha-diversity $i --m-metadata-file BNW-skin-metadata-16S.tsv --o-visualization ${i/.qza/.qzv}; done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime diversity alpha \\\n",
    "  --i-table BNW-skin-table-no-contaminants-at-least-12201-reads_Final.qza \\\n",
    "  --p-metric chao1 \\\n",
    "  --o-alpha-diversity BNW-skin-Core-metrics-NoCont-filtered-table-13084/chao1_index.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime diversity alpha-group-significance \\\n",
    "  --i-alpha-diversity BNW-skin-Core-metrics-NoCont-filtered-table-13084/chao1_index.qza \\\n",
    "  --m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "  --o-visualization BNW-skin-Core-metrics-NoCont-filtered-table-13084/chao1_index.qzv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime diversity beta-group-significance \\\n",
    "--i-distance-matrix BNW-skin-Core-metrics-NoCont-filtered-table-13084/bray_curtis_distance_matrix.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--m-metadata-column Severity \\\n",
    "--o-visualization BNW-skin-Core-metrics-NoCont-filtered-table-13084/bray_curtis_PERMANOVA.qzv \\\n",
    "--p-method permanova \\\n",
    "--p-pairwise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime diversity beta-group-significance \\\n",
    "--i-distance-matrix BNW-skin-Core-metrics-NoCont-filtered-table-13084/jaccard_distance_matrix.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--m-metadata-column Severity \\\n",
    "--o-visualization BNW-skin-Core-metrics-NoCont-filtered-table-13084/jaccard_PERMANOVA.qzv \\\n",
    "--p-method permanova \\\n",
    "--p-pairwise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unique feature analysis of severity groupings without contaminants - split into the three severity groupings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table filter-samples \\\n",
    "--i-table BNW-skin-table-no-contaminants-at-least-12201-reads_Final.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--p-where \"[Severity]='Confidently healthy'\" \\\n",
    "--o-filtered-table BNW-skin-table-NoCont-CONFHEALTHYonly.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-table-NoCont-CONFHEALTHYonly.qza \\\n",
    "--o-visualization BNW-skin-table-NoCont-CONFHEALTHYonly.qzv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table filter-samples \\\n",
    "--i-table BNW-skin-table-no-contaminants-at-least-12201-reads_Final.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--p-where \"[Severity]='Mangy'\" \\\n",
    "--o-filtered-table BNW-skin-table-NoCont-MANGYonly.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-table-NoCont-MANGYonly.qza \\\n",
    "--o-visualization BNW-skin-table-NoCont-MANGYonly.qzv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table filter-samples \\\n",
    "--i-table BNW-skin-table-no-contaminants-at-least-12201-reads_Final.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--p-where \"[Severity]='Severe mange'\" \\\n",
    "--o-filtered-table BNW-skin-table-NoCont-SEVEREonly.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-table-NoCont-SEVEREonly.qza \\\n",
    "--o-visualization BNW-skin-table-NoCont-SEVEREonly.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Provide overview of ASVs for each severity type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime tools export \\\n",
    "--input-path BNW-skin-table-NoCont-CONFHEALTHYonly.qza \\\n",
    "--output-path BNW-skin-table-NoCont-CONFHEALTHYonly.biom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "biom convert -i BNW-skin-table-NoCont-CONFHEALTHYonly.biom/feature-table.biom -o BNW-skin-table-NoCont-CONFHEALTHYonly.biom/feature-table.tsv --to-tsv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime tools export \\\n",
    "--input-path BNW-skin-table-NoCont-MANGYonly.qza \\\n",
    "--output-path BNW-skin-table-NoCont-MANGYonly.biom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "biom convert -i BNW-skin-table-NoCont-MANGYonly.biom/feature-table.biom -o BNW-skin-table-NoCont-MANGYonly.biom/feature-table.tsv --to-tsv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime tools export \\\n",
    "--input-path BNW-skin-table-NoCont-SEVEREonly.qza \\\n",
    "--output-path BNW-skin-table-NoCont-SEVEREonly.biom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "biom convert -i BNW-skin-table-NoCont-SEVEREonly.biom/feature-table.biom -o BNW-skin-table-NoCont-SEVEREonly.biom/feature-table.tsv --to-tsv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unique ASVs - min 3 samples due to severity groupings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table filter-features \\\n",
    "--i-table BNW-skin-table-NoCont-CONFHEALTHYonly.qza \\\n",
    "--o-filtered-table BNW-skin-table-NoCont-CONFHEALTHYonly-ms3.qza \\\n",
    "--p-min-samples 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-table-NoCont-CONFHEALTHYonly-ms3.qza \\\n",
    "--o-visualization BNW-skin-table-NoCont-CONFHEALTHYonly-ms3.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pre-filter\n",
    "Number of features: 2807\n",
    "Total frequency: 950,162\n",
    "\n",
    "# Post-filter\n",
    "Number of features: 280\n",
    "Total frequency: 738,131"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table filter-features \\\n",
    "--i-table BNW-skin-table-NoCont-MANGYonly.qza \\\n",
    "--o-filtered-table BNW-skin-table-NoCont-MANGYonly-ms3.qza \\\n",
    "--p-min-samples 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-table-NoCont-MANGYonly-ms3.qza \\\n",
    "--o-visualization BNW-skin-table-NoCont-MANGYonly-ms3.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pre-filter \n",
    "Number of features: 870 Total frequency: 374,474\n",
    "# Post-filter \n",
    "Number of features: 34 Total frequency: 151,855"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table filter-features \\\n",
    "--i-table BNW-skin-table-NoCont-SEVEREonly.qza \\\n",
    "--o-filtered-table BNW-skin-table-NoCont-SEVEREonly-ms3.qza \\\n",
    "--p-min-samples 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table summarize \\\n",
    "--i-table BNW-skin-table-NoCont-SEVEREonly-ms3.qza \\\n",
    "--o-visualization BNW-skin-table-NoCont-SEVEREonly-ms3.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prefilter \n",
    "Number of features: 496 Total frequency: 736,726\n",
    "# Post-filter \n",
    "Number of features: 96 Total frequency: 583,522"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Merge groups for taxa barplots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table group \\\n",
    "--i-table BNW-skin-table-NoCont-CONFHEALTHYonly-ms3.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--m-metadata-column wombatid \\\n",
    "--p-mode mean-ceiling \\\n",
    "--p-axis sample \\\n",
    "--o-grouped-table BNW-skin-table-NoCont-CONFHEALTHY-ms3-wombat-COLLAPSED.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table group \\\n",
    "--i-table BNW-skin-table-NoCont-MANGYonly-ms3.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--m-metadata-column wombatid \\\n",
    "--p-mode mean-ceiling \\\n",
    "--p-axis sample \\\n",
    "--o-grouped-table BNW-skin-table-NoCont-MANGY-ms3-wombat-COLLAPSED.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table group \\\n",
    "--i-table BNW-skin-table-NoCont-SEVEREonly-ms3.qza \\\n",
    "--m-metadata-file BNW-skin-metadata-16S.tsv \\\n",
    "--m-metadata-column wombatid \\\n",
    "--p-mode mean-ceiling \\\n",
    "--p-axis sample \\\n",
    "--o-grouped-table BNW-skin-table-NoCont-SEVERE-ms3-wombat-COLLAPSED.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qiime feature-table merge \\\n",
    "--i-tables BNW-skin-table-NoCont-CONFHEALTHY-ms3-wombat-COLLAPSED.qza BNW-skin-table-NoCont-MANGY-ms3-wombat-COLLAPSED.qza BNW-skin-table-NoCont-SEVERE-ms3-wombat-COLLAPSED.qza \\\n",
    "--o-merged-table BNW-skin-table-NoCont-wombat-COLLASPED.qza"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}