# part 1: import, demux, and denoising ### Background: 20 adult dairy cattle were sampled (via swabbing) to determine the microbial community composition between 5 different body sites (nasal/skin/udder/fecal/oral).  We want you to determine if and how the microbial composition changes between sites. We will first begin by downloading the raw sequencing data from the workshop website and then upload those files onto Alpine into a new directory (that you create). Workshop website/downloadable files: https://q2-workshop-432f96.gitlab.io/latest/day3/ and click the "Download the raw_reads.zip directory" **Steps:** 1. Log into Alpine using OnDemand and create a new directory for this new analysis in your _scratch_ directory. Hint: it should look something like: /scratch/alpine/$USER/cow/ 2. Move into your new directory using OnDemand 3. Create the following sub-directories using OnDemand: 1. slurm 2. taxonomy 3. tree 4. taxaplots 5. dada2 6. demux 7. metadata 4. Download the cow_barcodes.txt and cow_metadata.txt files from the website and upload them both to your metadata folder within your new cow directory. So, the filepath where you puts these files should look something like this: /scratch/alpine/$USER/cow/metadata 5. Upload the raw_reads directory and then unzip it using the following commands: Launch an interactive session and load qiime2 within your cow directory. ``` #launch an interactive session: sinteractive --reservation=microbiome --time=04:00:00 --partition=amilan --nodes=1 --ntasks=6 --qos=normal #insert your code here to activate qiime. Hint: there should be 2 things you add here ``` Unzip the raw_sequences directory and remove the .zip file so that you can free up space ``` unzip raw_reads.zip rm -r raw_reads.zip ``` Import the sequences/reads into a Qiime2-readable format (.qza). Note this might take 10-20 mins ``` qiime tools import \ --type EMPPairedEndSequences \ --input-path raw_reads \ --output-path cow_reads.qza ``` ## Demultiplex the reads by submitting a job. Note this may take ~30 mins a.    Go into your slurm subdirectory within your cow directory using OnDemand. Create a new file named **demux.sh** so you can submit a job that will demultiplex your sequences quicker. Fill in the lines that need editing (denoted by capital letters or hashes) to this demultiplexing command and add that to your new script. ``` #!/bin/bash #SBATCH --job-name=demux #SBATCH --nodes=1 #SBATCH --ntasks=12 #SBATCH --partition=amilan #SBATCH --time=04:00:00 #SBATCH --mail-type=ALL #SBATCH --output=slurm-%j.out #SBATCH --qos=normal #SBATCH --mail-user=ADD_YOUR_EMAIL@colostate.edu #What needs to go here in order to “turn on” qiime2? Hint: we do these 2 commands every time we activate qiime2! #change the following line if your file path looks different cd /scratch/alpine/$USER/cow/demux #Below is the command you will run to demultiplex the samples. qiime demux emp-paired \ --m-barcodes-file ../metadata/cow_barcodes.txt \ --m-barcodes-column barcode \ --p-rev-comp-mapping-barcodes \ --p-rev-comp-barcodes \ --i-seqs ../cow_reads.qza \ --o-per-sample-sequences demux_cow.qza \ --o-error-correction-details cow_demux_error.qza #visualize the read quality qiime demux summarize \ --i-data demux_cow.qza \ --o-visualization demux_cow.qzv ```  Run the script in your slurm directory as a job using:  ```  sbatch demux.sh  ``` ## Denoise (update denoise command for new version of qiime) Fill in the blank to denoise your samples based on what you think the reads should be trimmed at (from the front of the reads) or truncated at (from the ends of the reads) based on the demux_cow.qzv file. You can run this in the terminal or as a job. ``` cd /scratch/alpine/$USER/cow/dada2 qiime dada2 denoise-paired \ --i-demultiplexed-seqs ../demux/demux_cow.qza \ --p-trim-left-f 0 \ --p-trim-left-r 0 \ --p-trunc-len-f ADD TRUNCATION LENGTH \ --p-trunc-len-r ADD TRUNCATION LENGTH \ --p-n-threads 6 \ --o-representative-sequences cow_seqs_dada2.qza \ --o-denoising-stats cow_dada2_stats.qza \ --o-table cow_table_dada2.qza \ --o-base-transition-stats base-transition-stats_cow.qza #Visualize the denoising results: qiime metadata tabulate \ --m-input-file cow_dada2_stats.qza \ --o-visualization cow_dada2_stats.qzv qiime feature-table summarize \ --i-table cow_table_dada2.qza \ --o-feature-frequencies feature-frequencies.qza \ --o-sample-frequencies sample-frequencies.qza \ --o-summary ADD FILE NAME.qzv qiime feature-table tabulate-seqs \ --i-data cow_seqs_dada2.qza \ --o-visualization ADD FILE NAME.qzv ``` Stop and think: 1. What is the mean reads per sample? 2. How long are the reads? 3. What is the maximum length of all your sequences? 4. Which sample (not including extraction controls starting with EC) lost the highest % of reads? 5. Why did you chose to trim or truncate where you did? __________________________________________________________________ # Cow part 2 - taxonomy, and phylogenetic tree **Learning Objectives** 1. Practice independently running a partial Qiime2 workflow on new data using both the command line and slurm to submit jobs. 2. Practice recording commands and editing code to match your analysis. 3. Perform taxonomic assignments using greengenes2, filtering out unwanted taxonomy and large amplicons, run a nother job script for generating a SEPP tree. ## Remove long (300+ base pair) amplicons from the representative sequences file and the feature table Note: we did not do this for the decomp dataset, because this was not an issue for that dataset, but it is necessary here ``` # filter out any large amplicons from the seqs and table (because they are contaminates) cd /scratch/alpine/$USER/cow/dada2 qiime feature-table filter-seqs \ --i-data cow_seqs_dada2.qza \ --m-metadata-file cow_seqs_dada2.qza \ --p-where 'length(sequence) < 300' \ --o-filtered-data cow_seqs_dada2_filtered300.qza qiime feature-table tabulate-seqs \ --i-data cow_seqs_dada2_filtered300.qza \ --o-visualization cow_seqs_dada2_filtered300.qzv qiime feature-table filter-features \ --i-table cow_table_dada2.qza \ --m-metadata-file cow_seqs_dada2_filtered300.qza \ --o-filtered-table cow_table_dada2_filtered300.qza qiime feature-table summarize \ --i-table cow_table_dada2_filtered300.qza \ --o-feature-frequencies feature-frequencies-filtered.qza \ --o-sample-frequencies sample-frequencies-filtered.qza \ --o-summary cow_table_dada2_filtered300.qzv     ``` ## Classify taxonomy using GreenGenes2 First get the Greengenes2 database: ``` cd /scratch/alpine/$USER/cow/taxonomy ``` ``` wget --no-check-certificate https://ftp.microbio.me/greengenes_release/2024.09/2024.09.backbone.v4.nb.qza ``` Classify taxonomy using GreenGenes2 classify the ASVs (takes about 5 mins). ~={red}(1point)=~ ``` qiime feature-classifier classify-sklearn \ --i-reads ../dada2/cow_seqs_dada2_filtered300.qza \ --i-classifier 2024.09.backbone.v4.nb.qza \ --o-classification taxonomy_gg2_filtered.qza ``` - visualize the taxonomy of all the ASVs ``` qiime metadata tabulate \ --m-input-file taxonomy_gg2_filtered.qza \ --o-visualization taxonomy_gg2_filtered.qzv ``` - Filter mitochondria and chloroplast out to generate a filtered feature table, keep only ASVs with a class or lower taxonomy ``` qiime taxa filter-table \ --i-table ../dada2/cow_table_dada2_filtered300.qza \ --i-taxonomy taxonomy_gg2_filtered.qza \ --p-exclude INSERT WHAT TO FILTER OUT \ --p-include c__ \ --o-filtered-table ../dada2/table_nomitochloro_gg2_filtered300.qza ``` - Visualize the taxa bar plot ``` qiime taxa barplot \ --i-table ../dada2/table_nomitochloro_gg2_filtered300.qza \ --i-taxonomy taxonomy_gg2_filtered.qza \ --m-metadata-file ../metadata/cow_metadata.txt \ --o-visualization ../taxaplots/taxa_barplot_nomitochloro_gg2_filtered300.qzv ``` ### Filtered Taxa Bar Plot Questions 1. Look at your taxa bar plot, organized by cow sampling location (body_site) at the level 7 taxonomic level. What general trends do you notice?   2. What are the top 2 most abundant bacterial **classes** in the fecal samples? **_ 3. What highly abundant ASV is shared between both the udder and skin samples? 4. Which samples (still sorted by body_site) have higher alpha diversity in terms of observed features? 5. Do all samples contain archaea as well? 6. Why do we filter out sp004296775? 7. what is the difference between these two flags? --p-exclude mitochondria,chloroplast,sp004296775 \ --p-include c__ \ 8. do the positive controls look the same as each other? Yes or No? 9. Do the negative/extraction controls (Samples labeled as EC), look like the positive controls? Yes or no? 10. do the negative/extraction controls (Samples labeled as EC), look like the real samples? Yes or no? ## Phylogenetic tree Create a job script to run the phylogenetic tree building. Remember you must start a new terminal session, navigate to your slurm directory, and then submit the job. You do NOT need to start any other interactive sessions. This job will take about an hour. Go to OnDemand and create a new text file, name it: tree.sh ``` cd /scratch/alpine/$USER/cow/slurm tree.sh ``` - Put the commands into the file, and save it. ``` #!/bin/bash #SBATCH --job-name=tree #SBATCH --nodes=1 #SBATCH --ntasks=8 #SBATCH --partition=amilan #SBATCH --time=12:00:00 #SBATCH --mail-type=ALL #SBATCH --mail-user=$USER@colostate.edu #SBATCH --output=slurm-%j.out #SBATCH --qos=normal #Activate qiime module purge module load qiime2/2026.1_amplicon #Get reference tree (gg2) wget --no-check-certificate -P ../tree https://ftp.microbio.me/greengenes_release/2022.10/2022.10.backbone.sepp-reference.qza #Command qiime fragment-insertion sepp \ --i-representative-sequences ../dada2/cow_seqs_dada2_filtered300.qza \ --i-reference-database ../tree/2022.10.backbone.sepp-reference.qza \ --o-tree ../tree/tree_gg2.qza \ --o-placements ../tree/tree_placements_gg2.qza ``` - submit the job from the terminal ``` sbatch tree.sh ``` ------------------------------------------------------------------ # **Cow Site part 3 - rarefy and run core metrics ### Alpha Rarefaction Plot - Chose the input sequencings depths (min and max) for generating the alpha rarefaction plot: ``` qiime diversity alpha-rarefaction \ --i-table dada2/table_nomitochloro_gg2_filtered300.qza \ --m-metadata-file metadata/cow_metadata.txt \ --o-visualization alpha_rarefaction_curves_16S.qzv \ --p-min-depth 10 \ --p-max-depth INPUT SEQ DEPTH ``` ### Filter out controls before core metrics ``` # Get matadata with no controls See website, and upload into the correct file (hint: where have we put metadata in the past) ``` Go back into so you are in the cow directory ``` cd /scratch/alpine/$USER/cow ``` ``` qiime feature-table filter-samples \ --i-table dada2/table_nomitochloro_gg2_filtered300.qza \ --m-metadata-file metadata/cow_metadata_nocontrols.txt \ --o-filtered-table dada2/table_nomitochlorocontrols_gg2_filtered300.qza ``` ### Run Core Metrics ``` qiime diversity core-metrics-phylogenetic \ --i-table dada2/table_nomitochlorocontrols_gg2_filtered300.qza \ --i-phylogeny tree/INSERT TREE NAME.qza \ --m-metadata-file metadata/cow_metadata.txt \ --p-sampling-depth INSERT RAREFACTION DEPTH \ --output-dir core_metrics_results ``` ### Visualize alpha diversity plots - generate a plot to visualize the observed features ``` qiime diversity alpha-group-significance \ --i-alpha-diversity core_metrics_results/observed_features_vector.qza \ --m-metadata-file metadata/cow_metadata.txt \ --o-visualization core_metrics_results/observed_features_statistics.qzv ``` - generate a plot to visualize faiths PD ``` ## insert the entire code chunk for generating this visualization ``` ### Beta diversity significance - unweighted unifrac ``` qiime diversity beta-group-significance \ --i-distance-matrix core_metrics_results/unweighted_unifrac_distance_matrix.qza \ --m-metadata-file metadata/cow_metadata.txt \ --m-metadata-column body_site \ --o-visualization core_metrics_results/unweighted_unifrac_distance_matrix.qzv ``` - Weighted unifrac ``` # put the code chunk here that you would use for the beta diversity permanova ``` ## stop and think: 1. what is the name of the file you needed to use to figure out what min and max depths to use to generate the alpha rarefaction plot? (Hint: which file contains the sequencing depths for each sample) 2. what did you chose for the rarefaction depth (the input for core metrics -p-sampling-depth flag)? why? 3. Which cow body location had more observed features? 4. Which has the lowest? 5. What is the main difference between Faiths PD and Shannons alpha diversity metrics? 6. Which two body sites have the highest Faiths PD alpha diversity?  Are the groups significantly different? 7. Which diversity metrics produced by the core-metrics pipeline require phylogenetic information? 8. Does it seem like there are any groupings in the beta diversity? What are the groupings? 9. Why do you think these samples are grouping together? 10. What test can you run to determine if the groups are significantly different? 11. What command would you use to run that test? # Part 4 - ANCOM-BC2 and Machine Learning **Filter controls out of our table - Navigate into the cow tutorial and make a new ancombc2 directory for the ANCOM-BC2 analysis - Navigate into the ancombc2 directory ``` #Make an ancombc2 dir then move into that directory (hint: what commands have we used to make a directory in the past) ``` **Filter Samples ** - Choose the min frequency for sample filtering: ``` qiime feature-table filter-samples \ --i-table ../dada2/table_nomitochlorocontrols_gg2_filtered300.qza \ --p-min-frequency YOUR NUMBER HERE \ --o-filtered-table table_5k.qza ``` **Filter out low abundance and low prevalence ASVs ** ``` qiime feature-table filter-features \ --i-table INPUT TABLE.qza \ --p-min-frequency 50 \ --p-min-samples 4 \ --o-filtered-table table_5k_abund.qza ``` **Collapse features to genus level ** - We will collapse to the genus level to make it easier to interpret the results. (Hint: We used 7 for species, so think about which number you would use for genus.) ``` qiime taxa collapse \ --i-table table_5k_abund.qza \ --i-taxonomy ../taxonomy/taxonomy_gg2_filtered.qza \ --p-level GENUS NUMBER \ --o-collapsed-table table_5k_abund_GENUS NUMBER.qza ``` **Run ANCOM-BC2 ** ``` qiime composition ancombc2 \ --i-table table_5k_abund_GENUS NUMBER.qza \ --m-metadata-file ../metadata/cow_metadata_nocontrols.txt \ --p-fixed-effects-formula body_site \ --o-ancombc2-output ancombc2_results_bodysite_genus.qza ``` **Visualize the ANCOM-BC2 results ** - Generate a barplot to visualize the differentially abundant features. ``` qiime composition tabulate \ --i-data INPUT FILE \ --o-visualization ancombc2_bodysite_genus.qzv qiime composition ancombc2-visualizer \ --i-data INPUT FILE (SAME AS ABOVE) \ --o-visualization ancombc2_barplot_bodysite_genus.qzv ``` ## Generate a classification model to see how well we can predict cow body site ``` cd /scratch/alpine/$USER/cow/ mkdir ml cd ml ``` ``` qiime taxa collapse \ --i-table ../core_metrics_results/rarefied_table.qza \ --i-taxonomy ../taxonomy/taxonomy_gg2_filtered.qza \ --p-level 7 \ --o-collapsed-table rarefied_table_L7.qza ``` ``` qiime sample-classifier classify-samples \ --i-table rarefied_table_L7.qza \ --m-metadata-file ../metadata/cow_metadata_nocontrols.txt \ --m-metadata-column WHAT TO TEST \ --p-random-state 123 \ --p-n-jobs 1 \ --output-dir sample_classifier_results_bodysite ``` #### **Stop and Think** 1. When generating the filtered table for ANCOM-BC2, what value did you choose for `--p-min-frequency`? Which core metrics parameter should this match, and why do these values need to be the same? (Report your core metrics value here: ___) 2. Why do we filter out samples with low frequency and low abundance ASVs? 3. What was the most enriched genus in skin compared to fecal, and what was the most depleted genus in skin compared to fecal (make sure adjusted p is set to less than 0.05)? 4. Why might removing controls be important before downstream analysis? 5. What 2 features that are high in fecal samples? 6. What are 2 features that are low in nasal? 7. What is the accuracy of your model, and if the accuracy of the classifier is high, what does that suggest about the microbial compositions of each site? ### Cow R tutorial #### Reference if you want to work through some R on your own after the workshop **Set up the cow R analysis file structure** - Make a cow_r directory on your local computer, and inside the cow_r directory, make the following directories cow_r ├── 01_notes ├── 02_data ├── 03_metadata ├── 04_code └── 05_figures - Inside the 04_code directory, make the following directories 04_code ├── alpha_div ├── beta_div ├── taxonomy - Download the cow_metadata.txt, shannon.tsv, unweighted_unifrac.txt, tabulated_results.tsv, and cow_HW4_r.Rmd files from Canvas and put them in the correct directories. **What directory should the cow_HW4_r.Rmd file go in? - *Write the directory here:* #### Statistical analysis and figure generation in R - Now that we have set up the correct file structure and put our files in the correct directories, we can start our cow R analysis. - Open the cow_HW4_r.Rmd file and start working through the analysis. **Note that if you open the markdown file in your Downloads, the working directory will not be correct. Make sure to only open the markdown file after you have put it in the correct working directory.** **Read in metadata - Fill in the file path you used in the R Markdown to load the metadata. ``` metadata <- read_tsv("YOUR FILE PATH HERE") ``` **Read in alpha diversity data ** - Fill in the file path you used in the R Markdown to load the shannon data ``` shannon <- read_tsv("YOUR FILE PATH HERE") ``` **Read in beta diversity data ** - Fill in the file path you used in the R Markdown to load the unweighted unifrac data ``` uw_unifrac <- read_tsv("YOUR FILE PATH HERE") ``` **Load in tabulated results** - Fill in the file path you used in the R Markdown to load the tabulated_results.tsv ``` tabulated_results <- read_tsv("YOUR FILE PATH HERE") ```