The open source bioinformatics stack for analyzing your own molecular data

If you have your own sequencing data and want to work with it, the short answer is that almost everything you need is free, command-line, and installable in an afternoon. The tools are the following.

• samtools/bcftools for alignment and variant files
• bwa-mem2 or minimap2 for alignment
• DeepVariant or GATK HaplotypeCaller for small-variant calling
• Ensembl VEP for annotation
• salmon plus DESeq2 or limma for expression
• PLINK 2 for anything genotype-wide
• Nextflow or Snakemake with a conda or container environment to hold it together

The hard part is choosing among the several hundred tools that do nominally the same thing, and avoiding the handful of silent failure modes that produce plausible-looking output from wrong inputs. Finding tools is not the difficulty.

This page describes the stack we would run on a single person’s data, with the parameters that matter and the reasons behind each choice.

What “open source” buys you, and where the term gets fuzzy

The practical benefit of open source in this field is that you can read the code that produced a number you are about to act on. The ecosystem is unusually healthy by the normal measures of open development: package counts, release cadence, issue traffic, and contributor breadth have been used as proxies for real-world impact of bioinformatics software, and the distribution is heavily long-tailed, with a small number of tools carrying most of the usage.¹ A large fraction of the working set is redistributed as containers through BioContainers, which builds Docker and Singularity images directly from Bioconda recipes so that a tool and its dependency graph travel together.²

Be careful with licenses in proteomics and mass spectrometry, where “free” often means free-for-academic-use binaries rather than open source. MSFragger and DIA-NN are distributed as binaries under their own terms. MaxQuant is free but closed. OpenMS, ThermoRawFileParser, and msconvert from ProteoWizard are genuinely open. If you intend to build anything on top of your pipeline, read the license before you design around a tool.

Whole-genome sequencing: the path from FASTQ to an annotated VCF

Start by deciding on a reference build and never deviate from it. GRCh38 with the analysis set that includes decoy and alt contigs is the pragmatic default, because gnomAD, ClinVar, and nearly every annotation source are keyed to it. T2T-CHM13 is a better assembly and resolves regions GRCh38 cannot, but the annotation ecosystem lags. Mixing builds is the single most common way to get confidently wrong answers, since coordinates differ by millions of bases and nothing in the file format will warn you.

For alignment, bwa-mem2 mem -t 32 -R '@RG\tID:1\tSM:sample\tPL:ILLUMINA' is the standard for short reads, piped straight into samtools sort to avoid writing an intermediate SAM. Store the result as CRAM, not BAM: a 30x human genome is roughly 80–90 GB as BAM and closer to 15–25 GB as CRAM. CRAM’s compression is reference-based, which means the file is unreadable without the exact reference FASTA it was written against. Keep the reference and its .fai index alongside the CRAM, and record the MD5 of the FASTA.

For small variants we prefer DeepVariant over GATK for germline short-read data. Run it with --model_type=WGS and the matching model for your platform. It is a single step, it requires no base quality score recalibration, and it tends to produce fewer false-positive indels in homopolymer regions. GATK remains the right choice when you need joint calling across samples or the full gVCF workflow. Either way, expect roughly 4–5 million variants for one human genome against GRCh38, and normalize before you do anything else:

bcftools norm -m -any -f GRCh38.fa -Ou in.vcf.gz \
  | bcftools norm -d exact -Oz -o normalized.vcf.gz

Left-alignment and splitting of multi-allelic sites is the step people skip, and skipping it causes annotation joins against ClinVar or gnomAD to miss real matches because the same indel is represented two different ways.

Structural variants need separate callers. Manta is fast and accurate for deletions and duplications, Delly catches a different slice, and smoove wraps Lumpy with sensible filtering. Short reads miss a great deal of structural variation regardless, so treat a negative result as uninformative rather than reassuring. ExpansionHunter handles known repeat expansions given a variant catalog.

Annotate with Ensembl VEP in offline mode: vep --cache --offline --assembly GRCh38 --fork 8 --everything --vcf. Be deliberate about --pick and --pick_order. VEP reports a consequence per transcript, and the single “picked” consequence can change from missense to intronic depending on which transcript wins. For a personal genome, we keep all consequences and filter downstream rather than letting a flag make the call.

Anything that looks clinically meaningful in ClinVar or an ACMG secondary-findings gene should be confirmed in a CLIA-certified laboratory and interpreted with a genetic counselor or physician. Research-grade calls, including ours, are not diagnostic, and pathogenicity classification depends on family history and phenotype that a pipeline does not have.

RNA-seq and the problem of n = 1

For expression, salmon with a decoy-aware index built from the genome is the fastest reasonable choice: salmon quant -l A --gcBias --seqBias --validateMappings -p 16. Import to gene level with tximport, then use DESeq2 or limma-voom. If you plan to look at splicing or novel transcripts, you need genome alignment instead, so use STAR with two-pass mode.

The statistical reality of a single person’s longitudinal data is that standard differential expression models assume replication across samples. With one subject sampled over time, treat each timepoint as a repeated measure, model the subject-specific baseline, and report effect sizes with uncertainty rather than p-values borrowed from a group-comparison design. Whole-blood RNA is also dominated by cell composition, so run a deconvolution method such as granulator or ABIS before concluding that a gene changed. A shift in neutrophil fraction will move thousands of genes at once.

Holding the pipeline together

Install with mamba from the bioconda and conda-forge channels, pinned in an environment.yml with explicit versions. For multi-step work, use Nextflow with an nf-core pipeline (nf-core/sarek for germline variant calling, nf-core/rnaseq for expression) and -profile docker or -profile singularity so every process runs in a pinned container. Snakemake is the better choice if you prefer Python and file-target semantics. The point of either is that the workflow engine records what ran with which container digest, which is the only way to answer “why did this number change” six months later. Interoperating tools across environments has been a recognized problem since well before containers were standard, and the current answer is to make the environment part of the artifact.³²

If you write your own tooling, Rust is increasingly viable for the hot paths: rust-bio provides alignment, pattern matching, and format parsers with performance close to C and memory safety by construction.⁴

Questions people also ask

Is AI going to replace bioinformatics?

It is changing the interface rather than eliminating the work. Large language models are now good at writing the bcftools invocation, explaining a VEP consequence term, and scaffolding an analysis script, and the field’s history has been a progression of interfaces from command line to web portals to conversational agents.⁵ Deep learning has also become the method inside key tools, as in DeepVariant’s treatment of variant calling as image classification, and machine learning is now central to variant effect prediction and multi-omic integration.⁶ What does not go away is judgment about reference builds, batch effects, cell composition, and whether a result is real.

Is Python or R better for bioinformatics?

Use both, for different jobs. Python owns pipeline orchestration, file wrangling, single-cell analysis through scanpy and anndata, and anything touching machine learning. R owns statistical modeling of omics data, where Bioconductor’s DESeq2, limma, edgeR, and the associated annotation packages have no real Python equivalent in maturity. A reasonable split is Python for everything upstream of the count matrix, R from the count matrix onward.

What are the “big three” in bioinformatics?

The phrase usually refers to the three international sequence databases that mirror each other daily: GenBank at NCBI, the European Nucleotide Archive at EMBL-EBI, and the DNA Data Bank of Japan, which together form the INSDC. In a tooling context people sometimes mean the three foundational open projects instead: Biopython, Bioconductor, and the samtools/htslib family.

What is the highest paying job in bioinformatics?

Industry roles pay well above academic ones, and the top of the range is typically computational biology or machine learning positions at pharmaceutical and biotech companies, followed by quantitative roles at hedge funds and diagnostics startups that need genomic data engineering. The differentiator is usually strong software engineering paired with genuine domain knowledge rather than either alone.

Footnotes

1. Mikhail G. Dozmorov. GitHub Statistics as a Measure of the Impact of Open-Source Bioinformatics Software. Frontiers in Bioengineering and Biotechnology, 2018. https://doi.org/10.3389/fbioe.2018.00198 ↩

2. Jingwen Bai, Chakradhar Bandla, Jiaxin Guo, et al. BioContainers Registry: Searching Bioinformatics and Proteomics Tools, Packages, and Containers. Journal of Proteome Research, 2021. https://doi.org/10.1021/acs.jproteome.0c00904 ↩ ↩²

3. P. Romano, R. Giugno, A. Pulvirenti. Tools and collaborative environments for bioinformatics research. Briefings in Bioinformatics, 2011. https://doi.org/10.1093/bib/bbr055 ↩

4. Johannes Köster. Rust-Bio: a fast and safe bioinformatics library. Bioinformatics, 2015. https://doi.org/10.1093/bioinformatics/btv573 ↩

5. Van Q Truong, Marylyn D Ritchie. Eras of bioinformatics technologies from command-line interfaces to artificial intelligence (AI) chatbots. Briefings in Bioinformatics, 2026. https://doi.org/10.1093/bib/bbag256 ↩

6. David B. Olawade, Ayomikun Kade, Eghosasere Egbon, et al. Bioinformatics and artificial intelligence in genomic data analysis: current advances and future directions. Molecular Genetics and Genomics, 2025. https://doi.org/10.1007/s00438-025-02314-x ↩