Fastp Fastq Preprocessing

# Install with conda
conda install -c bioconda fastp

# Or download pre-compiled binary (Linux)
wget https://github.com/OpenGene/fastp/releases/download/v0.24.0/fastp
chmod +x fastp
./fastp --version
# fastp 0.24.0

# Verify
fastp --version

decisions:
  - id: D1
    param: -i / -I
    kind: required
    source: data
    ask: "Were these reads sequenced from single-end or paired-ends of each fragment?"
    default: "read off the file inventory - paired when R1/R2 mates are present"

  - id: D2
    param: library_prep
    kind: required
    source: user
    ask: "Which library prep produced these reads - 5'-, 3'-, small RNA, lnRNA or amplicon?"
    default: null

  - id: D3
    param: --trim_poly_x
    kind: optional
    source: user
    depends_on: [D2]
    ask: "Trim polyA/polyT tails left by the capture step?"
    default: "on for 3'-enriched RNA-seq, off otherwise"

  - id: D4
    param: --adapter_sequence / --detect_adapter_for_pe
    kind: required
    source: user
    depends_on: [D1]
    ask: "Auto-detect the adapter, or supply the kit's sequence explicitly?"
    default: "auto-detect"

  - id: D5
    param: --length_required
    kind: required
    source: user
    depends_on: [D2]
    ask: "Reads shorter than this after trimming are discarded - how short is too short to map?"
    default: "15 (raise to ~36 for standard RNA-seq; keep low for small RNA)"

  - id: D6
    param: --qualified_quality_phred
    kind: optional
    source: user
    ask: "How confident must a base call be before it counts as good quality?"
    default: "15 (20 = 1% error)"

  - id: D7
    param: --correction
    kind: optional
    source: user
    depends_on: [D1]
    ask: "Correct mismatched bases where the two mates overlap?"
    default: "off"
    skip_if: "single-end input - there is no overlap to compare"

  - id: D8
    param: --low_complexity_filter
    kind: optional_conditional
    source: data
    ask: "Drop low-complexity reads (homopolymer and short-repeat runs)?"
    default: "off"

  - id: D9
    param: --thread, --split
    kind: never_ask
    source: data
    reason: "Parallelism and output sharding; affect runtime and file layout, not the reads kept"
    default: "8 threads, no split"

# Paired-end adapter trimming with QC report
fastp \
    -i sample_R1.fastq.gz \
    -I sample_R2.fastq.gz \
    -o sample_R1.trimmed.fastq.gz \
    -O sample_R2.trimmed.fastq.gz \
    -h sample_qc.html \
    -j sample_qc.json \
    --thread 8

echo "Trimmed reads in: sample_R1.trimmed.fastq.gz"

# Single-end with auto adapter detection
fastp \
    -i sample.fastq.gz \
    -o sample.trimmed.fastq.gz \
    -h sample_qc.html \
    -j sample_qc.json \
    --thread 8 \
    --qualified_quality_phred 20 \
    --length_required 36

echo "Input reads:   $(zcat sample.fastq.gz | wc -l | awk '{print $1/4}')"
echo "Output reads:  $(zcat sample.trimmed.fastq.gz | wc -l | awk '{print $1/4}')"

# Paired-end with overlap-based adapter auto-detection
fastp \
    -i sample_R1.fastq.gz \
    -I sample_R2.fastq.gz \
    -o sample_R1.trimmed.fastq.gz \
    -O sample_R2.trimmed.fastq.gz \
    -h sample_qc.html \
    -j sample_qc.json \
    --thread 8 \
    --correction \
    --detect_adapter_for_pe \
    --qualified_quality_phred 20 \
    --length_required 36

# Specify adapters explicitly (if auto-detection fails)
# fastp -i R1.fq.gz -I R2.fq.gz \
#   --adapter_sequence AGATCGGAAGAGCACACGTCTGAACTCCAGTCA \
#   --adapter_sequence_r2 AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT \
#   -o R1.out.fq.gz -O R2.out.fq.gz

# Strict quality filtering (e.g., for variant calling)
fastp \
    -i sample_R1.fastq.gz \
    -I sample_R2.fastq.gz \
    -o sample_R1.filtered.fastq.gz \
    -O sample_R2.filtered.fastq.gz \
    -h sample_qc.html \
    -j sample_qc.json \
    --thread 8 \
    --qualified_quality_phred 25 \
    --unqualified_percent_limit 20 \
    --length_required 50 \
    --max_len1 150 \
    --max_len2 150 \
    --low_complexity_filter \
    --complexity_threshold 30

echo "Filtering complete. Check sample_qc.html for pass/fail rates."

# Remove polyA tails (QuantSeq 3′ mRNA-seq)
fastp \
    -i quantseq_R1.fastq.gz \
    -o quantseq_R1.trimmed.fastq.gz \
    -h quantseq_qc.html \
    -j quantseq_qc.json \
    --thread 8 \
    --trim_poly_x \
    --poly_x_min_len 10 \
    --qualified_quality_phred 20 \
    --length_required 25

# For Smart-seq2 paired-end with polyA
fastp \
    -i smartseq_R1.fastq.gz \
    -I smartseq_R2.fastq.gz \
    -o smartseq_R1.trimmed.fastq.gz \
    -O smartseq_R2.trimmed.fastq.gz \
    --trim_poly_x --poly_x_min_len 10 \
    --thread 8 \
    -h smartseq_qc.html -j smartseq_qc.json

import json
from pathlib import Path

def parse_fastp_json(json_path: str) -> dict:
    with open(json_path) as f:
        data = json.load(f)
    
    before = data["summary"]["before_filtering"]
    after = data["summary"]["after_filtering"]
    
    return {
        "total_reads_in":  before["total_reads"],
        "total_reads_out": after["total_reads"],
        "pct_passed":      after["total_reads"] / before["total_reads"] * 100,
        "q30_rate_before": before["q30_rate"] * 100,
        "q30_rate_after":  after["q30_rate"] * 100,
        "mean_len_before": before["read1_mean_length"],
        "mean_len_after":  after["read1_mean_length"],
        "adapter_trimmed": data["filtering_result"]["adapter_trimmed"],
    }

metrics = parse_fastp_json("sample_qc.json")
for key, val in metrics.items():
    print(f"{key:25s}: {val:.1f}" if isinstance(val, float) else f"{key:25s}: {val:,}")

# Quality gate: fail if < 70% reads pass filter
if metrics["pct_passed"] < 70:
    print("WARNING: Low pass rate — check raw data quality")

#!/bin/bash
# Batch paired-end preprocessing for multiple samples
SAMPLES=(ctrl_1 ctrl_2 treat_1 treat_2)
DATA="data"
OUT="trimmed"
QC="qc/fastp"
THREADS=8

mkdir -p "$OUT" "$QC"

for sample in "${SAMPLES[@]}"; do
    echo "=== Processing $sample ==="
    fastp \
        -i "$DATA/${sample}_R1.fastq.gz" \
        -I "$DATA/${sample}_R2.fastq.gz" \
        -o "$OUT/${sample}_R1.fastq.gz" \
        -O "$OUT/${sample}_R2.fastq.gz" \
        -h "$QC/${sample}.html" \
        -j "$QC/${sample}.json" \
        --thread $THREADS \
        --correction \
        --detect_adapter_for_pe \
        --qualified_quality_phred 20 \
        --length_required 36 \
        2>&1 | grep -E "Read[12]|Filtering|Adapter|passed"
done

# Aggregate QC metrics
python3 - << 'EOF'
import json, pandas as pd
from pathlib import Path

rows = []
for jf in sorted(Path("qc/fastp").glob("*.json")):
    with open(jf) as f: data = json.load(f)
    after = data["summary"]["after_filtering"]
    before = data["summary"]["before_filtering"]
    rows.append({
        "sample": jf.stem,
        "reads_in": before["total_reads"],
        "reads_out": after["total_reads"],
        "pct_passed": round(after["total_reads"]/before["total_reads"]*100, 1),
        "q30_after": round(after["q30_rate"]*100, 1),
    })
df = pd.DataFrame(rows)
print(df.to_string(index=False))
df.to_csv("fastp_summary.tsv", sep="\t", index=False)
EOF

# Run MultiQC to aggregate all fastp JSON reports
multiqc qc/fastp/ -o qc/ -n fastp_multiqc_report