Building a Bioinformatics Pipeline for Your Own Genome and Transcriptome

By the end of this guide you will have a single command that takes your raw sequencing files and produces a set of analysis-ready outputs. The outputs are an indexed CRAM alignment of your genome and a small-variant VCF restricted to high-confidence regions. You also get a gene-level expression count matrix from RNA sequencing, and quality control reports for every stage. The pipeline is written as a Nextflow workflow with a Conda or container environment, so each step is reproducible and resumable. You need a Linux machine with at least 32 GB of RAM and 16 cores, roughly 1 TB of free disk for a 30x human genome plus intermediates, and your FASTQ files. Everything below uses tools you can install in an afternoon. We assume you are comfortable on a command line and reading tool documentation, but not that you have run a variant caller before.

Two notes on scope before the steps. This pipeline measures and quantifies. It does not interpret clinically. Any variant you find that looks medically meaningful belongs in front of a genetic counselor or physician who can order a validated confirmatory test. Research-grade sequencing is not a diagnostic assay.

1. Set up a workflow manager instead of a shell script

The temptation is to write run_everything.sh. Resist it. A ten-step pipeline on human-scale data will fail partway through, and you will want to resume from the failure rather than repeat six hours of alignment. Workflow managers give you caching, parallelism, and provenance. They also make a pipeline portable across a laptop and a cloud batch queue without rewriting the logic. Surveys of the field consistently find that these systems are what separate a pipeline you can rerun from a pipeline you can only run once.¹

We recommend Nextflow. It has the largest body of existing human-genomics workflows (nf-core), handles CRAM and BAM streaming well, and its -resume flag is reliable in practice. Snakemake is a reasonable alternative if you prefer Python and file-pattern rules. The broader point is that the choice of system matters less than using one at all, and reuse of existing published workflows is where most of the time savings live.²

Install Nextflow and a container runtime:

curl -s https://get.nextflow.io | bash
sudo mv nextflow /usr/local/bin/
nextflow -version

# Docker or, on a shared machine without root, Apptainer
sudo apt-get install -y apptainer

Create the project skeleton:

mkdir -p ~/pmp/{ref,fastq,work,results,conf}
cd ~/pmp

Put a nextflow.config in place that pins resources and enables containers:

process {
  executor = 'local'
  container = 'quay.io/biocontainers/mulled-v2-...'   // per-process overrides below
  cpus = 8
  memory = 16.GB
  errorStrategy = 'retry'
  maxRetries = 1
}
apptainer.enabled = true
apptainer.autoMounts = true
workDir = '/home/you/pmp/work'

2. Get the reference genome right

Reference choice affects every downstream number, and a mistake here is expensive to discover late. Use GRCh38 with the analysis set that includes decoy and ALT contigs, and use a version that masks the false duplication in the segmental region around GRCh38 chr21. The no_alt_analysis_set from the Broad or the 1000 Genomes FTP is the practical choice for short-read alignment because it avoids reads scattering across ALT haplotypes.

cd ref
wget https://storage.googleapis.com/genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.fasta
wget https://storage.googleapis.com/genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.fasta.fai
wget https://storage.googleapis.com/genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.dict

samtools faidx Homo_sapiens_assembly38.fasta
bwa-mem2 index Homo_sapiens_assembly38.fasta     # ~80 GB RAM peak, or use bwa index

Download the accessory files you will need for variant filtering and for masking problem regions:

# Genome in a Bottle high-confidence regions, for sanity-checking your calls
wget https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/NA12878_HG001/latest/GRCh38/HG001_GRCh38_1_22_v4.2.1_benchmark.bed
# ENCODE blacklist, useful for RNA-seq coverage work
wget https://github.com/Boyle-Lab/Blacklist/raw/master/lists/hg38-blacklist.v2.bed.gz

Keep a ref/PROVENANCE.md with the exact URLs and download dates. Six months from now you will need to know which GRCh38 patch your VCF coordinates refer to.

3. Quality control the raw reads before you align anything

Alignment is the most expensive step, so spend two minutes checking the input first. fastp is our default because it does adapter detection, quality trimming, and reporting in one pass at high throughput.

fastp \
  -i fastq/WGS_R1.fastq.gz -I fastq/WGS_R2.fastq.gz \
  -o trimmed/WGS_R1.fq.gz -O trimmed/WGS_R2.fq.gz \
  --detect_adapter_for_pe \
  --cut_tail --cut_tail_mean_quality 20 \
  --length_required 50 \
  --thread 8 \
  --json qc/wgs_fastp.json --html qc/wgs_fastp.html

For whole-genome sequencing, be conservative with trimming. Aggressive quality trimming biases coverage and can hurt variant calling more than it helps, because bwa soft-clips low-quality ends anyway. Adapter removal and a minimum length filter are enough. For RNA sequencing, trimming matters more because adapter read-through is common with short inserts.

Numbers to check in the JSON output: duplication rate above roughly 20% on a PCR-free library suggests a low-complexity library, GC content far from 41% for human genome suggests contamination, and a sharp quality drop in the last 20 cycles is normal but should not extend into the middle of the read.

4. Align, mark duplicates, and write CRAM

Align with bwa-mem2, sort, and mark duplicates in a single piped pass so you never write an unsorted BAM to disk. Note the read group: variant callers need it, and getting it wrong causes confusing downstream failures.

RG='@RG\tID:WGS1\tSM:SUBJECT01\tLB:lib1\tPL:ILLUMINA'

bwa-mem2 mem -t 16 -R "$RG" -K 100000000 -Y \
    ref/Homo_sapiens_assembly38.fasta \
    trimmed/WGS_R1.fq.gz trimmed/WGS_R2.fq.gz \
  | samtools fixmate -m -u -@ 4 - - \
  | samtools sort -u -@ 4 -m 2G -T /tmp/srt - \
  | samtools markdup -@ 4 --reference ref/Homo_sapiens_assembly38.fasta \
      -f qc/markdup_stats.txt - results/SUBJECT01.cram

samtools index results/SUBJECT01.cram

Two flags worth understanding. -K 100000000 fixes the batch size so that output is deterministic regardless of thread count, which matters if you want bit-identical reruns. -Y uses soft clipping for supplementary alignments, which keeps sequence in the record and makes structural variant callers happier.

CRAM instead of BAM is a deliberate choice. A 30x human genome is roughly 90 GB as BAM and 45 GB as CRAM at default settings. The cost is that every tool reading the CRAM needs the exact same reference FASTA, so keep it next to the data.

Check coverage before moving on:

samtools coverage results/SUBJECT01.cram | head -25
mosdepth --by 1000 --fast-mode -t 4 qc/SUBJECT01 results/SUBJECT01.cram

For a 30x whole genome you should see mean depth between 28 and 35 with at least 90% of the callable genome above 20x. If mean depth is 30 but the 20x fraction is 70%, coverage is uneven and you likely have a GC bias problem from PCR amplification.

5. Call small variants

For a single sample with good coverage, use DeepVariant. It is a convolutional network over pileup images and it outperforms hand-tuned callers on Illumina short reads without any site-level filtering tuning on your part.

apptainer exec docker://google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
    --model_type=WGS \
    --ref=ref/Homo_sapiens_assembly38.fasta \
    --reads=results/SUBJECT01.cram \
    --output_vcf=results/SUBJECT01.dv.vcf.gz \
    --output_gvcf=results/SUBJECT01.dv.g.vcf.gz \
    --num_shards=16 \
    --intermediate_results_dir=/tmp/dv_tmp

Expect roughly 4 to 8 hours on 16 cores for a 30x genome. The gVCF output is worth generating even for a single sample, because it records confident reference calls. That distinction matters later when you want to say “this position is reference” rather than “this position was not called”.

The GATK path (HaplotypeCaller followed by variant quality score recalibration) remains the reference implementation and the one most exome protocols describe step by step.³ It is useful to know, and it is the right choice if you need joint calling across a cohort. For one person’s genome, DeepVariant gives you fewer knobs to get wrong.

Sanity-check the result:

bcftools stats results/SUBJECT01.dv.vcf.gz > qc/variant_stats.txt
grep -E "^SN|^TSTV" qc/variant_stats.txt

For a single human genome against GRCh38, expect roughly 4.0 to 5.0 million variants in total. About 3.4 to 4.0 million of those are SNVs, and the transition-to-transversion ratio should sit near 2.0 to 2.1. A Ti/Tv of 1.5 means you are calling a lot of false positives. A variant count of 6 million means either contamination or a filtering problem.

6. Annotate

Annotation converts coordinates into biology. Use Ensembl VEP with a local cache, because online mode is slow and rate-limited.

vep --cache --dir_cache ref/vep --offline \
  --assembly GRCh38 --species homo_sapiens \
  --fasta ref/Homo_sapiens_assembly38.fasta \
  --input_file results/SUBJECT01.dv.vcf.gz \
  --output_file results/SUBJECT01.vep.vcf.gz \
  --vcf --compress_output bgzip \
  --everything --pick_allele_gene \
  --plugin CADD,ref/cadd/whole_genome_SNVs.tsv.gz \
  --custom ref/gnomad/gnomad.genomes.v4.1.sites.vcf.bgz,gnomADg,vcf,exact,0,AF,AF_popmax \
  --fork 8 --stats_file qc/vep_summary.html

--pick_allele_gene gives one consequence per allele per gene, which keeps the output tractable. Without it, a single variant can produce twenty transcript-level rows and your downstream filtering becomes ambiguous. If you want the full transcript set for a specific gene, rerun VEP on a subset.

Then reduce to something you can read. A reasonable first pass filters on population frequency and predicted impact:

bcftools view -f PASS results/SUBJECT01.vep.vcf.gz \
  | filter_vep --filter "gnomADg_AF < 0.001 or not gnomADg_AF" \
               --filter "IMPACT is HIGH" \
  -o results/rare_high_impact.vcf

You will get a few hundred to a couple of thousand rows. Most are artifacts, misannotated transcripts, or benign loss-of-function in genes that tolerate it. Do not read this list as a health report. If a variant in it sits in a gene with an established clinical association, that is a reason to talk to a clinical geneticist about targeted confirmatory testing. The variant is not a finding in itself.

7. Quantify RNA expression

Run RNA sequencing quantification as a separate branch of the same workflow. For gene and transcript-level abundance, we use salmon in selective-alignment mode against a decoy-aware transcriptome index. It is fast, corrects for GC and positional bias, and avoids the cost of a genomic alignment you may not need.

# Build the decoy-aware index once
grep "^>" ref/Homo_sapiens_assembly38.fasta | cut -d " " -f 1 | sed 's/>//g' > ref/decoys.txt
cat gencode.v45.transcripts.fa ref/Homo_sapiens_assembly38.fasta > ref/gentrome.fa
salmon index -t ref/gentrome.fa -d ref/decoys.txt -p 16 -i ref/salmon_idx --gencode

salmon quant -i ref/salmon_idx -l A \
  -1 trimmed/RNA_R1.fq.gz -2 trimmed/RNA_R2.fq.gz \
  --validateMappings --seqBias --gcBias --posBias \
  --numBootstraps 30 \
  -p 16 -o results/rna_quant

The bootstraps matter. They give you per-transcript inferential variance, which tells you which estimates are trustworthy for lowly expressed or highly similar transcripts. Check the mapping rate in results/rna_quant/logs/salmon_quant.log. For a good poly-A library you want above 80%. Below 60% usually means rRNA contamination, genomic DNA carryover, or the wrong library type.

If you also want a genomic alignment for splice-junction and variant work in RNA, run STAR with two-pass mode:

STAR --genomeDir ref/star_idx \
  --readFilesIn trimmed/RNA_R1.fq.gz trimmed/RNA_R2.fq.gz \
  --readFilesCommand zcat --runThreadN 16 \
  --twopassMode Basic \
  --outSAMtype BAM SortedByCoordinate \
  --quantMode GeneCounts \
  --outFileNamePrefix results/rna_star_

Aggregate transcript estimates to genes in R with tximport, which handles the effective-length correction properly:

library(tximport); library(readr)
t2g <- read_tsv("ref/tx2gene.tsv", col_names = c("tx","gene"))
txi <- tximport("results/rna_quant/quant.sf", type = "salmon",
                tx2gene = t2g, countsFromAbundance = "lengthScaledTPM")
saveRDS(txi, "results/gene_level.rds")

A single RNA sample has no comparison group, so treat the output as a baseline. Its value is longitudinal. The same tissue with the same library prep, sampled again in six months, gives you a within-person contrast that is far more informative than comparing yourself to a public cohort.

8. Wrap it, log it, and make it rerunnable

Now convert the commands into Nextflow processes. The pattern is one process per tool. Each one has explicit input and output, plus a publishDir for anything you want to keep.

process ALIGN {
  tag "${sample}"
  cpus 16
  memory 64.GB
  container 'quay.io/biocontainers/bwa-mem2:2.2.1--he513fc3_0'
  publishDir "results/align", mode: 'copy'

  input:
    tuple val(sample), path(r1), path(r2)
    path ref
  output:
    tuple val(sample), path("${sample}.cram"), path("${sample}.cram.crai")
  script:
  """
  bwa-mem2 mem -t ${task.cpus} -K 100000000 -Y \\
      -R '@RG\\tID:${sample}\\tSM:${sample}\\tPL:ILLUMINA' \\
      ${ref}/Homo_sapiens_assembly38.fasta ${r1} ${r2} \\
    | samtools fixmate -m -u - - \\
    | samtools sort -u -m 2G -T tmp - \\
    | samtools markdup --reference ${ref}/Homo_sapiens_assembly38.fasta - ${sample}.cram
  samtools index ${sample}.cram
  """
}

Run it and keep the execution record:

nextflow run main.nf -profile apptainer \
  -with-report results/report.html \
  -with-trace results/trace.txt \
  -with-dag results/dag.svg \
  -resume

The trace file gives you real runtime and peak memory per process, which is how you size resources on the second run instead of guessing. Pin every container to a digest, not a tag, and commit the pipeline to git with the reference provenance file. Reproducibility comes from this combination of version-controlled workflow code and pinned execution environments rather than from documentation alone.⁴

Finally, collect all the QC into one page:

multiqc qc/ results/ -o results/multiqc

Common problems

Read group mismatches are the most common early failure. If DeepVariant or GATK complains about a missing sample name, your -R string was malformed, usually because the shell interpreted \t literally. Check with samtools view -H file.cram | grep '^@RG' and confirm you see real tab characters.

CRAM reference errors appear as garbled sequence or outright decode failures. This happens when a downstream tool cannot find the exact FASTA used at compression time. Set REF_PATH and REF_CACHE environment variables, or pass --reference explicitly to every samtools invocation.

Silent contig-order problems are worse than loud failures. If your VCF has contigs in a different order from the reference dictionary, some tools will process it and produce nonsense. Run bcftools view -h and compare against the .dict file order. Rebuild with bcftools reheader if they disagree.

Out-of-memory kills during bwa-mem2 indexing catch most people once. Indexing GRCh38 peaks near 80 GB. If your machine cannot do that, use classic bwa index (about 5 GB peak, slower alignment) or download a prebuilt index.

Low RNA mapping rates with a high salmon “mapping rate” discrepancy usually mean the wrong library type. Let -l A detect it, then check the inferred type in lib_format_counts.json. If it reports ISF when you expected ISR, your R1 and R2 are swapped.

Pipelines that only run on your machine are a recurring problem across the field, and the fix is the same everywhere: containerize each step and let the workflow manager handle staging and scheduling.¹ Reuse published workflows where they exist, because most pipeline bugs have already been found and fixed by someone else.²

Footnotes

1. Laura Wratten, Andreas Wilm, Jonathan Göke. Reproducible, scalable, and shareable analysis pipelines with bioinformatics workflow managers. Nature Methods, 2021. https://doi.org/10.1038/s41592-021-01254-9 ↩ ↩²

2. Marine Djaffardjy, George Marchment, Clémence Sebe, et al. Developing and reusing bioinformatics data analysis pipelines using scientific workflow systems. Computational and Structural Biotechnology Journal, 2023. https://doi.org/10.1016/j.csbj.2023.03.003 ↩ ↩²

3. Narendra Meena, Praveen Mathur, Krishna Medicherla, et al. A Bioinformatics Pipeline for Whole Exome Sequencing: Overview of the Processing and Steps from Raw Data to Downstream Analysis. BIO-PROTOCOL, 2018. https://doi.org/10.21769/bioprotoc.2805 ↩

4. Neha Kulkarni, Luca Alessandrì, Riccardo Panero, et al. Reproducible bioinformatics project: a community for reproducible bioinformatics analysis pipelines. BMC Bioinformatics, 2018. https://doi.org/10.1186/s12859-018-2296-x ↩