Getting Started With Raw Sequencing Reads
Most people coming into this field skip straight to calling variants because they watched a YouTube tutorial, then wonder why their results look like garbage. I've been doing this since we were working with 10x coverage human genomes and still lose sleep over base quality distributions. The gap between what the sequencer gives you and what you actually trust is where the real work lives. You start with FASTQ files. I know, that's obvious, but the mistake everyone makes is assuming the quality scores in there are reliable. They're not always. Base Quality Score Recalibration, commonly called BQSR in GATK terminology, exists for exactly this reason. The sequencer systematically overestimates the quality of certain contexts, particularly around homopolymers in Illumina data and in Kmer-specific patterns. Skipping BQSR is the single most common source of false positive variant calls I see in published work. After you've generated your FASTQ, the first real step is adapter trimming. You might think your library prep was clean, but Illumina's Nextera and TruSeq adapters routinely survive into the reads, especially in smaller fragments. I use fastp for this because it handles both adapter removal and quality filtering in one pass, usually cutting 2-5 percent of reads per sample depending on the library. That seems small until you're working with thirty samples and trying to maintain uniform coverage across a genome.
Alignment comes next. For human data, BWA-MEM2 is the current standard if you have the RAM for it. The original BWA-MEM algorithm will align a 30x human genome in roughly 45 minutes on a decent machine, while BWA-MEM2 gets that down to about 15-20 minutes with comparable accuracy. If you're working with non-human organisms, BWA-MEM works fine for moderate divergence, but once you hit more than 5 percent sequence divergence from your reference, you should switch to minimap2. It handles splice-aware alignment for RNA-seq data natively and deals with structural variation in long-read alignments without breaking a sweat. Post-alignment processing is where people lose entire weekends. Marking duplicates with Picard or sambamba is non-negotiable if you're doing variant calling. PCR duplicates inflate your allele frequency estimates and create the illusion of higher confidence at positions that are actually artifacts. When I process whole exome data, duplicate rates of 15-25 percent are typical. With whole genome sequencing, you're usually looking at 5-10 percent. Anything above that range suggests either low input DNA or a library prep issue that should be flagged before you proceed.
The Variant Calling Step
GATK's HaplotypeCaller operates in GVCF mode for cohort analysis. This is important because joint genotyping produces dramatically different results than calling samples in isolation. When I set up a cohort, I run HaplotypeCaller with the recommended GATK best practices, generating a GVCF for each sample, then genotype them together. Running genotyping on individual samples and merging the VCFs later introduces systematic biases that show up as batch effects in downstream analysis. The hard part is not calling variants. It's filtering them. GATK's Variant Quality Score Recalibration is technically superior to hard filtering, but it requires a large training set with known variant resources like HapMap, Omni, 1000 Genomes, and dbSNP. If you're working with a non-model organism without these resources, you fall back to hard filtering. The typical approach uses QD below 2.0, FS above 60.0, MQ below 40.0, and MQRankSum below -12.5 as rough cutoffs. These numbers are guidelines, not rules, and you should always visualize your distributions before applying them. I ran into a specific problem last year working with a panel of bacterial isolates that completely broke standard variant calling assumptions. The genome had a repetitive insertion sequence element creating false structural variants that confused GATK's local reassembly. The haplotypes it was constructing didn't match anything biological. I solved this by filtering out regions with mapping quality below 30 and using a custom BED file to mask the known repetitive elements before running the caller. If you're working with organisms that have high repeat content, building a repeat mask from RepeatMasker output beforehand can prevent weeks of troubleshooting.
Get the Full Details
RNA-Seq Specific Considerations
RNA-seq adds a completely different layer of complexity because spliced alignment changes the error model. Star and Hisat2 are the two mainstream aligners. Star is faster and generally more accurate for bulk RNA-seq, while Hisat2 uses less memory and handles novel splice sites better in non-model organisms. Quantification is where things get interesting. FeatureCounts and HTSeq both work, but Salmon and Kallisto provide transcript-level quantification that's substantially faster than full alignment-based approaches. In practice, I use Salmon in alignment-based mode when I need transcript-level resolution or isoform switching analysis, and featureCounts when I just need gene-level counts for differential expression. Batch effects in RNA-seq are brutal. They often exceed biological signal. If you're processing samples across multiple sequencing runs, different flow cell lanes, or different library prep dates, you need to include batch as a covariate in your differential expression model. I've seen labs publish findings that turned out to be entirely driven by the fact that all their treated samples were sequenced on one lane and controls on another. Combat and RUVseq are the standard correction methods, but they're not magic. They adjust for known and estimated batch factors, but they cannot recover signal that was never captured in the first place.
Long Read Analysis
PacBio and Oxford Nanopore data require a fundamentally different workflow. These platforms have higher raw error rates, though HiFi reads from PacBio now achieve 99.9 percent accuracy natively. For Nanopore, basecalling with Dorado or Guppy and then using tools like minimap2 for alignment and Medaka for polishing is the standard pipeline. Structural variant detection from long reads is where this technology really shines, identifying insertions, deletions, inversions, and translocations that short reads simply cannot resolve. The downside is computational cost. A single Nanopore run generating 50GB of FAST5 data takes significantly more processing than an equivalent Illumina run. Consensus polishing typically requires 3-5 rounds with tools like Pilon or Racon, and variant calling against a reference genome with long reads demands coverage of at least 30x for reliable SNV calling and 15-20x for structural variants. The accuracy improves non-linearly with coverage, so 50x is a safer target for clinical applications.
Practical Workflow Advice
The biggest time sink in any sequencing project is not the analysis itself, it's debugging environment issues and parameter mismatches. I use Nextflow or Snakemake for every project now, even small ones. A properly written pipeline takes maybe an hour to set up but saves hours every time you need to rerun something with new parameters or additional samples. Docker containers for tools like GATK, BWA, and Star eliminate the version drift problem that silently corrupts reproducibility. Storage is another issue nobody plans for adequately. A single 30x human whole genome produces roughly 150GB of raw FASTQ data, and after alignment, sorting, and indexing, your BAM file will be around 90-100GB. A GVCF for the same sample is another 20-30GB. If you're running a cohort study, your storage needs scale linearly and fast. I budget roughly 500GB per sample for complete raw and processed data over the lifecycle of a WGS project. Download links for the tools I mentioned are all available from their respective project pages. BWA-MEM2 from the GitHub repository, fastp from bioconda, GATK from Broad Institute's distribution site, Star from GitHub, and Salmon from the pattonlab GitHub. Using conda or bioconda for dependency management is strongly recommended, as these tools have many interdependent libraries that conflict with system installations.
Quality control at every stage is essential. MultiQC aggregates QC metrics from every tool in your pipeline into a single report, and it's the first thing I look at before making any decisions. FastQC reports alone are sufficient for initial assessment, but MultiQC lets you compare metrics across dozens of samples simultaneously and spot outliers that individual reports hide. A sample that looks normal in isolation but deviates from the cohort in aggregate metrics is almost always a problematic sample worth investigating before proceeding further.
