How to Analyze Your Own FASTQ Files

By the end of this guide you will have taken a pair of FASTQ files off disk and verified their integrity and encoding. You will have produced a quality report you can read, trimmed adapters, and aligned the reads to GRCh38. You will also have a sorted, indexed and duplicate-marked CRAM, with coverage statistics that tell you whether the data supports variant calling. You need a Linux or macOS machine with at least 32 GB of RAM and roughly 500 GB of free disk for a 30x human whole genome, plus seqkit or seqfu, fastp, bwa-mem2 or minimap2, samtools, and mosdepth. All of these install cleanly from bioconda. Expect eight to twenty hours of wall time for a human genome on sixteen cores, most of it in alignment. If you only have exome or RNA-seq data, every step below still applies, but the alignment tool and the coverage expectations change, and I note where.

1. Read the format before you trust any tool

A FASTQ file is four lines per read. The first begins with @ and holds the read identifier, the second holds the bases, the third begins with +, and the fourth holds one ASCII quality character per base. By contrast, a FASTA file has two parts per record, a > header and the sequence. It carries no quality information at all. That is the whole difference. FASTA is what you use for a reference genome or a protein sequence. FASTQ is what a sequencer emits, because it needs somewhere to put the per-base confidence.

Here is a single record from an Illumina paired-end run:

@A00123:456:HWTLMDSXX:1:1101:2013:1000 1:N:0:ATCACGTT+GGTACCAA
NTGGCTCAGCTTTAACGCGTGCTTTACGGTATCGCCGCTCCCGATTCGTAGTGCAGCTCG
+
#FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFF,FFFF,

The header after the space encodes the read number in the pair, whether the read was filtered, the control bits, and the index sequences. The quality line is Phred+33: subtract 33 from the ASCII code to get a Phred score Q, where the probability of the base being wrong is 10^(-Q/10). So F (ASCII 70) is Q37, roughly one error in five thousand. The leading # is Q2, Illumina’s marker for a base it does not stand behind, which is why the corresponding base is N.

This matters because the encoding has not always been the same. Solexa and early Illumina pipelines used Phred+64, and Solexa used a different score definition entirely that diverges from Phred at low quality values, which is exactly where the difference bites ¹. Anything sequenced after about 2011 is Phred+33. If you inherit old data or data from an unfamiliar source, check before you run a quality filter. A Phred+64 file read as Phred+33 will look like it has impossibly high quality, and your filters will pass everything.

A one-liner to inspect a gzipped file without decompressing it to disk:

zcat sample_R1.fastq.gz | head -n 8

Use gzcat on macOS. Never open a multi-gigabyte FASTQ in a text editor. If you want a paged view, zless sample_R1.fastq.gz works and does not load the file into memory.

2. Verify integrity and count reads

Truncated downloads are the single most common cause of a pipeline that fails four hours in. Check the gzip stream first, because this is fast and catches the majority of transfer errors:

gzip -t sample_R1.fastq.gz && echo OK

Then count records. The arithmetic approach works and is worth knowing, since every FASTQ record occupies exactly four lines:

echo $(( $(zcat sample_R1.fastq.gz | wc -l) / 4 ))

This fails silently if the file is malformed, so for anything you care about, use a parser that validates structure. seqkit stats gives you counts, total bases, and length distribution in one pass:

seqkit stats -a sample_R1.fastq.gz sample_R2.fastq.gz

The -a flag adds N50, Q20 and Q30 percentages, and average quality. For paired data, seqfu check is the tool I reach for, because it validates the four-line structure and confirms that R1 and R2 have matching identifiers in the same order, which is the specific failure that breaks aligners in confusing ways ². Read counts in R1 and R2 must be identical. If they are not, something upstream subsetted or filtered one file independently, and you should go back to the source rather than try to repair it.

Sanity-check the total base count against what you expect. A 30x human genome is about 90 gigabases, so at 2x150 bp you expect roughly 300 million read pairs. If seqkit stats reports 60 gigabases, you have a 20x genome, and you should know that before you interpret anything downstream.

3. Run FastQC and read the plots that matter

FastQC produces an HTML report per file in a few minutes:

fastqc -t 8 -o qc/ sample_R1.fastq.gz sample_R2.fastq.gz
multiqc qc/ -o qc/

MultiQC collapses multiple reports into one page, which is what you want if your data arrived as several lanes.

Four panels carry most of the information. Per-base sequence quality should sit above Q30 across the read, with a gradual decline toward the 3’ end that is normal for Illumina chemistry. A sharp cliff partway through the read suggests a run problem. Per-base sequence content should be flat after the first ten to fifteen bases, where random-priming bias in library prep reliably produces a wobble that you can ignore. Adapter content tells you how much of your library was shorter than the read length, which determines how aggressively you need to trim. Overrepresented sequences will often name the adapter directly, and for RNA-seq will surface highly expressed transcripts, which is expected rather than alarming.

Treat FastQC’s red and green icons as prompts rather than verdicts. A whole genome library will fail “per sequence GC content” if there is any contamination, but a normal amplicon panel will fail “sequence duplication levels” by design. Read the plot and do not rely on the badge.

4. Trim adapters and low-quality tails with fastp

I use fastp rather than Trimmomatic or cutadapt for routine work, because it detects adapters automatically from paired-end overlap, corrects mismatched overlapping bases, and emits its own HTML and JSON report in a single pass over the data at speeds that are typically I/O bound rather than CPU bound.

fastp \
  -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o trim_R1.fastq.gz -O trim_R2.fastq.gz \
  --detect_adapter_for_pe \
  --correction \
  --cut_tail --cut_tail_mean_quality 20 \
  --length_required 50 \
  --thread 16 \
  --json fastp.json --html fastp.html

A word on how hard to trim. For alignment-based work, quality trimming beyond removing adapters and clearly degenerate tails is usually counterproductive, because modern aligners soft-clip poorly matching ends and variant callers weight bases by quality anyway. Aggressive trimming shortens reads, which lowers mapping uniqueness in repetitive regions and costs you coverage where you least want to lose it. For de novo assembly or k-mer based methods, the calculus flips and harder trimming helps. The settings above are deliberately gentle.

Check the fastp JSON afterward. If more than about 20% of bases were removed, investigate the library rather than accepting the result. A very high adapter fraction means the insert size distribution was shorter than the read length, which limits the effective information content of the run regardless of how many reads you have.

5. Align to GRCh38 and produce a sorted CRAM

For short reads, bwa-mem2 is the current default: same algorithm and output as bwa mem, roughly two to three times faster. Use an analysis-set reference that includes alt contigs and the decoy sequence, because the decoy absorbs reads from sequences absent from the primary assembly and reduces false positive variant calls in the regions those reads would otherwise pile into.

bwa-mem2 index GRCh38_full_analysis_set_plus_decoy_hla.fa

bwa-mem2 mem -t 16 -K 100000000 -Y \
  -R '@RG\tID:run1\tSM:sample01\tPL:ILLUMINA\tLB:lib1' \
  GRCh38_full_analysis_set_plus_decoy_hla.fa \
  trim_R1.fastq.gz trim_R2.fastq.gz \
| samtools sort -@ 8 -m 2G -o sample.bam -

The -R read group is not optional in practice. GATK and most downstream callers refuse to run without it, and adding it later means rewriting the file. The -K 100000000 setting fixes the chunk size so output is deterministic across thread counts, which you want if you ever need to reproduce a result. The -Y flag uses soft clipping for supplementary alignments, which preserves the sequence for structural variant callers.

Then mark duplicates and convert to CRAM, which stores the same information referenced against the genome and typically costs 30 to 60 percent of the equivalent BAM:

samtools markdup -@ 8 --reference GRCh38_full_analysis_set_plus_decoy_hla.fa \
  sample.bam sample.cram
samtools index sample.cram

samtools markdup needs the file processed by samtools fixmate -m first if you sorted by coordinate directly. The simpler route is samtools collate then fixmate -m then sort then markdup. Keep the reference FASTA forever. CRAM is a reference-relative encoding, and without the exact reference you used, the file is not decodable.

For long reads from Oxford Nanopore or PacBio, substitute minimap2 -ax map-ont or -ax map-hifi. The tradeoffs differ substantially. Long reads resolve structural variants and repetitive regions that short reads cannot span, and whole-genome long-read sequencing has been used to characterize structural variation whose breakpoints and origins are inaccessible to short-read data ³. Long reads also carry higher per-base error, so your quality expectations from step 3 do not transfer. Long-read FASTQ files are also larger per base and compress differently, which is why specialized compressors reorder reads and model quality streams separately rather than relying on generic gzip ⁴.

6. Measure coverage and decide whether the data supports what you want to do

Alignment statistics are the point at which you learn whether the run succeeded.

samtools flagstat -@ 8 sample.cram
mosdepth --by 1000 --fast-mode -t 4 sample sample.cram

For a human whole genome, expect above 99% of reads mapped, above 95% properly paired, and a duplicate rate below 10% for a PCR-free library. mosdepth writes a sample.mosdepth.summary.txt with mean depth per chromosome and a sample.mosdepth.global.dist.txt from which you can read the fraction of the genome at or above any depth. The number I look at first is the fraction of callable genome at 10x or greater, because that is roughly where confident germline single-nucleotide variant calling begins. For a good 30x genome, that fraction should exceed 0.95.

Mean depth alone is misleading, because coverage is not uniform. GC-extreme regions, segmental duplications, and the first exons of many genes are systematically underrepresented. Two libraries with identical mean depth can differ substantially in how much of the genome is callable, which is why the distribution matters more than the average.

If chromosome-level mean depths in the mosdepth summary are grossly uneven, or if X and Y depths do not match the expected pattern for the sample, stop and work out why before proceeding.

7. Search a sequence against a database when you need to identify something

Occasionally you will pull an overrepresented sequence out of a FastQC report or an unmapped read out of a BAM and want to know what it is. Convert to FASTA and run a similarity search:

seqkit fq2fa unmapped.fastq.gz > unmapped.fasta

BLAST is the usual first stop, but the FASTA family of programs remains useful. It helps particularly for translated comparisons, where you are matching a DNA query against protein sequences and need to handle frameshifts. The programs fastx and tfastx handle them explicitly rather than by reporting disjoint high-scoring segments ⁵. The practical distinction is that BLAST optimizes for speed on large databases while FASTA’s algorithms produce full local alignments with statistical estimates that behave well for the multiple-alignment and homology-assessment tasks that follow a search ⁶. For a handful of sequences, run both and compare.

Interpretation stops here. Sequence identity tells you what a read resembles. It does not tell you what the read means for the person the sample came from. Any question about clinical significance belongs with a genetic counselor or physician working from a validated clinical assay. A research-grade pipeline you ran yourself is not a basis for that answer.

Common problems

Read counts differ between R1 and R2. Something filtered the files independently. Re-derive both from the original BCL or the original download rather than attempting to re-pair them, unless you have a tool that re-syncs by identifier and you accept the lost reads.

Quality scores look implausibly high and nothing fails QC. Check the encoding. A Phred+64 file interpreted as Phred+33 shifts every score up by 31, and the Solexa variant’s score definition differs from Phred below about Q10 in a way that no linear offset fixes ¹. seqkit stats -a will report an average quality above 40, which no real Illumina run achieves.

Alignment produces a huge SAM file that fills the disk. Always pipe bwa-mem2 mem output directly into samtools sort. Uncompressed SAM for a human genome runs to hundreds of gigabytes, and writing it to disk is pure waste.

Storage grows faster than expected. FASTQ compresses poorly with gzip because the quality stream is high-entropy, and it is often the dominant cost in a sequencing archive. Domain-specific compressors exploit the structure that gzip cannot see, reordering reads by similarity and modeling the two-dimensional redundancy across reads and positions, with meaningful gains over general-purpose tools ⁷. For long-read data the gains are larger still ⁴. In practice, my recommendation is to convert to CRAM after alignment, keep the reference alongside it, and retain FASTQ only if you expect to realign to a future reference build.

Very low mapping rate with high-quality reads. Either the reference is wrong or the library contains something unexpected. Take a thousand unmapped reads and run them through the search in step 7. Contamination, adapter dimers, and index hopping all produce distinctive signatures. Enrichment-based libraries have their own characteristic behavior, where the on-target fraction rather than the raw mapping rate is the number that describes whether capture worked ⁸.

Footnotes

1. Peter J. A. Cock, Christopher J. Fields, Naohisa Goto, et al. The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants. Nucleic Acids Research, 2009. https://doi.org/10.1093/nar/gkp1137 ↩ ↩²

2. Andrea Telatin, Piero Fariselli, Giovanni Birolo. SeqFu: A Suite of Utilities for the Robust and Reproducible Manipulation of Sequence Files. Bioengineering, 2021. https://doi.org/10.3390/bioengineering8050059 ↩

3. Akihiro Fujimoto, Jing Hao Wong, Yukiko Yoshii, et al. Whole-genome sequencing with long reads reveals complex structure and origin of structural variation in human genetic variations and somatic mutations in cancer. Genome Medicine, 2021. https://doi.org/10.1186/s13073-021-00883-1 ↩

4. Dohyeon Lee, Giltae Song. FastqCLS: a FASTQ compressor for long-read sequencing via read reordering using a novel scoring model. Bioinformatics, 2021. https://doi.org/10.1093/bioinformatics/btab696 ↩ ↩²

5. William R. Pearson, Todd Wood, Zheng Zhang, et al. Comparison of DNA Sequences with Protein Sequences. Genomics, 1997. https://doi.org/10.1006/geno.1997.4995 ↩

6. William R. Pearson. BLAST and FASTA Similarity Searching for Multiple Sequence Alignment. Methods in Molecular Biology, 2013. https://doi.org/10.1007/978-1-62703-646-7_5 ↩

7. Shiyang Long, Qingting Wei, Yanni Zou, et al. FastqCA: an effective FASTQ compressor through 2D spatial redundancy reduction. Frontiers in Genetics, 2026. https://doi.org/10.3389/fgene.2026.1827471 ↩

8. Meredith L. Carpenter, Jason D. Buenrostro, Cristina Valdiosera, et al. Pulling out the 1%: Whole-Genome Capture for the Targeted Enrichment of Ancient DNA Sequencing Libraries. The American Journal of Human Genetics, 2013. https://doi.org/10.1016/j.ajhg.2013.10.002 ↩