A method to estimate RT error rates.

1 Estimate RT error rates through sequencing

I want to first figure out what is in this data. In order to do that, I am thinking that I will need to do the following:

1. Combine the read1/read2 pairs.
2. Do simple grep searches for the end of the template sequences to see how many reads have them.
3. Count the 6nt indexes.

1.1 Combine read pairs

I am using flash for this. I copied the raw data into preprocessing/s1, s2, and s3 arbitrarily.

module add flash
cd preprocessing/s1
flash <(less *_1.fastq.gz) <(less *_2.fastq.gz)
cd ../s2
flash <(less *_1.fastq.gz) <(less *_2.fastq.gz)
cd ../s3
flash <(less *_1.fastq.gz) <(less *_2.fastq.gz)

Now I have a file in each directory named out.extendedFrags.fastq in each directory.

Let us find out how many reads have the end of the template.

Reading the template in the forward direction, we expect to find: TTGTAATACGACTCAC

In the reverse direction, we expect to find: GTGAGTCGTATTACAA

fwd="TTGTAATACGACTCAC"
rev="GTGAGTCGTATTACAA"
cd preprocessing/s1
reads=$(( $(wc -l out.extendedFrags.fastq | awk '{print $1}') / 4 ))
echo "There are ${reads} reads"
## There are 4,667,506 reads
fwd_reads=$(grep ${fwd} out.extendedFrags.fastq | wc -l)
echo $fwd_reads
## 2,285,256
rev_reads=$(grep ${rev} out.extendedFrags.fastq | wc -l)
echo $rev_reads
## 2,136,744
fwd_plus_rev=$(( ${fwd_reads} + ${rev_reads} ))
missing_reads=$(( ${reads} - ${fwd_plus_rev} ))
echo ${missing_reads}
## 245,506

cd ../s2
reads=$(( $(wc -l out.extendedFrags.fastq | awk '{print $1}') / 4 ))
echo "There are ${reads} reads"
## There are 1,847,635 reads
fwd_reads=$(grep ${fwd} out.extendedFrags.fastq | wc -l)
echo $fwd_reads
## 988,931
rev_reads=$(grep ${rev} out.extendedFrags.fastq | wc -l)
echo $rev_reads
## 747,696
fwd_plus_rev=$(( ${fwd_reads} + ${rev_reads} ))
missing_reads=$(( ${reads} - ${fwd_plus_rev} ))
echo ${missing_reads}
## 111,008

cd ../s3
reads=$(( $(wc -l out.extendedFrags.fastq | awk '{print $1}') / 4 ))
echo "There are ${reads} reads"
## There are 6,150,081
fwd_reads=$(grep ${fwd} out.extendedFrags.fastq | wc -l)
echo $fwd_reads
## 2,688,628
rev_reads=$(grep ${rev} out.extendedFrags.fastq | wc -l)
echo $rev_reads
## 3,030,038
fwd_plus_rev=$(( ${fwd_reads} + ${rev_reads} ))
missing_reads=$(( ${reads} - ${fwd_plus_rev} ))
echo ${missing_reads}
##

2 Extract index containing sequences

I wrote a short perl script which attempts to pull out the sequences which contain the template. For each sequence with the template, it then writes the reads as a fasta file where the ID of each sequence is the 14 nt. index. Thus we should be able to scan the actual sequences for mismatches.

cd jdestefa_160920/s1_CGATGT && \
    ../count_indexes.pl  out.extendedFrags.fastq.gz
cd jdestefa_160920/s2_CGATGT && \
    ../count_indexes.pl  out.extendedFrags.fastq.gz
cd jdestefa_160920/s3_CGATGT && \
    ../count_indexes.pl  out.extendedFrags.fastq.gz

The above three commands also print out a file ‘idx_count.txt’ which provides a count of how many times each random index was observed. Let us plot that information.

s1 <- read.table("jdestefa_160920/s1_CGATGT/idx_count.txt", sep=" ")
colnames(s1) <- c("s1_index", "s1_observations")
s1 <- data.table::as.data.table(s1)
s2 <- read.table("jdestefa_160920/s2_ATCACG/idx_count.txt", sep=" ")
colnames(s2) <- c("s2_index", "s2_observations")
s2 <- data.table::as.data.table(s2)
s3 <- read.table("jdestefa_160920/s3_TTAGGC/idx_count.txt", sep=" ")
colnames(s3) <- c("s3_index", "s3_observations")
s3 <- data.table::as.data.table(s3)
all <- merge(s1, s2, by.x="s1_index", by.y="s2_index", all.x=TRUE, all.y=TRUE)
all <- merge(all, s3, by.x="s1_index", by.y="s3_index", all.x=TRUE, all.y=TRUE)
colnames(all)[1] <- "index"

na_idx <- is.na(all)
all[na_idx] <- 0
all <- as.data.frame(all)
rownames(all) <- all[["index"]]
all[["index"]] <- NULL

library(ggplot2)
ggplot(s1, aes(x=s1_observations)) +
  geom_histogram(binwidth=1) +
  scale_x_continuous(limits=c(0,20))

## Warning: Removed 1 rows containing non-finite values (stat_bin).

## Warning: Removed 2 rows containing missing values (geom_bar).

ggplot(s2, aes(x=s2_observations)) +
  geom_histogram(binwidth=1) +
  scale_x_continuous(limits=c(0,20))

## Warning: Removed 30 rows containing non-finite values (stat_bin).

## Warning: Removed 2 rows containing missing values (geom_bar).

ggplot(s3, aes(x=s3_observations)) +
  geom_histogram(binwidth=1) +
  scale_x_continuous(limits=c(0,20))

## Warning: Removed 7 rows containing non-finite values (stat_bin).

## Warning: Removed 2 rows containing missing values (geom_bar).

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))