<jupyter_start><jupyter_text>Dada2 R Microbiome TutorialIf the required packages are not on your system please uncomment and run the first chunk.<jupyter_code>## Download Libraries
#source("https://bioconductor.org/biocLite.R")
#biocLite("dada2")
#biocLite("phyloseq")
#install.packages("RColorBrewer")
#install.packages("ggplot2")
## Load Required Packages
library(dada2)
library(phyloseq)
library(ggplot2)
library(RColorBrewer)
## Set path to forward and reverse reads
path = "./data"
list.files(path)
fnFS = sort(list.files(path, pattern = "R1.fastq"))
fnRs = sort(list.files(path, pattern = "R2.fastq"))

samp_names = sapply(strsplit(fnFS, "_"), `[`,2)

fnFS = file.path(path, fnFS)
fnRs = file.path(path, fnRs)<jupyter_output><empty_output><jupyter_text>Create new files for filtered reads<jupyter_code>filt_path = file.path(path, "filtered")
filtFs = file.path(filt_path, paste0(samp_names, "_F_filt.fastq.gz"))
filtRs = file.path(filt_path, paste0(samp_names, "_R_filt.fastq.gz"))<jupyter_output><empty_output><jupyter_text>Filter and Trim Function* The Filter and Trim Function allows a user to truncate reads by set lengths. In this demo, we are trimming to 230 base pairs for both forward and reverse reads. * rm.phix maps and removes reads that align with phiX contamination.* maxEE (Not used) sets the maximum number of expected errors in a read. [More Information](http://www.drive5.com/usearch/manual/expected_errors.html)<jupyter_code>out = filterAndTrim(fnFS, filtFs, fnRs, filtRs, truncLen = c(230, 230),maxN = 0,truncQ = 2,rm.phix = TRUE,compress = TRUE,multithread = TRUE)
head(out)<jupyter_output><empty_output><jupyter_text>Learn Error RateDada2 depends on a parametric error model within each amplicon dataset. The error rate is determined by alternating estimates of error rate and infering sample composition until they converge. The inital estimate of an error rate is that all sequences other than the most abundant are errors.<jupyter_code>errF = learnErrors(filtFs, multithread = T)
errR = learnErrors(filtRs, multithread = T)<jupyter_output>Initializing error rates to maximum possible estimate.
Sample 1 - 22515 reads in 15098 unique sequences.
   selfConsist step 2 
   selfConsist step 3 
   selfConsist step 4 
Convergence after  4  rounds.
Total reads used:  22515<jupyter_text>DereplicateThis function bins reads into unique sequence categories to reduce the computational overhead of variant comparisons.<jupyter_code>derepFs = derepFastq(filtFs, verbose = T)
derepRs = derepFastq(filtRs, verbose = T)

names(derepFs) = samp_names
names(derepRs) = samp_names<jupyter_output>Dereplicating sequence entries in Fastq file: ./data/filtered/mock_F_filt.fastq.gz
Encountered 9642 unique sequences from 22515 total sequences read.
Dereplicating sequence entries in Fastq file: ./data/filtered/mock_R_filt.fastq.gz
Encountered 15098 unique sequences from 22515 total sequences read.
Warning message in `names<-`(`*tmp*`, value = "mock"):
“derep-class objects cannot be renamed.”Warning message in `names<-`(`*tmp*`, value = "mock"):
“derep-class objects cannot be renamed.”<jupyter_text>Run Core Denoising Algorithm and MergingThe Dada2 algorithm uses a p-value threshold that compares each unique bin to another more prevelant variant. The learned error rate in conjunction with the abundance information of the bin is used to determine if a sequence is a real sequence variant or mearly an error/point mutation.<jupyter_code>dadaFs = dada(derepFs, err = errF, multithread = T)
dadaRs = dada(derepRs, err = errR, multithread = T)
mergers = mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose = T)
seqtab = makeSequenceTable(mergers)
dim(seqtab)
## Create table of merged sequence length
table(nchar(getSequences(seqtab)))<jupyter_output><empty_output><jupyter_text>Perform chimera filtering<jupyter_code>seqtab.nochim = removeBimeraDenovo(seqtab, method = "consensus", multithread=TRUE, verbose = T)
(sum(seqtab.nochim)/sum(seqtab))<jupyter_output>Identified 4 bimeras out of 18 input sequences.<jupyter_text>Follow read counts through pipeline<jupyter_code>getN <- function(x) sum(getUniques(x))
track <- cbind(out, getN(dadaFs), getN(mergers), rowSums(seqtab), rowSums(seqtab.nochim))
colnames(track) <- c("input", "filtered", "denoised", "merged", "tabled", "nonchim")
rownames(track) <- samp_names
track

## Use if more than one sample to create table
#getN <- function(x) sum(getUniques(x))
#track <- cbind(out, sapply(dadaFs, getN), sapply(mergers, getN), rowSums(seqtab), rowSums(seqtab.nochim))
#colnames(track) <- c("input", "filtered", "denoised", "merged", "tabled", "nonchim")
#rownames(track) <- samp_names

## Save file as csv
#write.csv(track, file = "data/Read_distro_thru_Processing.csv")<jupyter_output><empty_output><jupyter_text>Assign Taxonomy to VariantsThis chunk is comprised of two components: The first assigns variants taxanomic information from the Silva database to the genus level. The second expands upon the previous assignment to add the species level if distinguishable.<jupyter_code>taxa_table = assignTaxonomy(seqtab.nochim, "./data/database/silva_nr_v128_train_set.fa.gz", multithread = T)

## Add Species Information
tax_silva_species = addSpecies(taxa_table, "./data/database/silva_species_assignment_v128.fa.gz")<jupyter_output><empty_output><jupyter_text>Load Count Matrix and Taxonomy into Phyloseq Object<jupyter_code>ps = phyloseq(otu_table(seqtab.nochim,taxa_are_rows = FALSE), tax_table(tax_silva_species))
colnames(tax_table(ps))<jupyter_output><empty_output><jupyter_text>Count Number<jupyter_code>otu_glom<-tax_glom(ps, taxrank="Genus")

Percentage <- transform_sample_counts(ps, function(x) {(x/sum(x))})

percent = as.data.frame(t(Percentage@otu_table@.Data))
tax = as.data.frame(Percentage@tax_table)

percent$sa1 = percent$sa1*100
new_table = as.data.frame(tax$Genus)
new_table$percent = percent$sa1
colnames(new_table) = c("Genus", "Read Percentage")

plot_bar(ps, fill = "Genus")
new_table
tab = as.data.frame(table(tax$Genus))
colnames(tab) = c("Genus", "Frequency")
tab<jupyter_output><empty_output><jupyter_text>Merge Genus Level Counts<jupyter_code>### Genus Level
b1 = tax_glom(ps, taxrank = "Genus")

tax_table = as.data.frame(ps@tax_table)
genus_set = as.character(unique(tax_table$Genus))
genus_set = genus_set[!is.na(genus_set)]


Genus = sort(taxa_sums(b1), TRUE)[1:7]
Genus = as.data.frame(Genus)

rownames(Genus) = genus_set
colnames(Genus) = c("Read Counts")
Genus<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11260012322858903, 'fraction_numerical': 0.012630930375847197, 'mean_word_length': 7.00616522811344, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 46, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '648506', 'license': 'MIT', 'n_tokens_mistral': 2362, 'n_tokens_neox': 2213, 'n_words': 579, 'path': 'Dada2_R.ipynb', 'repo_name': 'glickmac/Microbiome_Demo', 'size': '6492'}