---
title: "Bacterial community analysis;  Gros et al. (doi 10.5194/bg-2022-150)"
---

This markdown describes the processing of 16S rRNA amplicon sequences derived from seawater samples along the TRANSSIZ expedition, published in https://doi.org/10.5194/bg-2022-150. The corresponding raw fastq files are available under ENA accession PRJEB50492. 

First we remove primers using *Cutadapt*  

```{console}

# copy fastq files into directory of choice
# Here adjusted to AWI server; using cutadapt bash-script

module load bio/cutadapt/3.2
bash ./../software/cutadapt.sh ./Original GTGYCAGCMGCCGCGGTAA CCGYCAATTYMTTTRAGTTT

# test rename
cd Clipped
for i in *fastq.gz; do nname=`echo $i | awk '{gsub(/bacV4V5_S[0-9]{1,2}_L001/,"clip");print}'`; echo -e $i $nname; done 

# if looking OK - execute:  
for i in *fastq.gz; do nname=`echo $i | awk '{gsub(/bacV4V5_S[0-9]{1,2}_L001/,"clip");print}'`; mv $i $nname; done

# write sampleNames for dada
ls -1 *R1_001.fastq.gz | sed 's/_R1_001\.fastq.gz//' > ../sampleNames.txt

```

*Now to DADA!*

```{r, eval=F}

# Done in RStudio within AWI-VM
# Provided IP address opened in browser
# Adjust paths etc. for your own system 

require(dada2)
require(ShortRead)
require(ggplot2)
require(gridExtra)

##########################################

# setwd 
setwd("/isibhv/projects/FRAMdata/FRAM_MicrObs/WaterCol_Autofim_Transsiz/")

# list files
path <- "/isibhv/projects/FRAMdata/FRAM_MicrObs/WaterCol_Autofim_Transsiz/Clipped"
fns <- list.files(path)
fns

# ensure fwd/rev reads  in same order
fnFs <- sort(list.files(path, pattern="_R1_001.fastq.gz"))
fnRs <- sort(list.files(path, pattern="_R2_001.fastq.gz"))

# Define sample names
sampleNames <- sort(read.table(
  "sampleNames.txt", 
  h=F, stringsAsFactors=F)$V1)

# Specify the full path to the fnFs and fnRs
fnFs <- file.path(path, fnFs)
fnRs <- file.path(path, fnRs)

#################################

# Quality check
QualityProfileFs <- list()
for(i in 1:length(fnFs)) {QualityProfileFs[[i]] <- list()
  QualityProfileFs[[i]][[1]] <- plotQualityProfile(fnFs[i])}
pdf("QualityProfileForward.pdf")
for(i in 1:length(fnFs)) {do.call("grid.arrange", 
    QualityProfileFs[[i]])}
dev.off()
rm(QualityProfileFs)

QualityProfileRs <- list()
for(i in 1:length(fnRs)) {
  QualityProfileRs[[i]] <- list()
  QualityProfileRs[[i]][[1]] <- plotQualityProfile(
    fnRs[i])}
pdf("QualityProfileReverse.pdf")
for(i in 1:length(fnRs)) {do.call("grid.arrange", 
  QualityProfileRs[[i]])}
dev.off()
rm(QualityProfileRs)

# Prepare for fastq filtering
filt_path <- file.path(path, "../Filtered")
if(!file_test("-d", filt_path)) dir.create(filt_path)
filtFs <- file.path(
  filt_path, paste0(sampleNames, "_F_filt.fastq"))
filtRs <- file.path(
  filt_path, paste0(sampleNames, "_R_filt.fastq"))

#################################

# Filter depending on expected overlap
# truncLen lowered based on QualityProfile (low quality esp. of rev-reads)
out <- filterAndTrim(
  fnFs, 
  filtFs, 
  fnRs, 
  filtRs,
  truncLen = c(230, 195),
  maxN = 0,
  minQ = 2,
  maxEE = c(3, 3), 
  truncQ = 0, 
  rm.phix = T,
  compress = F,
  multithread = 6)

head(out)
summary(out[, 2]/out[, 1])
# should be retaining >70%
# 0.8 here -- OK!

#################################

# Quality check 
QualityProfileFs.filt <- list()
for(i in 1:length(filtFs)) {
  QualityProfileFs.filt[[i]] <- list()
  QualityProfileFs.filt[[i]][[1]] <- plotQualityProfile(
    filtFs[i])}
pdf("QualityProfileForwardFiltered.pdf")
for(i in 1:length(filtFs)) {do.call("grid.arrange", 
    QualityProfileFs.filt[[i]])}
dev.off()
rm(QualityProfileFs.filt)

QualityProfileRs.filt <- list()
for(i in 1:length(filtRs)) {
  QualityProfileRs.filt[[i]] <- list()
  QualityProfileRs.filt[[i]][[1]] <- plotQualityProfile(
    filtRs[i])}
pdf("QualityProfileReverseFiltered.pdf")
for(i in 1:length(filtRs)) {  do.call("grid.arrange", 
    QualityProfileRs.filt[[i]])}
dev.off()
rm(QualityProfileRs.filt)

#################################

# Learn errors 
errF <- learnErrors(
  filtFs, multithread=6, 
  randomize=T, verbose=1, MAX_CONSIST=20)
errR <- learnErrors(
  filtRs, multithread=6, 
  randomize=T, verbose=1, MAX_CONSIST=20)

# Plot error profiles
pdf("ErrorProfiles.pdf")
plotErrors(errF, nominalQ = T)
plotErrors(errR, nominalQ = T)
dev.off()
# convergence after 5/6 rounds - ok!
# few outliers outside black line - ok!

# Dereplication 
derepFs <- derepFastq(filtFs, verbose=T)
derepRs <- derepFastq(filtRs, verbose=T)

# Rename by sampleNames
names(derepFs) <- sampleNames
names(derepRs) <- sampleNames

# Denoising
dadaFs <- dada(
  derepFs, err=errF, multithread=6, pool=T)
dadaRs <- dada(
  derepRs, err=errR, multithread=6, pool=T)

#################################

# Read merging
mergers <- mergePairs(
  dadaFs, 
  derepFs, 
  dadaRs,
  derepRs,
  minOverlap=10,
  verbose=T,
  propagateCol = c(
    "birth_fold", 
    "birth_ham"))

# create sequence table
seqtab <- makeSequenceTable(mergers)
dim(seqtab) # identified 10938 sequences
saveRDS(seqtab, "seqtab.rds")

# Removing chimeras 
# 7405 bimeras out of 10938
seqtab.nochim <- removeBimeraDenovo(
  seqtab, method="consensus", 
  multithread=12, verbose=T)

# 34 samples -- 10938 sequences
# approx. 90% kept - OK!
dim(seqtab.nochim)  
summary(rowSums(seqtab.nochim)/rowSums(seqtab))

# Determine amplicon size range 
table(rep(nchar(colnames(seqtab.nochim)), 
          colSums(seqtab.nochim)))

# Remove singletons and junk sequences
# "c" adjusted to size range of amplicons
seqtab.nochim2 <- seqtab.nochim[, nchar(
  colnames(seqtab.nochim)) %in% c(354:410) & 
    colSums(seqtab.nochim) > 1]

# Stats
dim(seqtab.nochim2) # 3125 sequences
summary(rowSums(seqtab.nochim2))
summary(rowSums(seqtab.nochim2)/rowSums(seqtab.nochim))

###################################################################################

## TAXONOMY ##

tax <- assignTaxonomy(
  seqtab.nochim2, 
  "../tax_db/silva_nr_v138_train_set.fa.gz", 
  tryRC = TRUE,
  multithread = 12)

#  Bacteria 3034 -- Archaea 61 -- Eukaryota 13 -- NA 17
summary(tax)

# select Bac/Arch
# Remove chloroplasts & mitochondria
# Remove NA on phylum level
table(tax[, 1])   
sum(is.na(tax[, 2]))   # result: 91
good <- tax[!is.na(tax[, 2]) & tax[, 1] %in% c(
  "Bacteria", "Archaea"),]
tax.good <- good[-c(
 grep("Chloroplast", good[, 4]), 
 grep("Mitochondria", good[, 5])), ]
seqtab.nochim2.good <- seqtab.nochim2[, rownames(tax.good)]
summary(rowSums(seqtab.nochim2.good))

# Format tables
seqtab.nochim2.print <- t(seqtab.nochim2.good)
tax.print <- tax.good
all.equal(rownames(seqtab.nochim2.print), 
          rownames(tax.print)) #TRUE
rownames(seqtab.nochim2.print) <- paste(
  "asv", 1:ncol(seqtab.nochim2.good), sep = "")
rownames(tax.print) <- rownames(seqtab.nochim2.print)

# Export
write.table(
  seqtab.nochim2.print,"seqtab.txt", 
  sep = "\t", quote=F)
write.table(
  tax.print,"tax.txt", 
  sep="\t", quote=F)
uniquesToFasta(
  seqtab.nochim2.good,
  "uniques.fasta")

# summary stats
getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(
  dadaFs, getN), sapply(mergers, getN), 
  rowSums(seqtab), rowSums(seqtab.nochim2))
colnames(track) <- c(
  "input","filtered","denoised",
  "merged","nochim","tabled")
rownames(track) <- sampleNames
track <- data.frame(track)
head(track)

write.table(track, 
  "dadastats.txt", 
  quote=F, sep="\t")

```

The resulting files - seqtab.txt and tax.txt - are then processed in context of environmental parameters using Rscripts 1_DataLoad.R, 2_Results.R and 3_Supplement.R. All necessary input files to reproduce analyses & figures are provided.