LRII Bioinformatics

/project/microbiome/data_queue/seq/LowReadII/rawdata

salloc --account=microbiome -t 0-06:00

mkdir -p /gscratch/grandol1/LowReadII/rawdata

cd /gscratch/grandol1/LowReadII/rawdata

unpigz --to-stdout /project/microbiome/data_queue/seq/LowReadII/rawdata/Low-Read-II_S1_L001_R1_001.fastq | split -l 1000000 -d --suffix-length=3 --additional-suffix=.fastq - LowReadII_R1_ ;unpigz --to-stdout /project/microbiome/data_queue/seq/LowReadII/rawdata/Low-Read-II_S1_L001_R2_001.fastq | split -l 1000000 -d --suffix-length=3 --additional-suffix=.fastq - LowReadII_R2_

//project/microbiome/data_queue/seq/LowReadII/rawdata/run_parse_count_onSplitInput.pl

cd /project/microbiome/data_queue/seq/LowReadII/rawdata

./run_splitFastq_fwd.sh

./run_splitFastq_rev.sh

./run_aggregate.sh

cd /project/microbiome/data_queue/seq/LowReadII/tfmergedreads

./run_slurm_mergereads.pl

cd /project/microbiome/data_queue/seq/LowReadII/otu

./run_slurm_mkotu.pl

salloc --account=microbiome -t 0-02:00 --mem=500G

module load swset/2018.05  gcc/7.3.0

module load vsearch/2.15.1

vsearch --sintax zotus.fa --db /project/microbiome/users/grandol1/ref_db/gg_16s_13.5.fa -tabbedout LRII.sintax -sintax_cutoff 0.8

awk -F "\t" '{OFS=","} NR==1 {print "OTU_ID","SEQS","SIZE","DOMAIN","KINGDOM","PHYLUM","CLASS","ORDER","FAMILY","GENUS","SPECIES"} {gsub(";", ","); gsub("centroid=", ""); gsub("seqs=", ""); gsub("size=", ""); match($4, /d:[^,]+/, d); match($4, /k:[^,]+/, k); match($4, /p:[^,]+/, p); match($4, /c:[^,]+/, c); match($4, /o:[^,]+/, o); match($4, /f:[^,]+/, f); match($4, /g:[^,]+/, g); match($4, /s:[^,]+/, s); print $1, d[0]=="" ? "NA" : d[0], k[0]=="" ? "NA" : k[0], p[0]=="" ? "NA" : p[0], c[0]=="" ? "NA" : c[0], o[0]=="" ? "NA" : o[0], f[0]=="" ? "NA" : f[0], g[0]=="" ? "NA" : g[0], s[0]=="" ? "NA" : s[0] }' LRII.sintax > LRIItaxonomy.csv

library(tidyverse)
library(ggplot2)

#read in otu table
LowReadOtu <- read.csv("/Volumes/Macintosh HD/Users/gregg/Desktop/WorkingFolder/otutable.txt", sep = "", header = TRUE, stringsAsFactors = FALSE)

#LowReadReads <- read.csv("/Volumes/Macintosh HD/Users/gregg/Desktop/WorkingFolder/filtermergestats.csv", header=FALSE, stringsAsFactors = FALSE, na.strings = "")

#LowReadReads <- separate(LowReadReads, 1, sep="[.]", into = c(NA, "samplename", NA, NA, NA, NA, NA, NA, NA))

#LowReadOtu <- as.data.frame(LowReadOtu)

#transpose
LowReadOtu<- LowReadOtu %>% 
  rownames_to_column() %>% 
  gather(variable, value, -rowname) %>% 
  spread(rowname, value)
colnames(LowReadOtu) <- LowReadOtu[1,]
LowReadOtu <- LowReadOtu[-1, ] 

#simplify naming
LowReadOtu$OTUID <- gsub("rna16S_[A-Z]+_[A-Z]+_","", LowReadOtu[,1])


#Sort by name
LowReadOtu <- arrange(LowReadOtu, OTUID)


#simplify otu names
colnames(LowReadOtu) <- gsub("centroid=", "", colnames(LowReadOtu))
LowReadOtu <- type_convert(LowReadOtu)
LowReadOtu[is.na(LowReadOtu)] <- 0

substrRight <- function(x, n){
  substr(x, nchar(x)-n+1, nchar(x))
}
LowReadOtu$Normalized <- gsub("[0-9]", "", substrRight(LowReadOtu$OTUID, 1))
LowReadOtu <- arrange(LowReadOtu, Normalized)
LowReadOtu$Volume <- sub('.*(\\d{1}).*', '\\1', LowReadOtu$OTUID)

#Total reads assigned to otus
LowReadOtu <- LowReadOtu %>%
  rowwise()%>%
  mutate(reads = sum(c_across(otu1:otu1967)))

#make sums first observation
LowReadOtu <- cbind(LowReadOtu[1], LowReadOtu[,3803], LowReadOtu[2:3802])
#Calculate Percent ISD for all samples
LowReadOtu$PercISD <- LowReadOtu$otu1 / LowReadOtu$reads
LowReadOtu <- LowReadOtu[,c(1, 3804, 2:3803)]

#add repeats together maintaining V and N categorizations
LRIIotu <- LowReadOtu %>%
  group_by(OTUID, Volume, Normalized) %>%
  summarise(TotalReads = sum(reads))
ggplot(LRIIotu, aes(TotalReads, Volume, color=OTUID, shape=Normalized))+
  geom_point(show.legend = FALSE)

#fix odd naming convention SAG_ vs SAG#
LRIIotu$OTUID <- gsub("SAG_", "SAG", LRIIotu$OTUID)
#remove V and N appendices from samplenames
LRIIotu$OTUID <- gsub("_.*$", "", LRIIotu$OTUID)
ggplot(LRIIotu[45:72,], aes(TotalReads, Volume, color=OTUID, shape=Normalized))+
  geom_point()+
  facet_wrap(~OTUID)
ggplot(LRIIotu[1:44,], aes(TotalReads, Volume, color=OTUID, shape=Normalized))+
  geom_point()+
  facet_wrap(~OTUID)
ggplot(LRIIotu, aes(TotalReads, Volume, color=OTUID, shape=Normalized))+
  geom_point()+
  facet_wrap(~OTUID)

#Explore how the blanks performed
BlanksOnly <- LowReadOtu[c(1:8, 73:80),] 
library(janitor)
BlanksOnly <- BlanksOnly %>%
  adorn_totals("row")
BlanksOnly <- BlanksOnly[which(BlanksOnly[17,] != 0)]
#Some contamination, but mostly ISD.  One replicate of Blank 1_2N is suspect at only 67% ISD.
#How many samples are more than 65% ISD?
LowReadOtu$OTUID[which(LowReadOtu$PercISD > 0.65)]

#Read in taxonomy data
LowReadSyntax <- read.csv("/Volumes/Macintosh HD/Users/gregg/Desktop/WorkingFolder/LRIItaxonomy.csv", header = TRUE, stringsAsFactors = FALSE)

#Simplify table
LowReadSyntax <- LowReadSyntax[,c(1, 3:11) ]


#Is Mock Community largely List mono, Pseu aer, Bac sub, E coli, Sal ent, Lac fer, Ent fae, and Sta aur ??? Isolate MC
MConly <- LowReadOtu[c(9:16, 81:88),]
#add totals per otu
MConly <- MConly %>%
  adorn_totals("row")
#remove empty otu columns
MConly <- MConly[which(MConly[17,] != 0)]

#isolate taxonomy data for MC samples only
MConlyTaxa <- LowReadSyntax[LowReadSyntax$OTU_ID %in% colnames(MConly)[4:65],]

#on a cursory examination, Family seems to be the most consistently called lowest taxonomic level for these samples at 2x151 sequencing.
MCtaxa <- unique(MConlyTaxa$FAMILY)
MCtaxa

#Mock community bacterial community should be:
#Listeria monocytogenes - 12%, Pseudomonas aeruginosa - 12%, Bacillus subtilis - 12%, Escherichia coli - 12%, Salmonella enterica - 12%, Lactobacillus fermentum - 12%, Enterococcus faecalis - 12%, Staphylococcus aureus - 12%, Saccharomyces cerevisiae - 2%, and Cryptococcus neoformans - 2%. 
#On the family level this is columns 3 through 9.  Salmonella seems to be the least resolved on this iSeq run
#3 through 9 are the expected families
MCtaxa <- MCtaxa[3:9]

#isolate the contaminants
MCcontaminantsTaxa <- subset(MConlyTaxa, !(FAMILY %in% MCtaxa))

#otu1 is not a contaminant; it is ISD.  Isolate the others.
MCcontaminantsTaxa <- MCcontaminantsTaxa[2:9,]

#isolate the relevant contaminant otus
MCcon <- MCcontaminantsTaxa$OTU_ID
MCcon

#Are these contaminants in the blanks?  4 of them are
McBlankContaminants <- BlanksOnly[, c("OTUID", "reads", "otu31", "otu18", "otu51", "otu2540")]

#Look at only Samples
Samples <- LowReadOtu[c(17:72, 89:144),]
Samples <- Samples %>%
  adorn_totals("row")
#Remove empty columns
Samples <- Samples[which(Samples[113,] != 0)]

#Cursorily, these are not mostly ISD or MC
#Lets remove Blank contaminant otus and see
SamplesnoBlank <- Samples %>%
  select(-one_of(colnames(BlanksOnly)[5:65]))
#That removes 61 otus.  How does it impact read counts? Redo calculations
SamplesnoBlank <- SamplesnoBlank %>%
  rowwise()%>%
  mutate(NewReads = sum(c_across(otu1:otu1967)))
SamplesnoBlank <- SamplesnoBlank[,c(1, 3696, 2:3695)] 
#Less than 5% change in read totals removing otus contaminant in blanks
SamplesnoBlank[113,2]/SamplesnoBlank[113,4]

#How about per sample?  
SamplesnoBlank$PercDiff <- 1 - (SamplesnoBlank$NewReads/SamplesnoBlank$reads)
SamplesnoBlank <- SamplesnoBlank[, c(1, 3697, 2:3696)]
SamplesnoBlank <- arrange(SamplesnoBlank, PercDiff)
#There are actually 10 samples that are more than ~10% the contaminants from the blanks with the largest at 33%
PoorSamples <- SamplesnoBlank[104:113,]
PoorSamples <- arrange(PoorSamples, OTUID)
PoorSamples[,1:6]