Bioinformatics Protocols
Battle-tested, copy-paste protocols. Every command is verified, with expected output and fixes for the errors you will actually hit.
Count the number of reads in a FASTQ file
A one-liner to count sequencing reads in a FASTQ file, why it works, and how to handle gzipped input.
Run a local BLAST search against a custom database
Build a BLAST database from a FASTA file and query it locally with blastn — no web submission required.
Generate a FastQC quality report for raw reads
Run FastQC on raw FASTQ reads to produce an HTML quality report, then optionally merge many reports into one summary with MultiQC.
Trim adapters and low-quality bases with fastp
Quality- and adapter-trim paired- or single-end FASTQ reads with fastp, producing gzipped outputs plus HTML and JSON QC reports.
Randomly subsample a FASTQ file with seqtk
Draw a reproducible random subset of reads from a FASTQ file with seqtk sample, keeping paired mates in sync via a shared seed.
Get sequence statistics for FASTA/FASTQ with seqkit
Use seqkit stats to print read and contig counts, length distribution, N50, GC, and quality metrics for any FASTA or FASTQ file.
Convert FASTQ to FASTA
Convert FASTQ reads to FASTA with seqkit fq2fa, plus a dependency-free awk one-liner that streams gzipped input.
Extract specific reads from a FASTQ by ID
Pull specific reads out of a FASTQ by matching their IDs against a list with seqkit grep, or invert the match to drop them.
Align short reads to a reference with BWA-MEM
Index a reference genome, align paired-end short reads with BWA-MEM, and produce a coordinate-sorted, indexed BAM in one pipe.
Sort and index a BAM file with samtools
Coordinate-sort a BAM file with samtools and build its .bai index so tools can randomly access any genomic region.
Summarize alignment stats with samtools flagstat
Read mapping rate, duplicate counts, and per-chromosome read totals from a BAM in seconds with samtools flagstat and idxstats.
Extract reads from a genomic region in a BAM
Pull every read overlapping a genomic locus out of an indexed BAM into a smaller region BAM with samtools view, then count reads in the interval.
Mark PCR duplicates with samtools markdup
Flag PCR duplicates in an aligned BAM using the four-step samtools collate, fixmate -m, sort, and markdup pipeline.
Align RNA-seq reads with HISAT2
Build a HISAT2 splice-aware index and align paired RNA-seq reads to a genome, piping straight into samtools for a sorted BAM.
Call variants with bcftools mpileup and call
Pile up read alignments and call SNPs and indels from a sorted BAM with bcftools, producing a compressed, indexed VCF.
Filter a VCF by depth and quality with bcftools
Drop low-confidence variants from a VCF by hard-filtering on QUAL and INFO/DP with bcftools, then compare record counts before and after.
Compress and index a VCF with bgzip and tabix
Block-compress a VCF with bgzip and build a tabix index so you can pull variants from any genomic region in milliseconds.
Find overlapping intervals with bedtools intersect
Use bedtools intersect to report, count, or exclude overlapping intervals between two BED/GFF/VCF files.
Assemble a bacterial genome with SPAdes
De novo assemble a bacterial isolate from paired Illumina reads with SPAdes, then sanity-check contig count, total length, and N50 using seqkit.
Assess assembly quality with QUAST
Run QUAST on a genome assembly to compute N50, total length, and contig counts, then add reference-based misassembly and genome-fraction metrics.
Check genome completeness with BUSCO
Run BUSCO against a lineage ortholog set to score genome assembly completeness, then read the C/S/D/F/M breakdown from the short summary.
Quantify transcript expression with Salmon
Build a Salmon transcriptome index and quantify paired-end RNA-seq reads alignment-free, with automatic library-type detection.
Align sequences with MAFFT
Build a multiple sequence alignment from a multi-FASTA with MAFFT, from the auto strategy to speed and accuracy tuning.
Build a maximum-likelihood tree with IQ-TREE
Infer a maximum-likelihood phylogeny from an alignment with IQ-TREE 2, using ModelFinder for model selection and 1000 ultrafast bootstraps.
Trim a multiple sequence alignment with trimAl
Strip gappy, poorly-aligned columns from a multiple sequence alignment with trimAl before you infer a phylogenetic tree.
Index a FASTA and extract a region with samtools faidx
Build a .fai index for a FASTA file with samtools faidx, then pull any sub-sequence by chromosome and coordinate.
Extract sequences from coordinates with bedtools getfasta
Turn BED intervals into FASTA sequences with bedtools getfasta, including strand-aware extraction that reverse-complements minus-strand features.
Download sequencing data from the SRA
Fetch public sequencing reads from the NCBI SRA by accession with prefetch and fasterq-dump, then split mates and compress to gzipped FASTQ.
Download a reference genome with NCBI datasets
Pull any reference genome assembly from NCBI by accession or taxon with the datasets CLI, then unzip it to a ready-to-use FASTA.
Calculate GC content of sequences with seqkit
Report per-sequence GC% and the overall mean GC of a FASTA/FASTQ file using seqkit fx2tab plus a short awk one-liner.