---
title: "02 - Read QC and trimming"
output: html_document
---

Run FastQC on the raw reads, trim adapters and low-quality bases with Trim Galore,
re-run FastQC on the trimmed reads, and summarize everything with MultiQC.

- Input: `data/raw-reads/SRR*.fastq.gz` (from `01-download-data`), `data/sample_sheet.csv`
- Output reads: `data/trimmed-reads/` (gitignored)
- Output reports: `output/02-qc-trim/`

Requires FastQC, Trim Galore and MultiQC. Create the environment once with:

```
conda env create -f envs/qc-trim.yml
```

```{r setup, include=FALSE}
# Paths are relative to code/, so chunks work the same run one at a time or knitted
knitr::opts_chunk$set(echo = TRUE, eval = TRUE)
```

# FastQC on raw reads

```{bash fastqc-raw}
set -euo pipefail
export PATH=/mmfs1/gscratch/srlab/sr320/miniforge3/envs/qc-trim/bin:$PATH

out=../output/02-qc-trim/fastqc-raw
mkdir -p "${out}"

fastqc \
  --threads 6 \
  --outdir "${out}" \
  ../data/raw-reads/SRR*.fastq.gz
```

# Trimming

Trim Galore trims adapters and bases below Q20, dropping reads shorter than 20 bp.
The reads are single-end 101 bp.

The adapters are given explicitly instead of auto-detected. With auto-detection, Trim Galore searches
only for the first 13 bp of the Illumina adapter (`AGATCGGAAGAGC`), which missed many adapter dimers:
1.6–6.5% of trimmed reads in most samples, and 22% in HB30, were still mostly adapter.
They were TruSeq adapter with sequencing errors in the first bases (`GACCGGAAGAGC…`) or with the
first four bases missing (`CGGAAGAGC…`). Two adapters catch both:

- the full 33 bp TruSeq adapter, which tolerates a few mismatches;
- the same adapter without its first four bases, for dimers that start mid-adapter.

On 500,000 HB30 reads this cut leftover adapter from 18% to 0.03%.
`--no_poly_g` turns off poly-G trimming, which Trim Galore 2.3 switches on automatically.
Poly-G tails are an artifact of 2-colour chemistry (NextSeq/NovaSeq); these reads are from a
4-colour HiSeq 2500, so G runs are real sequence.

`clip_5` and `clip_3` remove a fixed number of bases from each end of every read. They are 0 for now.
Set them after checking the Bismark M-bias plots in step 4, then re-run this chunk and alignment.
If the library turns out to be Zymo Pico Methyl-Seq, the Trim Galore docs suggest 10 on each end.

```{bash trim}
set -euo pipefail
export PATH=/mmfs1/gscratch/srlab/sr320/miniforge3/envs/qc-trim/bin:$PATH

clip_5=0
clip_3=0

reads_out=../data/trimmed-reads
report_out=../output/02-qc-trim/trim-galore
fastqc_out=../output/02-qc-trim/fastqc-trimmed
mkdir -p "${reads_out}" "${report_out}" "${fastqc_out}"

clip_args=()
if (( clip_5 > 0 )); then clip_args+=(--clip_R1 "${clip_5}"); fi
if (( clip_3 > 0 )); then clip_args+=(--three_prime_clip_R1 "${clip_3}"); fi

for fq in ../data/raw-reads/SRR*.fastq.gz; do
  trim_galore \
    --cores 20 \
    --quality 20 \
    --length 20 \
    --no_poly_g \
    --adapter AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
    --adapter CGGAAGAGCACACGTCTGAACTCCAGTCAC \
    "${clip_args[@]}" \
    --fastqc \
    --fastqc_args "--outdir ${fastqc_out}" \
    --output_dir "${reads_out}" \
    "${fq}"
done

mv "${reads_out}"/*_trimming_report.{txt,json} "${report_out}"/
ls -lh "${reads_out}"
```

## Check for leftover adapter

Count reads in the first 500,000 of each trimmed file that still contain the adapter core
(`GAAGAGCACACGTCTGAAC`). This should be well under 1% for every sample.

```{bash adapter-check}
set -euo pipefail
export PATH=/mmfs1/gscratch/srlab/sr320/miniforge3/envs/qc-trim/bin:$PATH

echo "run,reads_checked,pct_with_adapter" > ../output/02-qc-trim/leftover-adapter.csv
for fq in ../data/trimmed-reads/SRR*_trimmed.fq.gz; do
  run=$(basename "${fq}" _trimmed.fq.gz)
  # head closes the pipe early on purpose, so ignore pigz's SIGPIPE exit
  { pigz -dc "${fq}" || true; } | head -n 2000000 \
    | awk -v run="${run}" 'NR % 4 == 2 { n++; if (index($0, "GAAGAGCACACGTCTGAAC")) a++ }
        END { printf "%s,%d,%.2f\n", run, n, 100 * a / n }' \
    >> ../output/02-qc-trim/leftover-adapter.csv
done
cat ../output/02-qc-trim/leftover-adapter.csv
```

# MultiQC

One report covering raw FastQC, trimming, and trimmed FastQC.

```{bash multiqc}
set -euo pipefail
export PATH=/mmfs1/gscratch/srlab/sr320/miniforge3/envs/qc-trim/bin:$PATH

multiqc \
  --force \
  --title "SRP139854 read QC" \
  --outdir ../output/02-qc-trim/multiqc \
  ../output/02-qc-trim
```

# Reads kept after trimming

```{bash read-counts}
set -euo pipefail
export PATH=/mmfs1/gscratch/srlab/sr320/miniforge3/envs/qc-trim/bin:$PATH

echo "run,reads_trimmed" > ../output/02-qc-trim/read-counts-trimmed.csv
for fq in ../data/trimmed-reads/SRR*_trimmed.fq.gz; do
  run=$(basename "${fq}" _trimmed.fq.gz)
  n=$(( $(pigz -dc "${fq}" | wc -l) / 4 ))
  echo "${run},${n}" >> ../output/02-qc-trim/read-counts-trimmed.csv
done
```

```{r compare-counts}
sample_sheet <- read.csv("../data/sample_sheet.csv")
raw <- read.csv("../output/01-download-data/read-counts.csv")
trimmed <- read.csv("../output/02-qc-trim/read-counts-trimmed.csv")

counts <- Reduce(function(x, y) merge(x, y, by = "run"),
                 list(sample_sheet[, c("run", "sample", "treatment")], raw, trimmed))
counts$pct_kept <- round(100 * counts$reads_trimmed / counts$reads, 2)
counts <- counts[order(counts$treatment, counts$sample), ]

stopifnot(nrow(counts) == nrow(sample_sheet))

write.csv(counts, "../output/02-qc-trim/read-counts-summary.csv", row.names = FALSE, quote = FALSE)
counts
```