Genomics Pipelines

Genomics Pipelines

Identity

Reference System Usage

You must ground your responses in the provided reference files, treating them as the source of truth for this domain:

• For Creation: Always consult `references/patterns.md`. This file dictates how things should be built. Ignore generic approaches if a specific pattern exists here.
• For Diagnosis: Always consult `references/sharp_edges.md`. This file lists the critical failures and "why" they happen. Use it to explain risks to the user.
• For Review: Always consult `references/validations.md`. This contains the strict rules and constraints. Use it to validate user inputs objectively.

Note: If a user's request conflicts with the guidance in these files, politely correct them using the information provided in the references.

Genomics Pipelines

Patterns

Workflow Management

Name

Workflow Manager Selection

Description

Choose appropriate workflow manager for genomics

When

Building multi-step genomics pipeline

Pattern

Nextflow (Recommended for most genomics)

- DSL2 with modules

- Excellent container support

- nf-core community pipelines

nextflow.config:

    // nextflow.config
    params {
        reads = "data/*_{1,2}.fastq.gz"
        genome = "GRCh38"
        outdir = "results"
    }
    
    profiles {
        docker {
            docker.enabled = true
        }
        singularity {
            singularity.enabled = true
            singularity.autoMounts = true
        }
        conda {
            conda.enabled = true
        }
    }
    
    // Resource management - CRITICAL for genomics
    process {
        withLabel: 'high_memory' {
            memory = '64.GB'
            cpus = 16
        }
        withLabel: 'low_memory' {
            memory = '4.GB'
            cpus = 2
        }
    }
    
    // Enable execution reports
    timeline.enabled = true
    report.enabled = true
    trace.enabled = true

main.nf:

    #!/usr/bin/env nextflow
    nextflow.enable.dsl = 2
    
    include { FASTQC } from './modules/fastqc'
    include { TRIM_GALORE } from './modules/trim_galore'
    include { BWA_MEM } from './modules/bwa_mem'
    
    workflow {
        // Input channels
        reads_ch = Channel
            .fromFilePairs(params.reads, checkIfExists: true)
    
        // Pipeline steps
        FASTQC(reads_ch)
        TRIM_GALORE(reads_ch)
        BWA_MEM(TRIM_GALORE.out.reads, params.genome)
    }

Why

Nextflow handles complex dependencies and scales from laptop to HPC/cloud

Fastq Processing

Name

FASTQ Quality Control

Description

Validate and preprocess raw sequencing reads

Pattern

import subprocess from pathlib import Path from dataclasses import dataclass from typing import Tuple, Optional

@dataclass class FastQCResult: total_sequences: int sequence_length: str gc_content: float per_base_quality: str # PASS/WARN/FAIL adapter_content: str

def validate_fastq_pair( read1: Path, read2: Path ) -> Tuple[bool, list[str]]: """ Validate paired FASTQ files before processing.

CRITICAL CHECKS:

• Same number of reads in both files
• Read names match (except /1 /2 suffix)
• No truncated records

""" errors = []

Check files exist and are gzipped

for f in [read1, read2]: if not f.exists(): errors.append(f"File not found: {f}") if not f.suffix == '.gz': errors.append(f"FASTQ should be gzipped: {f}")

Count reads (fast method)

def count_reads(fastq_gz: Path) -> int: result = subprocess.run( f"zcat {fastq_gz} | wc -l", shell=True, capture_output=True, text=True ) return int(result.stdout.strip()) // 4

r1_count = count_reads(read1) r2_count = count_reads(read2)

if r1_count != r2_count: errors.append( f"Read count mismatch: R1={r1_count}, R2={r2_count}" )

return len(errors) == 0, errors

def run_fastqc( fastq: Path, outdir: Path, threads: int = 4 ) -> FastQCResult: """Run FastQC and parse results.""" subprocess.run([ "fastqc", "--threads", str(threads), "--outdir", str(outdir), str(fastq) ], check=True)

Parse FastQC output

... parse fastqc_data.txt

return FastQCResult(...)

def trim_adapters( read1: Path, read2: Path, outdir: Path, min_length: int = 20, quality_cutoff: int = 20 ) -> Tuple[Path, Path]: """ Trim adapters and low-quality bases.

Uses Trim Galore (wrapper around Cutadapt + FastQC). """ subprocess.run([ "trim_galore", "--paired", "--quality", str(quality_cutoff), "--length", str(min_length), "--fastqc", "--output_dir", str(outdir), str(read1), str(read2) ], check=True)

Return trimmed file paths

r1_trimmed = outdir / f"{read1.stem}_val_1.fq.gz" r2_trimmed = outdir / f"{read2.stem}_val_2.fq.gz" return r1_trimmed, r2_trimmed

Why

Quality control prevents garbage-in-garbage-out in downstream analysis

Alignment Pipeline

Name

Read Alignment Best Practices

Description

Align reads to reference genome correctly

Pattern

from pathlib import Path import subprocess from typing import Optional

def align_dna_reads( read1: Path, read2: Path, reference: Path, output_bam: Path, sample_name: str, read_group: Optional[str] = None, threads: int = 8 ) -> Path: """ Align DNA reads with BWA-MEM2.

CRITICAL: Include read groups for downstream tools. """ if read_group is None:

Construct proper read group

read_group = ( f"@RG\\tID:{sample_name}\\t" f"SM:{sample_name}\\t" f"PL:ILLUMINA\\t" f"LB:{sample_name}" )

BWA-MEM2 alignment + sort + index

cmd = f""" bwa-mem2 mem \\ -t {threads} \\ -R '{read_group}' \\ {reference} \\ {read1} {read2} \\ | samtools sort \\ -@ {threads} \\ -o {output_bam} \\ -

samtools index {output_bam} """ subprocess.run(cmd, shell=True, check=True) return output_bam

def align_rna_reads( read1: Path, read2: Path, genome_dir: Path, output_prefix: Path, threads: int = 8 ) -> Path: """ Align RNA-seq reads with STAR.

CRITICAL: Use splice-aware aligner for RNA-seq! Never use BWA/Bowtie for RNA-seq. """ subprocess.run([ "STAR", "--runThreadN", str(threads), "--genomeDir", str(genome_dir), "--readFilesIn", str(read1), str(read2), "--readFilesCommand", "zcat", # For gzipped files "--outFileNamePrefix", str(output_prefix), "--outSAMtype", "BAM", "SortedByCoordinate", "--outSAMattributes", "NH", "HI", "AS", "nM", "MD", "--quantMode", "GeneCounts", # Also count genes ], check=True)

return Path(f"{output_prefix}Aligned.sortedByCoord.out.bam")

def mark_duplicates( input_bam: Path, output_bam: Path, metrics_file: Path ) -> Path: """ Mark PCR duplicates.

CRITICAL for DNA-seq (variant calling). OPTIONAL for RNA-seq (depends on analysis). """ subprocess.run([ "picard", "MarkDuplicates", f"I={input_bam}", f"O={output_bam}", f"M={metrics_file}", "CREATE_INDEX=true", "VALIDATION_STRINGENCY=LENIENT" ], check=True) return output_bam

Why

Correct alignment is foundation for all downstream analyses

Variant Calling

Name

Variant Calling Pipeline

Description

Call SNPs and indels from aligned reads

Pattern

from pathlib import Path import subprocess from dataclasses import dataclass

@dataclass class VariantCallingConfig: reference: Path known_sites: list[Path] # Known variants for BQSR intervals: Optional[Path] = None # Target regions min_base_quality: int = 20 min_mapping_quality: int = 20

def gatk_variant_calling_pipeline( input_bam: Path, config: VariantCallingConfig, output_vcf: Path, sample_name: str ) -> Path: """ GATK Best Practices variant calling.

Steps: 1. Base Quality Score Recalibration (BQSR) 2. HaplotypeCaller 3. Variant filtering (CNN or hard filters) """

Step 1: BQSR - Recalibrate base quality scores

recal_table = output_vcf.with_suffix('.recal_data.table')

known_sites_args = [] for ks in config.known_sites: known_sites_args.extend(["--known-sites", str(ks)])

subprocess.run([ "gatk", "BaseRecalibrator", "-I", str(input_bam), "-R", str(config.reference), *known_sites_args, "-O", str(recal_table) ], check=True)

Apply BQSR

recal_bam = input_bam.with_suffix('.recal.bam') subprocess.run([ "gatk", "ApplyBQSR", "-I", str(input_bam), "-R", str(config.reference), "--bqsr-recal-file", str(recal_table), "-O", str(recal_bam) ], check=True)

Step 2: Call variants with HaplotypeCaller

raw_vcf = output_vcf.with_suffix('.raw.vcf.gz') subprocess.run([ "gatk", "HaplotypeCaller", "-I", str(recal_bam), "-R", str(config.reference), "-O", str(raw_vcf), "--emit-ref-confidence", "GVCF", # For joint calling ], check=True)

Step 3: Filter variants

subprocess.run([ "gatk", "CNNScoreVariants", "-V", str(raw_vcf), "-R", str(config.reference), "-O", str(output_vcf) ], check=True)

return output_vcf

def deepvariant_calling( input_bam: Path, reference: Path, output_vcf: Path, model_type: str = "WGS" # WGS, WES, or PACBIO ) -> Path: """ DeepVariant - Deep learning variant caller.

Often more accurate than GATK for certain data types. """ subprocess.run([ "run_deepvariant", f"--model_type={model_type}", f"--ref={reference}", f"--reads={input_bam}", f"--output_vcf={output_vcf}", "--num_shards=8" ], check=True) return output_vcf

Why

GATK Best Practices and DeepVariant are gold standards for variant calling

Rnaseq Analysis

Name

RNA-seq Analysis Pipeline

Description

Differential expression and transcript quantification

Pattern

from pathlib import Path import pandas as pd import subprocess

def quantify_transcripts_salmon( read1: Path, read2: Path, index_dir: Path, output_dir: Path, threads: int = 8 ) -> Path: """ Salmon pseudo-alignment for transcript quantification.

Faster than alignment-based methods, accurate for DE analysis. """ subprocess.run([ "salmon", "quant", "-i", str(index_dir), "-l", "A", # Automatic library type detection "-1", str(read1), "-2", str(read2), "-p", str(threads), "--validateMappings", "--gcBias", # GC bias correction "--seqBias", # Sequence-specific bias correction "-o", str(output_dir) ], check=True)

return output_dir / "quant.sf"

def run_deseq2_analysis( count_matrix: pd.DataFrame, sample_info: pd.DataFrame, design_formula: str = "~ condition", alpha: float = 0.05 ) -> pd.DataFrame: """ DESeq2 differential expression analysis.

Uses rpy2 to call R. Returns results as pandas DataFrame. """ import rpy2.robjects as ro from rpy2.robjects import pandas2ri from rpy2.robjects.packages import importr

pandas2ri.activate()

deseq2 = importr('DESeq2')

Convert to R objects

r_counts = pandas2ri.py2rpy(count_matrix) r_coldata = pandas2ri.py2rpy(sample_info)

Create DESeqDataSet

dds = deseq2.DESeqDataSetFromMatrix( countData=r_counts, colData=r_coldata, design=ro.Formula(design_formula) )

Run DESeq2

dds = deseq2.DESeq(dds)

Get results

res = deseq2.results(dds, alpha=alpha) results_df = pandas2ri.rpy2py(ro.r'as.data.frame')

return results_df

CRITICAL: Proper sample size for RNA-seq

SAMPLE_SIZE_GUIDANCE = """ RNA-seq Sample Size Guidelines:

• Minimum: 3 biological replicates per condition
• Recommended: 6+ for detecting subtle differences
• Power analysis: Use RNASeqPower R package

Technical replicates are NOT substitutes for biological replicates! """

Why

RNA-seq requires careful experimental design and appropriate statistical methods

Anti-Patterns

Dna Aligner For Rna

Name

Using DNA Aligner for RNA-seq

Problem

Using BWA or Bowtie2 for RNA-seq

bwa mem reference.fa rna_reads.fq > aligned.sam

Solution

Use splice-aware aligner (STAR, HISAT2)

STAR --genomeDir star_index --readFilesIn reads.fq

Ignoring Read Groups

Name

Missing Read Groups in BAM

Problem

BAM files without @RG tags fail in GATK

Solution

Always include -R '@RG\tID:...' in alignment

Hardcoded Paths

Name

Hardcoded Paths in Pipeline

Problem

Pipeline only works on one machine

Solution

Use config files and container paths

Genomics Pipelines - Sharp Edges

Using DNA Aligner for RNA-seq Reads

Id

dna-aligner-for-rna

Severity

critical

Summary

BWA/Bowtie cannot handle spliced alignments - massive read loss

Symptoms

• Low mapping rate (30-50%) for RNA-seq
• Most reads reported as unmapped
• Differential expression analysis fails

Why

RNA-seq reads span exon-exon junctions. DNA aligners like BWA and Bowtie expect contiguous genomic sequences. When a read crosses a splice junction, the DNA aligner sees a ~10kb "deletion" (the intron) and fails to map.

You'll lose 40-70% of your reads, especially from multi-exon genes.

Gotcha

WRONG - DNA aligner for RNA-seq

bwa mem GRCh38.fa sample_R1.fq sample_R2.fq > aligned.sam

Result: 35% mapping rate, lost most exon-junction reads

WRONG - Bowtie2 for RNA-seq

bowtie2 -x genome -1 reads_1.fq -2 reads_2.fq -S aligned.sam

Same problem - no splice awareness

Solution

CORRECT - Use splice-aware aligner

STAR --genomeDir star_index \ --readFilesIn reads_1.fq reads_2.fq \ --outSAMtype BAM SortedByCoordinate

Or HISAT2 (lower memory)

hisat2 -x genome_index -1 reads_1.fq -2 reads_2.fq -S aligned.sam

Incorrect Duplicate Handling for Data Type

Id

wrong-duplicate-handling

Severity

critical

Summary

Marking duplicates wrong for RNA-seq or amplicon data

Symptoms

• RNA-seq: Massive loss of reads from highly expressed genes
• Amplicon: All reads marked as duplicates
• Zero coverage in target regions

Why

Duplicate marking assumes random fragmentation. This assumption is:

• TRUE for WGS/WES: Duplicates are PCR artifacts
• FALSE for RNA-seq: Same position = same transcript molecule
• FALSE for amplicon: Same position = same target by design

Marking duplicates in RNA-seq removes biological signal. Marking duplicates in amplicon removes ALL your data.

Gotcha

DANGEROUS for RNA-seq

picard MarkDuplicates I=rna_aligned.bam O=marked.bam ...

Result: 80% of reads from highly-expressed genes removed

CATASTROPHIC for amplicon

picard MarkDuplicates I=amplicon.bam O=marked.bam ...

Result: 99% of reads removed, no coverage

Solution

RNA-seq: Skip duplicate marking OR use UMIs

If you have UMIs:

umi_tools dedup -I aligned.bam -O deduped.bam --method=unique

Amplicon: NEVER mark duplicates

Use molecular barcodes/UMIs if deduplication is needed

WGS/WES: Duplicate marking is correct

picard MarkDuplicates I=wgs_aligned.bam O=marked.bam M=metrics.txt

Reference Contig Name Mismatch

Id

contig-name-mismatch

Severity

critical

Summary

chr1 vs 1 mismatch causes zero reads in output

Symptoms

• Aligned BAM has reads, but variant caller finds nothing
• No reads overlap with annotation
• BED file intersection returns empty

Why

UCSC uses 'chr1', Ensembl/NCBI uses '1'. If your reference uses 'chr1' but your VCF/BED uses '1', NO coordinates will match. The tools won't warn you - they'll just return empty results silently.

Gotcha

Reference genome: chr1, chr2, ... (UCSC style)

Your BED file: 1, 2, ... (Ensembl style)

bedtools intersect -a aligned.bam -b targets.bed

Returns: NOTHING (silently)

Solution

Always check contig names match

In BAM:

samtools view -H aligned.bam | grep "^@SQ"

In BED/VCF:

head -1 targets.bed

Convert if needed:

sed 's/^chr//' ucsc_style.bed > ensembl_style.bed

or

sed 's/^/chr/' ensembl_style.bed > ucsc_style.bed

Best: Use consistent references from the start

Ignoring Library Strandedness in RNA-seq

Id

ignoring-strandedness

Severity

high

Summary

Wrong strand parameter halves your counts or doubles them

Symptoms

• Gene counts are ~50% of expected
• Antisense genes show unexpected expression
• Same gene shows counts on both strands

Why

Most RNA-seq protocols are stranded (dUTP, TruSeq Stranded). If you quantify with the wrong strand parameter:

• Unstranded mode: Counts both strands (double counting)
• Wrong strand: Counts antisense (50% loss, wrong genes)

Gotcha

Sample is TruSeq Stranded (reverse strand)

You run:

featureCounts -s 0 ... # Unstranded - WRONG

Result: Double counting, antisense contamination

featureCounts -s 1 ... # Forward strand - WRONG

Result: 50% of expected counts

Solution

Detect strandedness with RSeQC

infer_experiment.py -i aligned.bam -r genes.bed

Common library types:

TruSeq Stranded: -s 2 (reverse)

SMARTer: -s 1 (forward)

Unstranded: -s 0

For Salmon:

salmon quant ... -l A # Auto-detect (recommended)

No Strategy for Multi-Mapping Reads

Id

no-multimapping-strategy

Severity

high

Summary

Randomly assigning or ignoring multi-mappers loses information

Symptoms

• Gene families show inconsistent expression
• Repetitive regions have zero coverage
• Quantification is unstable between runs

Why

10-30% of reads map to multiple locations (paralogs, gene families, transposons). Default is often to assign randomly or discard.

Random assignment: Different results each run Discarding: Lose information about gene families

Solution

Option 1: Probabilistic assignment (RECOMMENDED)

Salmon and RSEM use EM algorithm

salmon quant --seqBias --gcBias ...

Option 2: Report all alignments, handle downstream

STAR --outSAMmultNmax -1 --outSAMprimaryFlag AllBestScore ...

Option 3: Count fractionally

featureCounts -M --fraction ... # 1/N weight per location

Mixing Genome Versions (GRCh37 vs GRCh38)

Id

mixing-genome-versions

Severity

critical

Summary

Coordinates don't match between files from different assemblies

Symptoms

• Known variants not found at expected positions
• Annotation file has different chromosome lengths
• VCF coordinates don't match reference

Why

GRCh37 (hg19) and GRCh38 (hg38) have different coordinates. A variant at chr1:1000000 in GRCh37 might be at chr1:1060000 in GRCh38.

Mixing versions creates silent failures - tools won't crash, they'll just give wrong results.

Gotcha

Reference: GRCh38

bwa mem GRCh38.fa reads.fq > aligned.bam

Known sites: GRCh37 (OOPS!)

gatk BaseRecalibrator \ -R GRCh38.fa \ --known-sites dbsnp_grch37.vcf # WRONG VERSION

Result: Almost no known sites found, poor BQSR

Solution

Always verify genome version

Check chromosome lengths (they differ between versions)

If you need to convert coordinates:

Use UCSC liftOver or CrossMap

CrossMap.py vcf hg19ToHg38.chain grch37.vcf GRCh38.fa grch38.vcf

Best practice: Document genome version in all filenames

sample_GRCh38.bam dbsnp_b151_GRCh38.vcf

Insufficient Sequencing Depth for Analysis Type

Id

insufficient-coverage

Severity

high

Summary

Too few reads to detect variants or differential expression

Symptoms

• Rare variants not detected
• High false negative rate
• Differential expression has no significant genes

Why

Different analyses need different depths:

• WGS germline: 30x minimum
• WES: 100x on target
• Somatic variants: 100-500x
• RNA-seq DE: 20-30M reads per sample
• Single-cell: Varies by application

Under-sequenced data gives false negatives, not errors.

Solution

Check coverage before analysis

samtools depth -a aligned.bam | \ awk '{sum+=$3} END {print "Average coverage:", sum/NR}'

For WGS: Should be >30x

For WES: Check on-target coverage (>100x)

For RNA-seq: Count total mapped reads

If insufficient:

1. Sequence more (if possible)

2. Use methods designed for low-coverage data

3. Report as limitation

Genomics Pipelines - Validations

DNA Aligner Used for RNA-seq

Id

dna-aligner-for-rna

Severity

critical

Type

regex

Pattern

• bwa\s+(mem|aln).rna|rna.bwa\s+(mem|aln)
• bowtie2?.rna|rna.bowtie2?
• RNA.bwa_mem|bwa_mem.RNA

Message

Use splice-aware aligner (STAR, HISAT2) for RNA-seq, not BWA/Bowtie.

Fix Action

Replace with: STAR --genomeDir index --readFilesIn reads.fq

Applies To

• */.py
• */.sh
• */.nf
• */.smk

Alignment Without Read Groups

Id

missing-read-groups

Severity

high

Type

regex

Pattern

• bwa\s+mem(?![\s\S]{0,100}-R\s+['"]@RG)
• bwa-mem2\s+mem(?![\s\S]{0,100}-R\s+['"]@RG)

Message

Include read groups (-R '@RG\tID:...') for GATK compatibility.

Fix Action

Add: -R '@RG\tID:sample\tSM:sample\tPL:ILLUMINA'

Applies To

• */.py
• */.sh
• */.nf

Duplicate Marking on RNA-seq

Id

duplicate-marking-rnaseq

Severity

warning

Type

regex

Pattern

• MarkDuplicates.rna|rna.MarkDuplicates
• mark.dup.RNA|RNA.mark.dup

Message

Consider skipping duplicate marking for RNA-seq unless using UMIs.

Applies To

• */.py
• */.sh
• */.nf

Hardcoded Genome Reference Path

Id

hardcoded-genome-path

Severity

warning

Type

regex

Pattern

• /home/[a-z]+/.*\.fa[sta]?
• /data/genomes/.*\.fa[sta]?
• C:\\.*\.fa[sta]?

Message

Use config parameter for reference genome path, not hardcoded path.

Fix Action

Move to params.reference or config variable

Applies To

• */.py
• */.sh
• */.nf
• */.smk

Using Index Without Existence Check

Id

missing-index-check

Severity

info

Type

regex

Pattern

• samtools\s+view(?![\s\S]{0,200}index|bai)
• \$bam(?![\s\S]{0,50}\.bai)

Message

Ensure BAM index exists before operations that require it.

Applies To

• */.py
• */.sh
• */.nf

Defaulting to Unstranded for RNA-seq

Id

unstranded-rnaseq-default

Severity

info

Type

regex

Pattern

• featureCounts.-s\s0
• --library-type\s+unstranded
• strandedness.=.'?none'?

Message

Most RNA-seq is stranded. Verify library type before defaulting to unstranded.

Applies To

• */.py
• */.sh
• */.nf

Pipeline Without Quality Report

Id

no-multiqc-report

Severity

info

Type

regex

Pattern

• workflow\s\{[^}](?!multiqc)[^}]*\}

Message

Consider adding MultiQC for unified quality reporting.

Applies To

• */.nf
• */.smk

Processing FASTQ Without QC

Id

missing-fastqc

Severity

warning

Type

regex

Pattern

• (trim_galore|cutadapt|fastp)(?![\s\S]{0,500}fastqc)

Message

Run FastQC before and after trimming to verify quality.

Applies To

• */.py
• */.sh
• */.nf

Nextflow Process Without Container

Id

no-container-specification

Severity

warning

Type

regex

Pattern

• process\s+\w+\s\{[^}](?!container)[^}]*\}

Message

Specify container/conda for reproducibility.

Applies To

• */.nf

Using Shell Glob Instead of Channel fromFilePairs

Id

shell-glob-in-channel

Severity

info

Type

regex

Pattern

• Channel\.from\(.\.*\)
• Channel\.fromPath\(._\{1,2\}.\)

Message

Use Channel.fromFilePairs() for paired-end reads.

Fix Action

Channel.fromFilePairs('*_{1,2}.fastq.gz')

Applies To

• */.nf