Which Bioinformatics Tools You Need for Your Own Data

If you have your own whole-genome FASTQs, RNA-seq reads, a proteomics matrix, and a few years of labs, the tool list you need is short.

• bwa-mem2 or dragen-os for alignment
• samtools for everything BAM-shaped
• DeepVariant for small variants
• manta plus CNVkit for structural variants and copy number
• bcftools and VEP for annotation
• salmon or STAR+featureCounts for RNA quantification
• DESeq2 or limma-voom for expression comparisons
• nextflow to hold it together

Everything else in the standard “top 50 bioinformatics tools” listicle is either a specialized instrument for a problem you don’t have or a web front end to one of the above.

The hard part is reference genome choices, coordinate conventions, and knowing which outputs are trustworthy at which coverage. Tool selection is no part of it. That is what the rest of this covers.

Start with the file formats, because they constrain everything

Raw reads arrive as gzipped FASTQ (.fq.gz), typically paired: sample_R1.fastq.gz, sample_R2.fastq.gz. A 30× human WGS run is roughly 100 GB in FASTQ, around 45–55 GB as a CRAM. Check md5sum against the manifest before you do anything else. Truncated transfers are the single most common failure and they produce alignments that look fine until coverage drops off a cliff in the last chromosomes processed.

Aligned reads live in BAM or CRAM. CRAM is reference-compressed, so a CRAM without the exact reference FASTA it was written against is unreadable. Keep the reference next to the CRAM, and keep its .fai and .dict. Variants live in VCF or its binary form BCF. Expression counts live as TSV or in an R SummarizedExperiment. Proteomics from most platforms arrives as a wide matrix: protein or aptamer identifiers by sample, plus a normalization annotation you need to read carefully.

The one convention that causes the most silent damage: BED is 0-based half-open, VCF and GFF are 1-based inclusive. A one-base shift will not error, it will just give you a different gene.

Reference genome: GRCh38, with the analysis set

Use GRCh38 and specifically the “no-alt analysis set” unless you have a caller that is alt-aware. Alt contigs split reads across loci and drive mapping quality to zero in regions like the MHC, which quietly deletes variants from your call set. T2T-CHM13 is a better assembly and resolves centromeres and segmental duplications that GRCh38 gets wrong. Almost every annotation database you want to join against is nonetheless keyed to GRCh38 coordinates. That covers ClinVar, gnomAD, dbSNP builds and most GWAS summary statistics. Our recommendation: call on GRCh38 as the primary, and run a second pass on T2T only if you care about a specific region it improves.

Do not mix. A VCF lifted over from hg19 with CrossMap or picard LiftoverVcf will lose 0.1–1% of sites and silently flip reference and alternate alleles at some others.

The alignment and small-variant path

bwa-mem2 mem -t 32 -R '@RG\tID:s1\tSM:s1\tPL:ILLUMINA' \
  GRCh38_no_alt.fa s1_R1.fq.gz s1_R2.fq.gz \
  | samtools sort -@ 8 -o s1.bam -
samtools index s1.bam

Then duplicates. For PCR-free libraries, marking duplicates changes very little, but mark them anyway with samtools markdup or gatk MarkDuplicates so downstream callers behave predictably.

For small variants (SNVs and indels), we use DeepVariant. It is a convolutional model over pileup images, it needs no base quality score recalibration, and on 30× Illumina WGS it reaches SNV F1 above 0.999 against GIAB benchmarks. GATK HaplotypeCaller is the alternative and is better documented for joint calling across cohorts, which is not your situation with n=1. Run DeepVariant through its Docker image with --model_type=WGS and give it --num_shards equal to your core count.

Validate before you interpret anything. Download the GIAB HG002 truth set and run hap.py against the high-confidence BED. If your pipeline reproduces the published numbers on HG002, it is working. If it does not, you have a reference or a region-filter problem, and no amount of downstream annotation will fix it. The validation step is the difference between a pipeline and a guess, and clinical genetics workflows treat it as mandatory for exactly that reason.¹

Coverage matters more than caller choice. At 30×, heterozygous SNV sensitivity is high. Indels longer than about 20 bp in homopolymer regions are unreliable on short reads at any depth. If a result hinges on an indel in a repeat, it needs orthogonal confirmation, and any variant you plan to act on needs confirmation in a CLIA laboratory and a conversation with a genetic counselor or clinician. Research-grade calls are not clinical results.

Structural variants and copy number

Short-read SV calling is the weakest part of the stack. manta is fast and gives good deletion and duplication calls with breakpoint resolution. CNVkit in whole-genome mode gives copy number segments but needs a panel of normals or a flat reference to behave. Expect a high false positive rate in segmental duplications. Intersect calls with a repeat mask and a gnomAD-SV frequency filter before you look at anything.

RNA-seq

Two paths, and the choice depends on what you want.

Selective alignment with salmon against a decoy-aware transcriptome index is fast (roughly 10 minutes for 30M reads on 16 cores) and gives transcript-level TPMs. Build the index with the full genome as decoy, otherwise reads from unannotated regions get force-assigned to transcripts.

salmon index -t gentrome.fa.gz -d decoys.txt -i salmon_idx -k 31
salmon quant -i salmon_idx -l A -1 r1.fq.gz -2 r2.fq.gz \
  --validateMappings --gcBias --seqBias -p 16 -o quant

STAR plus featureCounts gives you a genomic BAM, which you need for splice junction analysis, allele-specific expression, or looking at coverage over a specific exon in IGV. It costs about 32 GB of RAM for the human index.

For n=1 longitudinal data, differential expression is a different problem than the cohort case DESeq2 was built for. With repeated samples from one person over time, you have a time series. Fit limma-voom with a design matrix encoding time or condition, or use DESeq2 with your own timepoints as replicates. What you cannot do is compare one sample against a public reference cohort and call the differences meaningful. Batch effects between library preps will dominate.

Annotation, and where AI enters

Annotate with Ensembl VEP or SnpEff, then join against ClinVar, gnomAD v4, and dbNSFP. VEP’s --pick flag collapses to one consequence per variant, which is convenient and lossy. We prefer to keep all consequences and filter later.

The machine learning in bioinformatics is mostly inside the tools, not around them: DeepVariant’s classifier, AlphaFold and ESMFold for structure, AlphaMissense and CADD and REVEL for variant effect prediction. This is a long-running direction rather than a new one, with machine learning applied to sequence classification and structure prediction for decades.² The recent expansion is in scale and in applications like drug repurposing built on the same representations.³ Language models are genuinely useful as an interface to your own data: pointing an agent at a filtered VCF and asking it to summarize consequence classes works well. Asking it to interpret a variant’s clinical significance does not, because it will confabulate frequencies and inheritance patterns. Keep the model on the mechanical work and the databases on the facts.

Orchestration

Run everything under Nextflow or Snakemake with containers pinned by digest. nf-core/sarek covers the germline WGS path and nf-core/rnaseq covers RNA. Both are Nextflow pipelines with published test data and version-pinned containers. You will spend an afternoon learning the config system and save weeks of irreproducibility. The open-source toolkit ecosystem is what makes this possible at all, and its value has always been in composability rather than in any single program.⁴ The EBI’s hosted tool framework is a reasonable fallback for one-off jobs like a quick alignment or a BLAST search when you do not want local infrastructure.⁵

Store intermediate outputs. Recomputing a 30× alignment is 6–12 CPU-hours. Keeping the CRAM costs a few dollars a year.

Questions people also ask

What are some common software tools used in bioinformatics? For sequence data, the tools are as follows.

• BWA, Bowtie2, STAR and minimap2 (alignment)
• samtools and bcftools (manipulation)
• GATK and DeepVariant (variant calling)
• VEP and SnpEff (annotation)
• salmon and kallisto (RNA quantification)
• DESeq2 and limma (statistics)
• IGV (visualization)

For orchestration: Nextflow, Snakemake, Docker and Singularity. Next-generation sequencing pushed this toolset into its current shape, together with the databases the tools read from.⁶

Does bioinformatics require AI? No. Alignment, variant calling from pileups, and count-based differential expression are classical algorithms and statistics: Burrows-Wheeler indexing, hidden Markov models, negative binomial regression. Many standard methods long predate the current wave of machine learning.⁷ Neural networks now outperform hand-built models on specific tasks like variant classification and protein structure, so the best tool for some steps happens to be a model.

What are examples of bioinformatics applications? Examples include the following.

• Variant calling and annotation
• transcript quantification
• protein structure prediction
• phylogenetics
• metagenomic classification of microbiome reads
• epitope prediction
• mass-spec peptide identification
• copy number analysis
• methylation calling from bisulfite or nanopore data
• pharmacogenomic genotyping

The same core methods carry across genomics, proteomics, agriculture, and microbiology.⁸

Are there tools that do this without the command line? Yes. Galaxy, Geneious Prime, and web platforms like Sangerbox wrap common analyses in a GUI.⁹ They are fine for exploration. For a longitudinal dataset you intend to reanalyze as tools improve, a scripted pipeline you can rerun with one command is worth the learning cost.

Footnotes

1. Rute Pereira, Jorge Oliveira, Mário Sousa. Bioinformatics and Computational Tools for Next-Generation Sequencing Analysis in Clinical Genetics. Journal of Clinical Medicine, 2020. https://doi.org/10.3390/jcm9010132 ↩

2. Z EZZIANE. Applications of artificial intelligence in bioinformatics: A review. Expert Systems with Applications, 2006. https://doi.org/10.1016/j.eswa.2005.09.042 ↩

3. Unknown. Artificial Intelligence in Bioinformatics and Drug Repurposing: Methods and Applications. Frontiers Research Topics, 2022. https://doi.org/10.3389/978-2-88974-881-5 ↩

4. J. E. Stajich. Open source tools and toolkits for bioinformatics: significance, and where are we?. Briefings in Bioinformatics, 2006. https://doi.org/10.1093/bib/bbl026 ↩

5. M. Goujon, H. McWilliam, W. Li, et al. A new bioinformatics analysis tools framework at EMBL-EBI. Nucleic Acids Research, 2010. https://doi.org/10.1093/nar/gkq313 ↩

6. H. C. Lee, K. Lai, M. T. Lorenc, et al. Bioinformatics tools and databases for analysis of next-generation sequence data. Briefings in Functional Genomics, 2011. https://doi.org/10.1093/bfgp/elr037 ↩

7. IZET M. KAPETANOVIC, SIMON ROSENFELD, GRANT IZMIRLIAN. Overview of Commonly Used Bioinformatics Methods and Their Applications. Annals of the New York Academy of Sciences, 2004. https://doi.org/10.1196/annals.1310.003 ↩

8. Muhammad Aamer Mehmood. Use of Bioinformatics Tools in Different Spheres of Life Sciences. Journal of Data Mining in Genomics & Proteomics, 2014. https://doi.org/10.4172/2153-0602.1000158 ↩

9. Weitao Shen, Ziguang Song, Xiao Zhong, et al. Sangerbox: A comprehensive, interaction‐friendly clinical bioinformatics analysis platform. iMeta, 2022. https://doi.org/10.1002/imt2.36 ↩