Stepping through the original stamp protocol.

1 Stepping through the original STAMP protocol

In this document I am going to rerun the stamp process as close to exactly as possible to the original. I will put in this the actual commands I ran and then copy/paste the R as provided by the authors with modifications required to make it work.

2 Raw Read Trimming

The trimming process is via cutadapt, the original authors used reaper, but that is just a wrapper for cutadapt. Here is an example of my trimming:

The last adapter is the most important for this. The others may likely be just removed.

mkdir -p outputs/cutadapt && \
  less r1.fastq.xz | cutadapt -  \
    -a ACAGGTTGGATGATAAGTCCCCGGTCTGACACATCTCCCTAT -a AGATCGGAAGAGCACACGTCTGAAC -b AGATCGGAAGAGCACACGTCTGAAC  -a GTCATAGCTGTTT \
    -e 0.1 -n 3 \
    -m 8 -M 42 \
    --too-short-output=outputs/cutadapt/r1_tooshort.fastq \
    --too-long-output=outputs/cutadapt/r1_toolong.fastq \
    --untrimmed-output=outputs/cutadapt/r1_untrimmed.fastq \
     2>outputs/cutadapt.err 1>r1-trimmed_ca.fastq

## r1.fastq.xz: No such file or directory

3 Qiime usage to count reads/index

This step baffles me a little. The paper states that qiime was used with a 99% identity filter. However, after filtering the remaining index reads are only 9<=x<=13 nucleotides. It seems to me a simpler and more precise method to just take the first 9 nt of each index read and count the reads which are identical. With that in mind, here is an implementation of their qiime workflow:

## qiime needs a fasta input, here is a quick bash function to do it from:
## https://www.ecseq.com/support/ngs-snippets/convert-fastq-to-fasta-on-the-command-line
fq2fa() {
    paste - - - - < ${1} | cut -f 1,2 | sed 's/^@/>/' | tr "\t" "\n" > ${2}
}

fq2fa r1-trimmed_ca.fastq a1.fasta
## Repeat that for each sample, I dumped all the fasta files into one directory
for i in *.fasta; do
    pick_otus.py -i ${i} -o qiime
done
## Then I made a directory 'txt' into which I moved all the qiime text files.

## Error in running command bash

4 Begin their R scripts

4.1 Tally the qiime results

This is from Qiime_Tally_Data_STAMP_CAUTI.R. I am making changes to it to make it a bit less offensive to my eyes.

This is such an abomination. Also, I don’t think it can work on any R which is newer than ~ 5 years old because of the way it is calling read.table(). Its intention is to read in the .txt files produced by pick_otus.py and output a matrix with two columns per sample, the first of which contains the sequence of the index followed by the number of OTUs detected by qiime.

4.1.1 Original

## This R script will take the output files containing representative STAMP barcode sequences and
## their individual abundance in the sequencer output from all samples  that
## you sequenced and copies them in a single .csv file in the format that is
## expected from the Nb estimation script

txt_dir <- "txt"
output_file <- "qiime_tally.csv"

files <- list.files(txt_dir)
names <- files

v.length <- rep(NA, length(names))

for(f in 1:length(files)) {
  file <- file.path(txt_dir, files[f])
  mtrx <- read.table(file, sep="\t", header=FALSE, fill=TRUE)
  colnames(mtrx) <- c(paste0(files[f], "_seq"), files[f])

  ## eval(parse(text=paste("Data=as.matrix(read.table('",file,"', fill=TRUE, sep='\t'))",sep="")))
  ## colnames(Data)=c(paste(names[f],"_seq",sep=""),names[f])
  ## name=paste("Matrix",names[f],sep="")
  ## assign(name,Data)
  ## name=paste("Matrix",names[f],sep="")
  ## assign(name,Data)
  ## v.length[f]=nrow(Data)
}
max <- max(v.length) + 1
fill.up <- max - v.length

FillUpMatrix <- matrix(data=NA, nrow=fill.up[1], ncol=2)

## Error in matrix(data = NA, nrow = fill.up[1], ncol = 2): invalid 'nrow' value (too large or NA)

eval(parse(text=paste("OutputMatrix=rbind(Matrix",names[1],",FillUpMatrix)",sep="")))

## Error in eval(quote(list(...)), env): object 'Matrixa1_otus.txt' not found

for(ff in 2:length(files)) {
  eval(parse(text=paste("MatrixX=Matrix", names[ff], sep="")))
  FillUpMatrix <- matrix(data=NA, nrow=fill.up[ff], ncol=2)
  Matrix <- rbind(MatrixX, FillUpMatrix)
  OutputMarix <- cbind(OutputMatrix,Matrix)
}

## Error in eval(parse(text = paste("MatrixX=Matrix", names[ff], sep = ""))): object 'Matrixa2_otus.txt' not found

write.csv(OutputMatrix, file=outputfilename)

## Error in is.data.frame(x): object 'OutputMatrix' not found

4.1.2 Reimplemented

One caveat, I am using a more recent qiime than the authors, which provides the output of pick_otus.py in a somewhat different format (either that or they used a different mode, though I tested out 3 of them and did not get the input format they expected from any of them). Thus my names are not going to be the index sequences, but instead just ‘denovo0’ through ‘denovoNNN’.

txt_dir <- "txt"
output_file <- "qiime_tally.csv"
files <- list.files(txt_dir)
## These two lines recreate the column names used by the authors.
count_names <- gsub(x=files, pattern="_otus", replacement="")
seq_names <- paste0(count_names, "_seq")
final_df <- data.frame()
for (f in 1:length(files)) {
  filename <- files[f]
  message("Starting to read: ", filename, ".")
  count_name <- gsub(x=filename, pattern="_otus", replacement="")
  sequence_name <- paste0(count_name, "_seq")
  column_names <- c(sequence_name, count_name)
  input <- file.path(txt_dir, filename)
  lines <- readLines(input)
  pair <- as.data.frame(matrix(nrow=length(lines), ncol=2))
  colnames(pair) <- column_names
  for (l in 1:length(lines)) {
    line <- lines[l]
    split <- strsplit(x=line, split="\t")[[1]]
    pair[l, 1] <- split[1]
    pair[l, 2] <- length(split) - 1
  }
  rownames(pair) <- pair[[1]]
  if (ncol(final_df) == 0) {
    final_df <- pair
  } else {
    final_df <- merge(final_df, pair, by="row.names", all=TRUE)
    rownames(final_df) <- final_df[["Row.names"]]
    final_df[["Row.names"]] <- NULL
  }
}

## Starting to read: a1_otus.txt.

## Starting to read: a2_otus.txt.

## Starting to read: a3_otus.txt.

## Starting to read: a4_otus.txt.

## Starting to read: a5_otus.txt.

## Starting to read: a6_otus.txt.

## Starting to read: a7_otus.txt.

## Starting to read: a8_otus.txt.

## Starting to read: b2_otus.txt.

## Starting to read: b3_otus.txt.

## Starting to read: b4_otus.txt.

## Starting to read: b5_otus.txt.

## Starting to read: b6_otus.txt.

## Starting to read: b7_otus.txt.

## Starting to read: b8_otus.txt.

## Starting to read: c1_otus.txt.

## Starting to read: c2_otus.txt.

## Starting to read: c3_otus.txt.

## Starting to read: c4_otus.txt.

## Starting to read: c5_otus.txt.

## Starting to read: c6_otus.txt.

## Starting to read: c7_otus.txt.

## Starting to read: c8_otus.txt.

The above block produces a matrix in the same format as the qiime_tally_data… script. The following block will provide a similar matrix using a different method: I chose to take the first 9 nucleotides of each barcode and write out the number of times each sequence occurs with a simple perl script. Here is an example loop which iterates over every sample directory. Note that my trimmer produces xz compressed output files, I use a less preprocessor to read all of my inputs.

cd preprocessing
for d in $(find . -type d); do
    cd ${d}
    ../count_idx.pl r1-trimmed_ca.fastq.xz
done

## Error in running command bash

The above script produces text files named ‘idx_count.txt’ which I compressed. The next block reads these and creates a matrix similar to what is above with a couple of interesting differences. My intention is to dump these blocks into a small R package to make processing this data easier. Thus, it is using a sample sheet in a format similar to what we use for other data types in our lab, and processing it in a similar fashion.

setwd("/mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_stamp_2020")
meta <- extract_metadata(metadata="sample_sheets/all_samples.xlsx")

sample_df <- data.frame()
all_counts <- data.table::data.table()
cutoff <- 3
meta[["observed"]] <- 0
meta[["passed"]] <- 0
meta[["rejected"]] <- 0
meta[["max"]] <- 0
for (s in rownames(meta)) {
  table <- meta[s, "indextable"]
  sampleid <- meta[s, "sampleid"]
  element <- read.csv(file=table, sep="\t", header=FALSE)
  colnames(element) <- c("idx", "number")
  ## This tmp allocation is a little dumb, but for the moment
  ## I will use it to create a full matrix of barcodes.
  tmp <- data.table::as.data.table(element)
  colnames(tmp) <- c("idx", sampleid)
  if (ncol(all_counts) == 0) {
    all_counts <- tmp
  } else {
    all_counts <- merge(all_counts, tmp, by="idx", all=TRUE)
  }
  element[["dilution"]] <- meta[s, "dilution"]
  element[["sample"]] <- meta[s, "condition"]
  element[["sampleid"]] <- sampleid
  element[["idx"]] <- NULL
  meta[s, "observed"] <- nrow(element)
  meta[s, "reads"] <- sum(element[["number"]])
  meta[s, "max"] <- max(element[["number"]])
  keep_idx <- element[["number"]] >= cutoff
  meta[s, "passed"] <- sum(keep_idx)
  meta[s, "rejected"] <- sum(!keep_idx)
  sample_df <- rbind(sample_df, element)
}

write.csv(file="my_tag_counts.csv", x=all_counts)

all_mtrx <- as.data.frame(all_counts)
rownames(all_mtrx) <- all_mtrx[["idx"]]
all_mtrx[["idx"]] <- NULL
na_idx <- is.na(all_mtrx)
all_mtrx[na_idx] <- 0
all_mtrx <- edgeR::cpm(all_mtrx)
## Here is an aggressive possibility for filtering
null_rows <- rowSums(all_mtrx) < ncol(all_mtrx)
aggressive <- all_mtrx[!null_rows, ]
## Here is an arbitrary possibility for filtering
null_rows <- rowSums(all_mtrx) < 1
arbitrary <- all_mtrx[!null_rows, ]

meta[["aggressive"]] <- 0
meta[["arbitrary"]] <- 0
for (col in colnames(all_mtrx)) {
  agg_idx <- aggressive[, col] > 0
  meta[col, "aggressive"] <- sum(agg_idx)
  arb_idx <- arbitrary[, col] > 0
  meta[col, "arbitrary"] <- sum(arb_idx)
}


library(ggplot2)
ggplot(data=sample_df, mapping=aes_string(fill="sampleid")) +
  geom_violin(mapping=aes_string(x="sampleid", y="number")) +
  scale_y_log10()

observed <- ggplot(data=meta,
                   mapping=aes_string(x="dilution", y="observed",
                                      fill="condition", colour="condition")) +
  scale_shape_manual(values=21) +
  geom_point(size=3) +
  scale_x_log10()
observed

passed <- ggplot(data=meta,
                 mapping=aes_string(x="dilution", y="passed",
                                    fill="condition", colour="condition")) +
  scale_shape_manual(values=21) +
  geom_point(size=3) +
  scale_x_log10()
passed

arb <- ggplot(data=meta,
                 mapping=aes_string(x="dilution", y="arbitrary",
                                    fill="condition", colour="condition")) +
  scale_shape_manual(values=21) +
  geom_point(size=3) +
  scale_x_log10() +
  geom_smooth(method="lm")
arb

## `geom_smooth()` using formula 'y ~ x'

5 Describe this relationship

In the following block, I am feeding the tag and dilution data by sample to glm and looking at the result.

meta[["log"]] <- log10(meta[["dilution"]])
lm_result <- glm(formula=log ~ arbitrary, data=meta)
lm_result

## 
## Call:  glm(formula = log ~ arbitrary, data = meta)
## 
## Coefficients:
## (Intercept)    arbitrary  
##    0.731329     0.000832  
## 
## Degrees of Freedom: 22 Total (i.e. Null);  21 Residual
## Null Deviance:       113 
## Residual Deviance: 33.1  AIC: 79.6

summary(lm_result)

## 
## Call:
## glm(formula = log ~ arbitrary, data = meta)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.880  -0.535  -0.128   0.290   4.391  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.731329   0.608720    1.20     0.24    
## arbitrary   0.000832   0.000117    7.13  4.9e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 1.575)
## 
##     Null deviance: 113.217  on 22  degrees of freedom
## Residual deviance:  33.071  on 21  degrees of freedom
## AIC: 79.62
## 
## Number of Fisher Scoring iterations: 2

coef(lm_result)

## (Intercept)   arbitrary 
##   0.7313287   0.0008322

6 Predict new bacterial loads from made-up tags

Now, let us pretend that we performed the experiment on two actual samples, and in them we got 1500 and 3000 tags respectively. How many bacteria do we think there were in each sample?

test_data <- data.frame(row.names=c("t1", "t2"))
test_data[["arbitrary"]] <- c(1500, 3000)
num_found <- predict.glm(object=lm_result, newdata=test_data, type="response")

10 ^ num_found

##      t1      t2 
##   95.43 1690.67

pander::pander(sessionInfo())
message(paste0("This is hpgltools commit: ", get_git_commit()))
this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to ", this_save))
tmp <- sm(saveme(filename=this_save))