Mapping reads

Mapping finds where each read comes from in a reference sequence. minimap2, Bowtie 2, BWA and BWA-MEM2 map DNA reads, HISAT2 maps RNA-seq reads across introns; samtools sorts, indexes, filters and counts the alignments, and mosdepth measures the depth.

The examples work in a folder of your own:

The reference

The example reference has two pieces of human chromosome 10 (GRCh37), named human_CYP2C19 and human_CYP2C9. samtools faidx indexes a FASTA file and cuts regions out of it:

cp /data/examples/reads/reference.fa .
samtools faidx reference.fa
cut -f 1,2 reference.fa.fai
samtools faidx reference.fa human_CYP2C19:11590-11640

minimap2

minimap2 -ax sr reference.fa /data/examples/reads/NA12878_R1.fastq /data/examples/reads/NA12878_R2.fastq > mapped.sam
samtools sort -o mapped.bam mapped.sam
samtools index mapped.bam
samtools flagstat mapped.bam

-ax sr is the preset for short paired reads (map-ont and map-pb for long reads, asm5 for an assembly against a reference). The two steps can be one pipe: minimap2 -ax sr … | samtools sort -o mapped.bam -. A read group names the sample in the BAM file – variant callers use it: -R '@RG\tID:NA12878\tSM:NA12878'.

Bowtie 2

bowtie2-build reference.fa refindex > /dev/null
bowtie2 -x refindex -1 /data/examples/reads/NA12878_R1.fastq -2 /data/examples/reads/NA12878_R2.fastq 2> bowtie2.log | samtools sort -o bt2.bam -
samtools index bt2.bam
cat bowtie2.log

bowtie2-build makes the index (six .bt2 files) once. The summary of the mapping goes to the standard error – here into bowtie2.log. Presets: --very-fast … --very-sensitive, and --local for local alignment.

BWA

bwa index reference.fa
bwa mem -R '@RG\tID:NA12878\tSM:NA12878' reference.fa /data/examples/reads/NA12878_R1.fastq /data/examples/reads/NA12878_R2.fastq 2> bwa.log | samtools sort -o bwa.bam -
samtools index bwa.bam
samtools flagstat bwa.bam

bwa index makes the index once: five files (.amb .ann .bwt .pac .sa) next to the FASTA file – so the FASTA file must be in a folder of your own, not in the read-only examples. bwa mem is BWA's aligner for reads of 70 bases and more; it reads .fastq.gz files too. Its messages go to the standard error (here bwa.log); the CPU times in them are not measured in the browser – the real times are.

BWA-MEM2

bwa-mem2 index reference.fa > /dev/null 2>&1
bwa-mem2 mem -R '@RG\tID:NA12878\tSM:NA12878' reference.fa /data/examples/reads/NA12878_R1.fastq /data/examples/reads/NA12878_R2.fastq 2> bwamem2.log | samtools sort -o bwamem2.bam -
samtools index bwamem2.bam
samtools flagstat bwamem2.bam

BWA-MEM2 is bwa mem's algorithm, rewritten to use the vector instructions of modern processors: on a server it maps several times faster than BWA, with the same alignments as bwa mem 0.7.17. It needs an index of its own – bwa-mem2 index writes .0123 .amb .ann .bwt.2bit.64 .pac; one made by bwa index does not do. In the browser it is not faster than BWA, and it takes a few hundred MB of memory even for a small reference. Compare the two: samtools view bwa.bam | cut -f 1-6 | head and the same for bwamem2.bam.

HISAT2, for RNA-seq reads

A read from a spliced RNA can span an intron: its first part matches one exon, the rest the next. HISAT2 splits such reads across the intron. The example reads are 256 pairs of yeast RNA-seq reads; the reference is chromosome I of yeast.

hisat2-build /data/examples/rnaseq/yeast_chrI.fa yeast_chrI > /dev/null 2>&1
hisat2 -x yeast_chrI -1 /data/examples/rnaseq/SRR6357070_1.fastq -2 /data/examples/rnaseq/SRR6357070_2.fastq 2> hisat2.log | samtools sort -o rnaseq.bam -
samtools index rnaseq.bam
cat hisat2.log
samtools view rnaseq.bam | awk '$6 ~ /N/ {print $3, $4, $6}' | head -3

hisat2-build writes eight .ht2 files; -x names them without the ending. A spliced read has an N in its CIGAR string: 94M113N7M is 94 bases on the genome, an intron of 113 bases skipped, then 7 more. Known introns can help: --known-splicesite-infile /data/examples/rnaseq/yeast_chrI.splice_sites.tsv; for DNA reads use --no-spliced-alignment. hisat2 unpacks .gz read files first, as its wrapper script does; --un and --al (also done by the wrapper) are not available – samtools fastq -f 4 gets the unmapped reads. The RNA-seq recipe goes on to differential expression.

samtools

Command	What it does
`samtools view -H mapped.bam`	the header (reference names, read groups, the programs that made the file)
`samtools view mapped.bam \| head`	alignments as text: name, flag, position, MAPQ, CIGAR …
`samtools view -b -q 20 -F 4 -o good.bam mapped.bam`	keep mapped reads with MAPQ ≥ 20
`samtools view -c mapped.bam human_CYP2C19:11000-12000`	count the reads in a region (needs the index)
`samtools flagstat mapped.bam`	mapped, properly paired, duplicates …
`samtools stats mapped.bam \| grep ^SN`	summary numbers: error rate, insert size, read lengths
`samtools idxstats mapped.bam`	reads per reference sequence
`samtools depth -a -r human_CYP2C19:11610-11620 mapped.bam`	depth at each position
`samtools tview -d T -p human_CYP2C19:11600 mapped.bam reference.fa`	the reads as text around a position

samtools view mapped.bam | head -3 | cut -f 1-9
samtools view -c mapped.bam human_CYP2C19:11000-12000
samtools depth -a -r human_CYP2C19:11610-11620 mapped.bam

To see the reads, igv mapped.bam opens them in the Genome browser.

Depth with mosdepth

samtools depth prints a line for every position; mosdepth summarises the depth of a whole BAM file at once – per sequence, in windows, or in the regions of a BED file – and writes compressed BED files that the genome browser can show. It needs the index of the BAM file (samtools index).

mosdepth cyp2c mapped.bam
cat cyp2c.mosdepth.summary.txt
mosdepth -n --by /data/examples/annotation/cyp2c_exons.bed --thresholds 1,10,20 exons mapped.bam
zcat exons.regions.bed.gz | head -5
zcat exons.thresholds.bed.gz | head -5

The summary has the length of each sequence, the bases covered by reads, and the mean, smallest and largest depth. cyp2c.mosdepth.global.dist.txt tells which share of the bases has at least each depth – the usual "90% of the exome at 20×". -n leaves out the per-base file; --by 1000 gives windows of 1000 bases instead of a BED file; -Q 20 counts only reads with a mapping quality of 20 or more. mosdepth keeps an array as long as the longest sequence in the BAM file: yeast chromosomes and the example slices are fine, a whole human chromosome is too big for a browser.

SAM flags

The second column of a SAM line is a sum of bits: 1 paired, 2 properly paired, 4 unmapped, 16 reverse strand, 256 secondary, 1024 duplicate, 2048 supplementary. -f keeps reads that have all the given bits, -F drops reads that have any of them.

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

minimap2: Li H (2018) Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34(18): 3094–3100. doi:10.1093/bioinformatics/bty191; Li H (2021) New strategies to improve minimap2 alignment accuracy. Bioinformatics 37(23): 4572–4574. doi:10.1093/bioinformatics/btab705
Bowtie 2: Langmead B, Salzberg SL (2012) Fast gapped-read alignment with Bowtie 2. Nature Methods 9(4): 357–359. doi:10.1038/nmeth.1923
BWA (bwa mem): Li H (2013) Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997 [q-bio.GN]. arxiv.org/abs/1303.3997; Li H, Durbin R (2009) Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics 25(14): 1754–1760. doi:10.1093/bioinformatics/btp324
BWA-MEM2: Vasimuddin M, Misra S, Li H, Aluru S (2019) Efficient Architecture-Aware Acceleration of BWA-MEM for Multicore Systems. 2019 IEEE International Parallel and Distributed Processing Symposium (IPDPS), pp. 314–324. doi:10.1109/IPDPS.2019.00041
HISAT2: Kim D, Paggi JM, Park C, Bennett C, Salzberg SL (2019) Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nature Biotechnology 37(8): 907–915. doi:10.1038/s41587-019-0201-4
SAMtools: Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. (2021) Twelve years of SAMtools and BCFtools. GigaScience 10(2): giab008. doi:10.1093/gigascience/giab008; Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. (2009) The Sequence Alignment/Map format and SAMtools. Bioinformatics 25(16): 2078–2079. doi:10.1093/bioinformatics/btp352
mosdepth: Pedersen BS, Quinlan AR (2018) Mosdepth: quick coverage calculation for genomes and exomes. Bioinformatics 34(5): 867–868. doi:10.1093/bioinformatics/btx699

RNA-seq: counts and abundances

After mapping, RNA-seq reads are counted. featureCounts counts the reads that fall in each gene; StringTie assembles the transcripts that the alignments show and estimates their abundance. kallisto and Salmon need no alignments to the genome: they estimate the abundance of each transcript from the reads and the transcript sequences.

The examples use the 256 pairs of yeast RNA-seq reads and chromosome I of yeast, and work in a folder of their own:

The alignments

featureCounts and StringTie work on a BAM file sorted by position. --dta makes HISAT2 report alignments for transcript assembly (it wants longer anchors across introns, as StringTie does).

hisat2-build /data/examples/rnaseq/yeast_chrI.fa yeast_chrI > /dev/null 2>&1
hisat2 --dta -x yeast_chrI -1 /data/examples/rnaseq/SRR6357070_1.fastq -2 /data/examples/rnaseq/SRR6357070_2.fastq 2> hisat2.log | samtools sort -o rnaseq.bam -
samtools index rnaseq.bam

featureCounts: reads per gene

featureCounts -p --countReadPairs -a /data/examples/rnaseq/yeast_chrI.gtf -o counts.txt rnaseq.bam 2> featurecounts.log
cat counts.txt.summary | awk '$2 > 0'
tail -n +3 counts.txt | sort -k7,7nr | head -5 | cut -f 1,6,7

counts.txt has a line per gene: its ID, where its exons are, its length and the count – a column for each BAM file, so featureCounts … -o counts.txt s1.bam s2.bam s3.bam makes the table of a whole experiment, the input of DESeq2 or edgeR. The summary says why reads were not counted: Unassigned_NoFeatures (no gene there), Unassigned_Ambiguity (more than one gene), Unassigned_MultiMapping (the read maps to several places; -M counts it anyway). Its report goes to the standard error – here featurecounts.log.

-p --countReadPairs: the reads are pairs, and each pair (each fragment) counts once. featureCounts 2.1 stops with "Paired-end reads were detected in single-end read library" if -p is missing for paired reads.
-s 1 or -s 2: a stranded library, where the reads come from one strand of the RNA. Most RNA-seq kits today are reverse-stranded (-s 2); Salmon tells which (below).
-t exon -g gene_id (the defaults): count over the exons, and add them up by gene. -f counts each exon by itself.

featureCounts counts on several threads – -T 2 to -T 4 here; one if not asked for more. That needs a cross-origin isolated page, which the workbench makes its page when it starts – see programs on several threads.

StringTie: transcripts and their abundance

stringtie rnaseq.bam -G /data/examples/rnaseq/yeast_chrI.gtf -o stringtie.gtf -A gene_abundance.tsv
awk '$3 == "transcript"' stringtie.gtf | wc -l
sort -t $'\t' -k9,9gr gene_abundance.tsv | head -5 | cut -f 1,2,7-9

StringTie builds the transcripts that the reads show, guided by the annotation (-G). Each is called STRG.N.M; one that matches an annotated transcript names it (reference_id, ref_gene_name). For each it gives the coverage, FPKM and TPM (in the GTF file); -A writes a table per gene. -e leaves out new transcripts and only estimates the known ones; -B adds the tables that the R package Ballgown reads. Assemblies of several samples are joined with stringtie --merge -G genes.gtf -o merged.gtf s1.gtf s2.gtf.

kallisto and Salmon: abundances without a BAM file

Both need the sequences of the transcripts. gffread cuts them out of the genome, as the annotation describes them:

gffread -w transcripts.fa -g /data/examples/rnaseq/yeast_chrI.fa /data/examples/rnaseq/yeast_chrI.gtf
grep -c '>' transcripts.fa

kallisto finds, for each read pair, which transcripts it could come from (pseudoalignment: compatible k-mers, no base-by-base alignment), and shares the reads out among the transcripts with an EM algorithm:

kallisto index -i transcripts.idx transcripts.fa > /dev/null 2>&1
kallisto quant -i transcripts.idx -o kallisto_out /data/examples/rnaseq/SRR6357070_1.fastq /data/examples/rnaseq/SRR6357070_2.fastq 2> kallisto.log
sort -t $'\t' -k4,4gr kallisto_out/abundance.tsv | head -5

Salmon maps the reads to the transcripts (selective alignment), learns the lengths of the fragments and the strand of the reads, and shares the reads out among the transcripts:

salmon index -t transcripts.fa -i salmon_index > /dev/null 2>&1
salmon quant -i salmon_index -l A -1 /data/examples/rnaseq/SRR6357070_1.fastq -2 /data/examples/rnaseq/SRR6357070_2.fastq -o salmon_out 2> salmon.log
sort -t $'\t' -k5,5gr salmon_out/quant.sf | head -5
grep expected_format salmon_out/lib_format_counts.json

abundance.tsv (kallisto) and quant.sf (Salmon) give for each transcript its length, its effective length, the estimated number of reads (est_counts, NumReads) and TPM. -l A lets Salmon find out how the library was made: here ISR – inward pairs, read 1 from the reverse strand. That is a reverse-stranded library: featureCounts -s 2, stringtie --rf, kallisto quant --rf-stranded. Try the counts with -s 2, and with -s 1.

Do the four agree? Join the counts of the genes – yeast genes here have one transcript each, named like the gene:

tail -n +3 counts.txt | cut -f 1,7 | sort > fc.tsv
tail -n +2 kallisto_out/abundance.tsv | cut -f 1,4 | sort > kallisto.tsv
tail -n +2 salmon_out/quant.sf | cut -f 1,5 | sort > salmon.tsv
join fc.tsv kallisto.tsv | join - salmon.tsv | sort -k2,2nr | head -8

Good to know

Counts are what DESeq2 and edgeR need: featureCounts' table as it is; for kallisto and Salmon, R's tximport sums the transcripts of each gene. TPM is for comparing genes within a sample, not for testing.
Indexes belong to the program that made them: bwa index, hisat2-build, kallisto index and salmon index each make their own. An index that Salmon made on a 64-bit computer cannot be read here, nor one made here there.
With thousands of reads the numbers here are an exercise, not biology: real samples have millions.
Reads of single cells carry a cell barcode and a UMI as well: kallisto bus and bustools count them per cell – see Single-cell RNA-seq.

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

featureCounts: Liao Y, Smyth GK, Shi W (2014) featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30(7): 923–930. doi:10.1093/bioinformatics/btt656
StringTie: Shinder I, Pertea G, Hu R, Rudnick Z, Pertea M (2026) StringTie3 improves total RNA-seq assembly by resolving nascent and mature transcripts. Nature Methods 23(6): 1126–1137. doi:10.1038/s41592-026-03080-3; Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL (2015) StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nature Biotechnology 33(3): 290–295. doi:10.1038/nbt.3122; Kovaka S, Zimin AV, Pertea GM, Razaghi R, Salzberg SL, Pertea M (2019) Transcriptome assembly from long-read RNA-seq alignments with StringTie2. Genome Biology 20: 278. doi:10.1186/s13059-019-1910-1
kallisto: Bray NL, Pimentel H, Melsted P, Pachter L (2016) Near-optimal probabilistic RNA-seq quantification. Nature Biotechnology 34(5): 525–527. doi:10.1038/nbt.3519
Salmon: Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C (2017) Salmon provides fast and bias-aware quantification of transcript expression. Nature Methods 14(4): 417–419. doi:10.1038/nmeth.4197
gffread: Pertea G, Pertea M (2020) GFF Utilities: GffRead and GffCompare. F1000Research 9: 304. doi:10.12688/f1000research.23297.2

Single-cell RNA-seq: bustools

In single-cell RNA-seq each cell is in a droplet of its own, and every read carries two tags: the barcode of its droplet – its cell – and a UMI, the tag of the RNA molecule it was copied from. kallisto bus pseudoaligns the reads and writes what each one says as a BUS record. bustools sorts the records, corrects the barcodes and counts the molecules of each gene in each cell.

The examples work in a folder of your own:

The reads and the BUS file

/data/examples/singlecell/ has data of mouse cells, prepared with 10x Chromium v2 (SRA run SRR8599150): 2,500 of its read pairs, and a small BUS file that kallisto bus wrote for reads of the experiment, with the files of the reference it used.

head -2 /data/examples/singlecell/R1.fastq
head -2 /data/examples/singlecell/R2.fastq
bustools inspect /data/examples/singlecell/input.bus

Read 1 has 26 bases: the cell's barcode (16) and the UMI (10). Read 2, 98 bases, is the cDNA – the part that says which transcript the molecule came from. kallisto bus (kallisto bus -i INDEX -o DIR -x 10xv2 R1.fastq R2.fastq) needs an index of the mouse transcriptome, which is too big for the browser; input.bus is what it wrote. Each record holds a barcode, a UMI, an equivalence class – the set of transcripts the read fits, numbered as in matrix.ec – and how many reads said exactly that:

bustools text -p /data/examples/singlecell/input.bus | head -5
head -3 /data/examples/singlecell/matrix.ec
head -3 /data/examples/singlecell/transcripts.txt

A line of matrix.ec is an equivalence class and its transcripts, by their line in transcripts.txt (counted from 0).

Sort, correct, sort again

bustools sort -o sorted.bus /data/examples/singlecell/input.bus
bustools correct -w /data/examples/singlecell/whitelist.txt -o corrected.bus sorted.bus
bustools sort -o corrected.sorted.bus corrected.bus

Sorting puts the records of the same barcode and UMI together, and merges those that are the same (34 records become 31). correct compares each barcode with the list of the technology – 10x v2 has 737,280 barcodes, and whitelist.txt holds those near the barcodes of this file: a barcode in the list stays (28 records), one that is a single letter from exactly one barcode of the list takes that barcode (1), and the rest are left out (2) – sequencing errors, or droplets that were not made by the kit. Correction changes barcodes, so the records are sorted again.

Count: a matrix of cells and genes

mkdir -p counts
bustools count --genecounts -g /data/examples/singlecell/transcripts-to-genes.txt -e /data/examples/singlecell/matrix.ec \
    -t /data/examples/singlecell/transcripts.txt -o counts/cells_x_genes corrected.sorted.bus
ls counts
head -5 counts/cells_x_genes.mtx

count counts, for each barcode and each gene, the different UMIs – molecules, not reads: the copies of one molecule that PCR made count once. -g gives each transcript its gene, --genecounts counts genes rather than equivalence classes. The matrix is in Matrix Market format, which lists only what is not zero: after the header, its size – barcodes, genes, entries that are not zero – then one line per entry: the barcode's row, the gene's column, the count. cells_x_genes.barcodes.txt and cells_x_genes.genes.txt name the rows and the columns. Here each barcode has one molecule: 34 records are too few for more.

More cells, and which barcodes are cells

large/ holds 10,235 records of 358 barcodes of the same experiment, already corrected and sorted, with the files of their reference:

L=/data/examples/singlecell/large
bustools inspect -e $L/matrix.ec -w $L/whitelist.txt $L/input.bus | head -12
mkdir -p large_counts
bustools count --genecounts -g $L/transcripts-to-genes.txt -e $L/matrix.ec -t $L/transcripts.txt -o large_counts/cells_x_genes $L/input.bus
head -3 large_counts/cells_x_genes.mtx

358 barcodes, 3,733 genes. Most droplets hold no cell – only a little free RNA of the cells that broke – so most barcodes have very few molecules, and a few have many: those are the cells. Count the molecules of each barcode in Python, and draw them from the most to the fewest (a "knee plot"):

import os
import numpy as np
import scipy.io
import matplotlib.pyplot as plt

folder = os.path.expanduser("~/sc/large_counts")
m = scipy.io.mmread(f"{folder}/cells_x_genes.mtx").tocsr()
umis = np.asarray(m.sum(axis=1)).ravel()   # molecules per barcode
genes = np.asarray((m > 0).sum(axis=1)).ravel()   # genes seen per barcode
print(m.shape, "– molecules per barcode, the most first:", np.sort(umis)[::-1][:16].astype(int))
print("barcodes with 100 molecules or more:", (umis >= 100).sum(), "– median molecules per barcode:", np.median(umis))
plt.loglog(np.arange(1, len(umis) + 1), np.sort(umis)[::-1], marker=".", linestyle="none")
plt.xlabel("barcodes, by their molecules"); plt.ylabel("molecules (UMIs)")
plt.show()

Fourteen barcodes have 100 molecules or more – from 151 to 1,230 –, then the count falls to a few dozen and the median barcode has one: the cells, and then the empty droplets. Where the curve bends (the "knee") is where cells are told from empty droplets; tools such as DropletUtils (R) do this with statistics. These 358 barcodes are the first of the experiment's in alphabetical order, not a sample of it: they show how it is done, not what the experiment found.

More: bustools capture keeps (or, with -x, leaves out) the records of a list of transcripts or barcodes, bustools allowlist makes a barcode list from the data, bustools umicorrect corrects UMIs, bustools compress and decompress pack a BUS file; bustools COMMAND without arguments lists its options. In the browser bustools runs on one thread – sort -t N becomes -t 1, with a note – and sorts in 64 MB of memory unless told otherwise (1 GB at most: -m 4G becomes -m 64M).

The same steps are in the Galaxy tab: the tools bustools sort, correct, count, inspect and text (section Single-cell) – upload input.bus, whitelist.txt, matrix.ec, transcripts.txt and transcripts-to-genes.txt, then sort, correct, sort and count. They are the workbench’s own: Galaxy servers have no bustools tools and process single-cell reads with others (RNA STARsolo, Alevin). bustools is compiled to WebAssembly from its released source by the genomics-wasm project, and gives the same records and matrices as bustools built for Linux (see Programs and versions).

References

The papers to cite when you use this program in your work – with the version you used (Programs and versions):

bustools (and kallisto bus): Melsted P, Booeshaghi AS, Liu L, Gao F, Lu L, Min KH, et al. (2021) Modular, efficient and constant-memory single-cell RNA-seq preprocessing. Nature Biotechnology 39(7): 813–818. doi:10.1038/s41587-021-00870-2; Melsted P, Ntranos V, Pachter L (2019) The barcode, UMI, set format and BUStools. Bioinformatics 35(21): 4472–4473. doi:10.1093/bioinformatics/btz279; Sullivan DK, Min KH, Hjörleifsson KE, Luebbert L, Holley G, Moses L, et al. (2025) kallisto, bustools and kb-python for quantifying bulk, single-cell and single-nucleus RNA-seq. Nature Protocols 20(3): 587–607. doi:10.1038/s41596-024-01057-0

Variant calling

bcftools finds the positions where the reads differ from the reference, decides the genotypes, filters the calls and turns VCF files into tables. FreeBayes is a second caller, which looks at haplotypes rather than single positions. bgzip and tabix compress and index the VCF files.

The examples work in a folder of your own:

From a BAM file to variants

First a sorted BAM file (see Mapping reads):

minimap2 -ax sr -R '@RG\tID:NA12878\tSM:NA12878' /data/examples/reads/reference.fa \
    /data/examples/reads/NA12878_R1.fastq /data/examples/reads/NA12878_R2.fastq 2> /dev/null \
  | samtools sort -o mapped.bam -
samtools index mapped.bam

Then the two steps of bcftools: mpileup collects the evidence at each position, call decides the genotypes:

bcftools mpileup -f /data/examples/reads/reference.fa mapped.bam | bcftools call -mv -Oz -o calls.vcf.gz
bcftools index calls.vcf.gz
bcftools view -H calls.vcf.gz | cut -f 1-6 | head

-m is the multiallelic caller, -v keeps only variant sites, -Oz writes compressed VCF (-Ov plain VCF, -Ob BCF).

Filter and count

bcftools filter -i 'QUAL>=20 && INFO/DP>=5' -Oz -o filtered.vcf.gz calls.vcf.gz
bcftools index filtered.vcf.gz
bcftools view -H filtered.vcf.gz | wc -l
bcftools view -H -v snps filtered.vcf.gz | wc -l
bcftools view -H -v indels filtered.vcf.gz | wc -l
bcftools stats filtered.vcf.gz | grep -E "^SN|^TSTV"

The expression language: QUAL, INFO/DP, FORMAT/GT="1/1", TYPE="snp", with &&, ||. -e excludes instead of including.

A table of the calls

bcftools query -H -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t%INFO/DP[\t%GT]\n' filtered.vcf.gz > calls.tsv
head calls.tsv

A table like this reads straight into pandas or R (pd.read_csv("calls.tsv", sep="\t")). A known position, such as the CYP2C19*2 variant (rs4244285, G>A) at position 11,616 of human_CYP2C19 – NA12878 carries it on one chromosome (0/1):

bcftools view -H -r human_CYP2C19:11616 filtered.vcf.gz | cut -f 1-10

FreeBayes

FreeBayes reads the same sorted BAM file and the reference (which needs its .fai index – samtools faidx; FreeBayes makes one if it is missing) and writes a VCF file:

freebayes -f /data/examples/reads/reference.fa mapped.bam > freebayes.vcf
bcftools view -H freebayes.vcf | wc -l
bcftools filter -i 'QUAL>=20' -Oz -o freebayes.filtered.vcf.gz freebayes.vcf
bcftools index freebayes.filtered.vcf.gz
bcftools view -H -r human_CYP2C19:11616 freebayes.filtered.vcf.gz | cut -f 1-6

FreeBayes also writes the sites it considered and does not believe – QUAL near 0, genotype 0/0 – so filter on QUAL. Its INFO fields differ from bcftools': DP depth, RO and AO the reads with the reference and the other allele, TYPE snp, mnp, ins, del or complex. Useful options: -r human_CYP2C19:11000-12500 (one region; needs the BAM index), -p 1 (a haploid organism), -C 3 (at least three reads for an allele), --min-mapping-quality 20. Which calls do the two callers share?

bcftools isec -n=2 -w1 filtered.vcf.gz freebayes.filtered.vcf.gz | grep -vc "^#"
bcftools isec -C -w1 filtered.vcf.gz freebayes.filtered.vcf.gz | grep -v "^#" | cut -f 1-5

-n=2 keeps the sites in both files, -C the sites of the first file only. Positions where the callers disagree are worth a look in the genome browser.

bgzip and tabix

A VCF file that tools read by region must be compressed with bgzip (not plain gzip) and indexed:

bcftools view filtered.vcf.gz > plain.vcf
bgzip -c plain.vcf > again.vcf.gz
tabix -p vcf again.vcf.gz
tabix again.vcf.gz human_CYP2C19:11000-12000 | cut -f 1-5

Also useful: bcftools norm (normalise indels against the reference), bcftools isec (compare two call sets), bcftools annotate. bcftools COMMAND --help lists the options of a command. The plugins (bcftools +NAME) cannot run in a browser.

igv filtered.vcf.gz mapped.bam shows the calls next to the reads in the Genome browser.

Variants of a tumour that its matched normal sample does not have – somatic variants – need callers that compare the two: see Somatic variants. One of them, lofreq call, also calls the variants of a single sample at low frequencies, as in a virus population.

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

BCFtools: Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. (2021) Twelve years of SAMtools and BCFtools. GigaScience 10(2): giab008. doi:10.1093/gigascience/giab008
FreeBayes: Garrison E, Marth G (2012) Haplotype-based variant detection from short-read sequencing. arXiv:1207.3907 [q-bio.GN]. arxiv.org/abs/1207.3907
HTSlib (bgzip, tabix): Bonfield JK, Marshall J, Danecek P, Li H, Ohan V, Whitwham A, et al. (2021) HTSlib: C library for reading/writing high-throughput sequencing data. GigaScience 10(2): giab007. doi:10.1093/gigascience/giab007

Somatic variants: tumour and normal

A somatic variant is one that the tumour has and the rest of the body does not. The callers here compare the reads of a tumour with the reads of a matched normal sample of the same person: SomaticSniper, LoFreq (lofreq somatic) and Strelka2.

The examples work in a folder of your own:

The example data

/data/examples/somatic/ holds real reads of the breast-cancer cell line HCC1395 (tumor.bam) and of HCC1395BL, a cell line made from the blood of the same donor (normal.bam). They come from the SEQC2 consortium, which sequenced the two cell lines many times over and published the somatic variants it is confident of. To keep them small, the reads are those of two windows of chromosome 21 (GRCh38), and the windows are the reference: chr21_10572001_10582000 (10 kb) and chr21_21404501_21423500 (19 kb). A position plus 10,572,000 – or plus 21,404,500 – is the position on chromosome 21.

ls /data/examples/somatic/
samtools view -H /data/examples/somatic/tumor.bam | grep -E "^@(SQ|RG)"
grep -v "^##" /data/examples/somatic/truth.vcf | cut -f 1-7

truth.vcf has the four SEQC2 variants in these windows: two SNVs and two deletions of one base. The SNV at chr21_21404501_21423500:14045 is a rare one (about 7 % of the tumour cells): in this particular set of reads no read shows it, and no caller can find it. The BAM files are sorted and indexed, and both name the two windows – the callers check that tumour, normal and reference agree.

Look at the first SNV yourself, tumour and normal together:

igv --genome /data/examples/somatic/reference.fa /data/examples/somatic/tumor.bam /data/examples/somatic/normal.bam --locus chr21_10572001_10582000:4,659-4,699

About a quarter of the tumour's reads have a G at position 4,679; none of the normal's do. That is what a somatic SNV looks like. (A germline variant would be in both, at about half or all of the reads.)

SomaticSniper

bam-somaticsniper -F vcf -f /data/examples/somatic/reference.fa \
    /data/examples/somatic/tumor.bam /data/examples/somatic/normal.bam somaticsniper.vcf
grep -v "^##" somaticsniper.vcf | cut -f 1,2,4,5,9-11

The tumour BAM comes first, then the normal, then the output file. Each line has the two samples, NORMAL then TUMOR, in the fields that FORMAT names. Two of them matter most: SS, the somatic status (0 wild type, 1 germline, 2 somatic, 3 loss of heterozygosity – LOH, 5 unknown), and SSC, the somatic score (higher is surer). Only the first call, the true SNV, has a high SSC; the others are LOH calls and weak candidates. To keep the somatic calls with a score of 40 or more:

bam-somaticsniper -F vcf -Q 40 -L -G -f /data/examples/somatic/reference.fa \
    /data/examples/somatic/tumor.bam /data/examples/somatic/normal.bam somaticsniper.strict.vcf
grep -v "^#" somaticsniper.strict.vcf | cut -f 1,2,4,5

-Q is the lowest somatic score reported, -L leaves out LOH calls, -G gains of reference. SomaticSniper calls SNVs only. -n and -t name the normal and tumour samples in the VCF – -t is not a number of threads.

LoFreq: lofreq somatic

LoFreq calls variants at any frequency the reads support. lofreq somatic calls the normal and the tumour, and keeps the tumour's calls that the normal does not have. For indels LoFreq needs indel qualities in the reads (lofreq indelqual); then the BAM files must be indexed again:

lofreq indelqual --dindel -f /data/examples/somatic/reference.fa -o tumor.iq.bam /data/examples/somatic/tumor.bam
lofreq indelqual --dindel -f /data/examples/somatic/reference.fa -o normal.iq.bam /data/examples/somatic/normal.bam
samtools index tumor.iq.bam && samtools index normal.iq.bam
lofreq somatic -n normal.iq.bam -t tumor.iq.bam -f /data/examples/somatic/reference.fa --call-indels -o lofreq_
ls lofreq_*
zcat lofreq_somatic_final.snvs.vcf.gz | grep -v "^##"
zcat lofreq_somatic_final.indels.vcf.gz | grep -v "^##"

Everything starts with the prefix given with -o. lofreq_somatic_final.snvs.vcf.gz and …indels.vcf.gz are the results; the other files are the steps on the way (relaxed and stringent calls of each sample, and the raw somatic calls). INFO has the allele frequency in the tumour (AF), the depth (DP) and a strand-bias score (SB). On the example LoFreq finds the SNV and both deletions. It warns that no dbSNP file was given: with human data, -d dbsnp.vcf.gz removes known germline variants. LoFreq's authors call this workflow experimental. lofreq call -f REF.fa -o calls.vcf sample.bam calls the variants of one sample, at low frequencies too.

Strelka2

Strelka2 is used in two steps, as on any computer. The first configures a run in a folder of its own and checks the input; the second runs it:

configureStrelkaSomaticWorkflow.py --normalBam /data/examples/somatic/normal.bam --tumorBam /data/examples/somatic/tumor.bam \
    --referenceFasta /data/examples/somatic/reference.fa --runDir strelka_run
strelka_run/runWorkflow.py -m local -j 1
ls strelka_run/results/variants/
bcftools view -H strelka_run/results/variants/somatic.snvs.vcf.gz | cut -f 1-7
bcftools view -H -f PASS strelka_run/results/variants/somatic.indels.vcf.gz | cut -f 1-7

Strelka writes every candidate it scored, and marks in FILTER whether it passed: PASS, or the filters it failed (LowEVS, a low empirical variant score; LowDepth). bcftools view -f PASS keeps the passing calls: here the SNV and both deletions. The normal and tumour columns hold the counts of each base (AU, CU, GU, TU) and, for indels, of reads with and without the indel (TAR, TIR); QSS and QSI are the somatic qualities. strelka_run/results/stats/runStats.tsv summarises the run.

Options of the first step: --exome (exome or panel data: no filter for unusually high depth), --callRegions regions.bed.gz (call only in these regions – compressed with bgzip and indexed with tabix -p bed; /data/examples/somatic/regions.bed.gz is one), --region chr:start-end, --outputCallableRegions; --allHelp lists them all. Each run needs a new folder. Strelka's workflow scripts are Python 2, which is not here: the terminal runs them itself, with Strelka's own programs. -j does not matter (the tasks run one after another), and there is no sge mode.

Compare the callers with the truth

grep -v "^#" /data/examples/somatic/truth.vcf | cut -f 1,2,4,5 > truth.tsv
grep -v "^#" somaticsniper.strict.vcf | cut -f 1,2,4,5 > sniper.tsv
zcat lofreq_somatic_final.snvs.vcf.gz lofreq_somatic_final.indels.vcf.gz | grep -v "^#" | cut -f 1,2,4,5 > lofreq.tsv
bcftools view -H -f PASS strelka_run/results/variants/somatic.snvs.vcf.gz | cut -f 1,2,4,5 > strelka.tsv
bcftools view -H -f PASS strelka_run/results/variants/somatic.indels.vcf.gz | cut -f 1,2,4,5 >> strelka.tsv
for f in sniper lofreq strelka; do echo "$f: $(grep -c -F -x -f truth.tsv $f.tsv) of $(wc -l < truth.tsv) true variants, $(wc -l < $f.tsv) calls"; done

Two windows chosen around known variants are no benchmark – they show how the callers report, not how good they are. Real data needs more: a panel of normals, known germline variants, a look at each candidate in the genome browser.

In Galaxy

Galaxy (the Galaxy tab) has Strelka Somatic, Insert indel qualities and Call variants with LoFreq, which run the same programs. Upload the two BAM files and the reference to a history first.

LoFreq, SomaticSniper and Strelka2 are compiled to WebAssembly from their released sources by the genomics-wasm project, and give the same calls as the same programs built for Linux (see Programs and versions).

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

LoFreq: The second paper for the indel qualities of lofreq indelqual --dindel, which are Dindel's.; Wilm A, Aw PPK, Bertrand D, Yeo GHT, Ong SH, Wong CH, et al. (2012) LoFreq: a sequence-quality aware, ultra-sensitive variant caller for uncovering cell-population heterogeneity from high-throughput sequencing datasets. Nucleic Acids Research 40(22): 11189–11201. doi:10.1093/nar/gks918; Albers CA, Lunter G, MacArthur DG, McVean G, Ouwehand WH, Durbin R (2011) Dindel: accurate indel calls from short-read data. Genome Research 21(6): 961–973. doi:10.1101/gr.112326.110
SomaticSniper: Larson DE, Harris CC, Chen K, Koboldt DC, Abbott TE, Dooling DJ, et al. (2012) SomaticSniper: identification of somatic point mutations in whole genome sequencing data. Bioinformatics 28(3): 311–317. doi:10.1093/bioinformatics/btr665
Strelka2: Kim S, Scheffler K, Halpern AL, Bekritsky MA, Noh E, Källberg M, et al. (2018) Strelka2: fast and accurate calling of germline and somatic variants. Nature Methods 15(8): 591–594. doi:10.1038/s41592-018-0051-x

Metagenomics: who is in the sample

The DNA of a sample of soil, gut or water comes from many organisms at once. Kraken2 gives each read the taxon – species, genus, family … – whose DNA it matches in a database of known genomes. Bracken then estimates how many reads each species has, including the reads that Kraken2 could place only on a genus or a family.

The examples work in a folder of your own:

The example data

/data/examples/metagenomics/ holds 4,096 pairs of real Illumina reads of a mock community: a mixture made to a known recipe, the ZymoBIOMICS Microbial Community DNA Standard – the DNA of eight bacteria, 12 % each, and of two yeasts, 2 % each. Because the answer is known, a mock community shows what a method gets right and what it gets wrong. db/ is a small Kraken2 database: the first million bases of three bacterial chromosomes – Escherichia coli, Salmonella enterica and Shigella flexneri – with their taxonomy.

ls /data/examples/metagenomics/ /data/examples/metagenomics/db/
zcat /data/examples/metagenomics/reads_1.fastq.gz | head -4
cat /data/examples/metagenomics/db/seqid2taxid.map

E. coli and Salmonella are in the community; Shigella is not – it is there as a close relative of E. coli. The other six bacteria and the yeasts are not in the database at all. seqid2taxid.map gives each sequence of the database its NCBI taxonomy ID: 562 is E. coli, 28901 S. enterica, 623 S. flexneri.

Kraken2: a taxon for every read

kraken2 --db /data/examples/metagenomics/db --paired --report sample.report --output sample.kraken \
    /data/examples/metagenomics/reads_1.fastq.gz /data/examples/metagenomics/reads_2.fastq.gz
head -5 sample.kraken

Kraken2 cuts each read into k-mers (31-base minimizers of 35-base windows, here) and looks them up in the database, which knows for each one the taxon it belongs to – the species, or, for a k-mer that several species share, their common ancestor. It puts the read on the taxon whose path from the root collects most of its k-mers. With --paired the two reads of a pair are classified together, so the 4,096 "sequences" are read pairs: 231 are classified, 3,865 are not.

A line of sample.kraken: C or U (classified or not), the name of the read, the taxonomy ID it was given (0: none), the lengths of the two reads, and its k-mers – 562:113 are 113 k-mers of taxon 562, 0:4 four k-mers that are not in the database, and |:| separates the first read from the second. A read of 151 bases has 117 k-mers of 35.

The report

cat sample.report

One line per taxon: the percentage of all reads in its clade (the taxon and everything below it), the reads in the clade, the reads given to this taxon itself, its rank (U unclassified, R root, D domain, K kingdom, P phylum, C class, O order, F family, G genus, S species; R1, S1 … are ranks in between), its taxonomy ID and its name, indented by depth. This report is what Bracken reads.

Most reads are unclassified. Only E. coli and Salmonella of the community are in the database, and only a fifth of their genomes (one megabase each). They are 24 % of the DNA, so about 5 % of the reads can be classified – and Kraken2 classifies 5.6 %. Kraken2 can find only what its database holds: a read of an organism that is not there stays unclassified, or goes to the nearest relative that is.
31 reads stop at the family, Enterobacteriaceae: their k-mers are shared by species of different genera, and Kraken2 does not guess between them.
9 reads go to Shigella flexneri, which is not in the community. Shigella is E. coli in all but name – the Shigella species are lineages within E. coli – and these are E. coli reads with k-mers that, of these three genomes, only Shigella's has. A database with hundreds of E. coli genomes would put most of them on E. coli. A name in a Kraken2 report is a hypothesis, not proof.

What the database knows decides all of this. kraken2-inspect shows how its minimizers are shared out:

kraken2-inspect --db /data/examples/metagenomics/db

A quarter of them belong to the family only – mostly what E. coli and Shigella have in common. Salmonella has 42 % to itself, E. coli and Shigella 18 % and 16 %: a read from the part of E. coli that Shigella shares can say no more than "Enterobacteriaceae".

Bracken: how much of each species

Bracken takes the reads that Kraken2 left on a genus or a family and shares them out among the species below, in proportion to how often each species gives reads that end there. It knows that from the database: database151mers.kmer_distrib, made by bracken-build for reads of 151 bases – the reads' length, which you give with -r.

bracken -d /data/examples/metagenomics/db -i sample.report -o sample.bracken -r 151 -l S -t 10
cat sample.bracken

-l S asks for species (G for genera, F for families …), and -t 10 leaves out species with fewer than 10 reads of their own. The table has, per species, Kraken2's reads, the reads added, the new estimate and its share of all the estimates. E. coli gets 30 of the 31 family reads – Salmonella's reads seldom end at the family, E. coli's often do, because of Shigella. Shigella, with 9 reads, is below the threshold and is left out. The two species of the community that the database holds come out at 45 % and 55 %; in the standard they are equal. With -t 1 every species counts:

bracken -d /data/examples/metagenomics/db -i sample.report -o sample_t1.bracken -r 151 -l S -t 1
cat sample_t1.bracken
cat sample_bracken_species.report

Now Shigella has 11 reads – a species that is not in the sample. The threshold is a guard against such names; it cannot tell a true rare species from a false one. Bracken also writes a report in Kraken2's format with the new counts (-w names it; by default it is the report's name with _bracken_species), which shows only the classified reads: the percentages there are of those, not of all reads.

Your own database

A Kraken2 database is built from three things: the sequences, in library/ (FASTA files ending in .fna); seqid2taxid.map, the taxonomy ID of each sequence; and taxonomy/ with NCBI's nodes.dmp (each taxon's parent and rank) and names.dmp (the names). The example has all three. Copy it, remove what was built from it, and build it again:

cp -r /data/examples/metagenomics/db ~/mydb
rm ~/mydb/*.k2d ~/mydb/*.kmer_distrib
ls -R ~/mydb
kraken2-build --build --db ~/mydb
kraken2-inspect --db ~/mydb --skip-counts
cmp ~/mydb/hash.k2d /data/examples/metagenomics/db/hash.k2d && echo "the same hash table"

kraken2-build --build estimates the size of the hash table, then fills it: hash.k2d (the minimizers and their taxa), taxo.k2d (the taxonomy) and opts.k2d (the settings: --kmer-len, --minimizer-len …). Then Bracken's distribution for reads of 151 bases:

bracken-build -d ~/mydb -t 1 -k 35 -l 151
cmp ~/mydb/database151mers.kmer_distrib /data/examples/metagenomics/db/database151mers.kmer_distrib && echo "the same distribution"

bracken-build classifies the database's own sequences with Kraken2 (database.kraken, kept for the next read length), cuts them into every possible read of the length given with -l, and counts where the reads of each genome end. -k is the database's k-mer length (35, Kraken2's default). Make one distribution for each read length you have.

On a computer, databases for real samples are made with kraken2-build --standard or --download-library, which download genomes and NCBI's taxonomy – tens of gigabytes for the standard database – or downloaded ready-made; Galaxy servers have them installed. Neither is possible in the browser, which holds a few hundred megabytes at most. Here you can build small databases of your own sequences: put them in library/, write seqid2taxid.map, and take the lineages from NCBI's taxonomy into nodes.dmp and names.dmp.

Options worth knowing

--confidence 0.1: a read is placed only where at least that share of its k-mers agrees – it moves up the tree, or stays unclassified (default 0). --minimum-hit-groups 3 asks for more evidence before a call (default 2). --use-names writes names in --output instead of numbers. --classified-out 'cl#.fq' --unclassified-out 'un#.fq' writes the reads themselves (the # becomes _1 and _2 for pairs) – the unclassified ones are what to look at next. --report-zero-counts lists every taxon of the database, --use-mpa-style writes the report in MetaPhlAn's format (d__Bacteria|k__Pseudomonadati|…). --output - writes no classifications. man kraken2, man bracken, kraken2-build --help and bracken -h have the rest.

In the browser Kraken2 and Bracken run on one thread (--threads and -t of bracken-build are left out, with a note), reads compressed with bzip2 cannot be read (gzip can), and nothing is downloaded. Bracken's two Python programs start a Python of their own for each command, which takes a few seconds.

In Galaxy

Galaxy (the Galaxy tab) has Kraken2 and Bracken in its Metagenomics section, with the example database. Upload the two read files, run Kraken2 with Print a report set to Yes, and give the report to Bracken.

Kraken2 and Bracken are compiled to WebAssembly from their released sources by the genomics-wasm project, and give the same classifications, reports and estimates as the same programs built for Linux (see Programs and versions).

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

Kraken 2: Wood DE, Lu J, Langmead B (2019) Improved metagenomic analysis with Kraken 2. Genome Biology 20: 257. doi:10.1186/s13059-019-1891-0; Lu J, Rincon N, Wood DE, Breitwieser FP, Pockrandt C, Langmead B, et al. (2022) Metagenome analysis using the Kraken software suite. Nature Protocols 17(12): 2815–2839. doi:10.1038/s41596-022-00738-y
Bracken: Lu J, Breitwieser FP, Thielen P, Salzberg SL (2017) Bracken: estimating species abundance in metagenomics data. PeerJ Computer Science 3: e104. doi:10.7717/peerj-cs.104

Genes and protein domains: Prodigal and HMMER

A new bacterial genome is a long string of letters. Prodigal finds the genes in it that code for proteins. HMMER then compares proteins – or DNA – with profile HMMs: models of a whole family, built from an alignment of its members, which find relatives that a search with a single sequence can miss.

The examples work in a folder of your own:

Prodigal: the genes of a bacterium

/data/examples/bacteria/ has the first 60,000 bases of the chromosomes of Escherichia coli K-12 and Salmonella enterica Typhimurium, from NCBI RefSeq.

prodigal -i /data/examples/bacteria/ecoli.fna -o ecoli.gff -f gff -a ecoli.faa -d ecoli.ffn
head -5 ecoli.gff
grep -c '>' ecoli.faa

Prodigal first learns the genome – how rich in G and C its coding regions are, which start codons and ribosome binding sites it uses – and then finds its genes. It warns that it would like 100,000 bases or more to learn from; 60,000 are enough to see how it works. It finds 56 genes. A line of ecoli.gff is one gene: the sequence, the program, the type (CDS, coding sequence), start and end, a score, the strand, and in the last column what Prodigal knows about it – partial (10: it runs off the left end of the sequence), the start codon, the ribosome binding site (rbs_motif) and its confidence. The second gene, from 337 to 2,799 on the + strand, is thrA, just where the annotation of the genome has it. The first runs off the start of the window, so Prodigal cannot know where it begins.

ecoli.faa holds the proteins and ecoli.ffn the genes' DNA, one record per gene, named after the sequence and the gene's number:

head -2 ecoli.faa
prodigal -i /data/examples/bacteria/salmonella.fna -o salmonella.gff -f gff -a salmonella.faa -q
prodigal -p meta -i /data/examples/bacteria/combined.fna -o meta.gff -f gff -q
grep -c CDS ecoli.gff salmonella.gff meta.gff

-q keeps Prodigal quiet but for warnings. -p meta does not learn: it uses 50 models that come with Prodigal, for genomes of different GC content – right for a metagenome, where the contigs come from many organisms, and for sequences too short to learn from. combined.fna holds both windows. Other options: -c (no genes running off the ends), -g 4 (another genetic code: 11 is that of bacteria and archaea), -s FILE (every possible start, with its scores), -t FILE (save what it learnt, and use it again).

HMMER: a profile from an alignment

/data/examples/hmmer/ holds the tutorial files of HMMER 3.4, under their own names, so the tutorial of the HMMER User's Guide can be followed here as it is written. First a profile of the globins, built from an alignment of four of them, and a search of 45 globins with it:

hmmbuild globins4.hmm /data/examples/hmmer/globins4.sto
hmmsearch globins4.hmm /data/examples/hmmer/globins45.fa > globins.out
head -24 globins.out

The profile has 149 match states (M=149): for each column of the alignment, how likely each amino acid is there, and how likely a gap. hmmsearch ranks the sequences by E-value – the number of hits this good that a search of a database of this size would give by chance – and its score in bits. All 45 globins are found, with E-values of 1e-36 or smaller, although the profile was built from only four. The output then shows each sequence's domains and their alignments to the profile. For tables, one line per sequence and one per domain:

hmmsearch --tblout globins.tbl --domtblout globins.domtbl --noali globins4.hmm /data/examples/hmmer/globins45.fa > /dev/null
grep -v '^#' globins.tbl | head -3

Which domains does a protein have? hmmscan

The other way round: one protein against a collection of profiles, as Pfam's are used. The tutorial joins three profiles – globins, fibronectin type III (fn3) and protein kinase – into a small collection, which hmmpress prepares for hmmscan, and scans the Sevenless receptor of Drosophila (7LESS_DROME, a UniProt record, which HMMER reads as it is):

cat /data/examples/hmmer/globins4.hmm /data/examples/hmmer/fn3.hmm /data/examples/hmmer/Pkinase.hmm > minifam.hmm
hmmpress minifam.hmm
hmmscan --domtblout 7less.domtbl minifam.hmm /data/examples/hmmer/7LESS_DROME > 7less.out
grep -v '^#' 7less.domtbl | awk '{print $1, $20, $21, $13}'

Nine fibronectin type III domains and one kinase domain, with where each lies on the protein (columns 20 and 21: from, to) and its E-value (i-Evalue, column 13); a globin, none. Two of the fn3 domains are weak, with E-values of 0.5 and 0.19: a domain like that is a candidate to look at, not a finding. The InterPro card of the Databases tab shows what InterPro, with the whole of Pfam, finds in the same protein (UniProt P13368).

One sequence as the query: phmmer and jackhmmer

phmmer --tblout hbb.tbl /data/examples/hmmer/HBB_HUMAN /data/examples/hmmer/globins45.fa > hbb.out
grep -v '^#' hbb.tbl | head -3
jackhmmer -N 3 /data/examples/hmmer/HBB_HUMAN /data/examples/hmmer/globins45.fa | grep '^@@'

phmmer searches with one protein, as BLASTP does, with HMMER's scores and E-values. jackhmmer searches again and again: it builds a profile from what it found and searches with that, until nothing new comes in – here after two rounds, as all 45 globins are in the first.

DNA: nhmmer

nhmmer --tblout made1.tbl /data/examples/hmmer/MADE1.hmm /data/examples/hmmer/dna_target.fa > made1.out
grep -v '^#' made1.tbl | cut -c1-140

MADE1 is a short DNA transposon of the human genome, from Dfam. nhmmer searches both strands: two copies in this piece of chromosome 1, near 174,500 and 302,400, are each found on both strands, as MADE1 reads almost the same backwards on the other strand (its ends are inverted repeats); the fifth hit, with an E-value above 1, is noise.

From genes to proteins to relatives

The proteins that Prodigal predicted can be searched like any others. E. coli's ThrA, the second gene, against Salmonella's proteins:

awk '/^>/{p=($1==">NC_000913.3_1_60000_2")} p' ecoli.faa > thrA.faa
phmmer --tblout thrA.tbl thrA.faa salmonella.faa > /dev/null
grep -v '^#' thrA.tbl | cut -c1-100
hmmscan --tblout ecoli_domains.tbl minifam.hmm ecoli.faa > /dev/null
grep -v '^#' ecoli_domains.tbl | wc -l

One hit, the first gene of the Salmonella window – its ThrA – with an E-value of 0, too small for the computer to write. None of E. coli's 56 proteins has a globin, fibronectin or kinase domain of these three profiles: the table has no lines but its comments. To annotate a genome, proteins are scanned against all of Pfam's families – more than 20,000 – which is too big for the browser: use the InterPro card, or a computer.

More: hmmalign globins4.hmm globins45.fa aligns sequences to a profile; hmmfetch takes one profile out of a collection, hmmstat lists them, hmmemit writes sequences a profile would produce, hmmconvert converts the format; NAME -h lists a program's options. Here HMMER runs on one thread – --cpu is left out, with a note – and cannot read files compressed with gzip (.gz): unpack them first.

In Galaxy

Galaxy (the Galaxy tab) has Prodigal Gene Predictor and hmmbuild, hmmsearch, hmmscan, phmmer, nhmmer and hmmalign in its Annotation section, which run the same programs.

Prodigal and HMMER are compiled to WebAssembly from their released sources by the genomics-wasm project, and give the same genes, scores and alignments as the same programs built for Linux (see Programs and versions).

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

Prodigal: Hyatt D, Chen GL, LoCascio PF, Land ML, Larimer FW, Hauser LJ (2010) Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics 11: 119. doi:10.1186/1471-2105-11-119
HMMER: The HMMER User’s Guide asks you to cite the web site, hmmer.org, and the version you used; where a journal wants a paper, Eddy (2011).; HMMER: biosequence analysis using profile hidden Markov models. hmmer.org; Eddy SR (2011) Accelerated Profile HMM Searches. PLoS Computational Biology 7(10): e1002195. doi:10.1371/journal.pcbi.1002195

Genotypes and association: PLINK 2

PLINK works with the genotypes of many people at once – thousands of people, millions of variants. It checks their quality, counts allele frequencies, measures linkage disequilibrium and relatedness, and tests each variant for association with a trait: a genome-wide association study, or GWAS. Every PLINK command has the same shape: the data (--vcf, --pfile, --bfile), filters, what to compute, and --out, the name of the files it writes.

The examples work in a folder of your own:

The data, and PLINK's own files

/data/examples/plink2/1kg-phase3-chr21.vcf.gz holds the genotypes of the 2,504 people of the 1000 Genomes Project (Phase 3) at the first 600 variants of chromosome 21 – real public data, and the test file of PLINK 2 itself. A VCF file has a column for each person:

zcat /data/examples/plink2/1kg-phase3-chr21.vcf.gz | grep -v '^##' | head -3 | cut -f 1-5,10-12
plink2 --vcf /data/examples/plink2/1kg-phase3-chr21.vcf.gz --make-pgen --out chr21
ls chr21.*

A genotype such as 0|1 gives the alleles of the two copies of the chromosome – 0 the reference allele (REF), 1 the alternative (ALT); | means phased: which copy carries which is known. --make-pgen writes PLINK 2's own files: chr21.pgen, the genotypes, packed tightly; chr21.pvar, the variants (the first columns of the VCF); chr21.psam, the people; and chr21.log, a record of what was done. --pfile chr21 reads the three back – much quicker than the VCF for every new command. PLINK 1's files, .bed, .bim and .fam, are still common: --make-bed writes them and --bfile reads them.

The log, and what PLINK prints, begin with the options in effect: besides yours, --threads 1 and --memory 640 require, which the terminal adds (see the end of this page).

Quality control: frequencies, missing genotypes, Hardy–Weinberg

plink2 --pfile chr21 --freq --missing --hardy --out qc
head -4 qc.afreq
awk 'NR > 1 && ($5 < 0.01 || $5 > 0.99)' qc.afreq | wc -l
tail -n +2 qc.hardy | sort -g -k10,10 | head -4

qc.afreq has the frequency of each ALT allele and the number of alleles counted – 5,008, two for each person. Most variants are rare: for 566 of the 600 the less common allele has a frequency below 1%, and hundreds were seen only once among the 5,008 chromosomes (a frequency of 0.0002). qc.smiss and qc.vmiss count the missing genotypes of each person and of each variant: none here, as the 1000 Genomes release gives every genotype – in data of your own, a person or a variant with many missing genotypes is a warning.

qc.hardy holds the Hardy–Weinberg exact test: are there as many heterozygotes as the allele frequency predicts? O(HET_A1) is the share of heterozygotes observed, E(HET_A1) the share expected, P the p-value. A few variants are very far from it, always with too many heterozygotes – 0.88 where 0.50 is expected. These variants lie on the short arm of chromosome 21, near its centromere, among repeated sequences – most likely, reads from two copies of a repeat were aligned to one place, and the difference between the copies looks like a heterozygous genotype in almost everyone. Such errors are what the Hardy–Weinberg filter is for.

Filters

plink2 --pfile chr21 --maf 0.01 --hwe 1e-6 --make-pgen --out chr21_qc

Filters act before anything else. --maf 0.01 keeps the variants whose minor allele frequency is 1% or more; --hwe 1e-6 removes those with a Hardy–Weinberg p-value below 10⁻⁶; --geno 0.02 would remove variants with more than 2% of their genotypes missing, and --mind 0.02 such people. The log says what each did: 25 variants removed by --hwe (all with too many heterozygotes), 566 by --maf, 9 remaining. --keep FILE and --remove FILE choose people by their IDs, --extract FILE and --exclude FILE variants, --chr with --from-bp and --to-bp a region.

Linkage disequilibrium

plink2 --pfile chr21_qc --indep-pairwise 50 5 0.2 --out ld
cat ld.prune.in
plink2 --pfile chr21_qc --r2-unphased --ld-window-r2 0.2 --out ld
cat ld.vcor

Variants close to each other are inherited together, so their genotypes are correlated: linkage disequilibrium, measured as r². --indep-pairwise 50 5 0.2 moves along the chromosome in windows of 50 variants, 5 at a time, and of every pair with r² above 0.2 keeps one: ld.prune.in lists the variants kept (6 of the 9), ld.prune.out the others. Analyses that expect independent variants – principal components, relatedness – use the pruned set (--extract ld.prune.in). --r2-unphased lists the pairs: rs559962037 and rs556687069, 229 bases apart, have r² = 0.96 – they are nearly always inherited together.

Association: does a variant go with a trait?

The 1000 Genomes data has no trait to study, so these examples use synthetic data made for PLINK's own tests: 603 invented people, 30 variants (v1 … v30) and three phenotypes – BINARY (case 2 or control 1) and QUANT (a number), both drawn at random, and FIRTH, which makes every carrier of v1's ALT allele a case – with two covariates, C1 and C2 (in a real study: age, sex, principal components).

D=/data/examples/plink2
head -3 $D/synthetic_phenotypes.tsv
plink2 --vcf $D/synthetic_glm.vcf --pheno $D/synthetic_phenotypes.tsv --covar $D/synthetic_covariates.tsv --glm hide-covar --out assoc
ls assoc.*
cut -f 3,9,11-16 assoc.QUANT.glm.linear | head -5

--glm tests every variant with every phenotype: linear regression for a quantitative trait (assoc.QUANT.glm.linear), logistic regression for cases and controls (.glm.logistic.hybrid), with the covariates in each model (hide-covar: their own lines are left out). A line gives the variant, the allele tested (A1) and its frequency, the number of people, the effect – BETA, the change of the trait for each copy of A1, or OR, the odds ratio – its standard error, the test statistic and the p-value. QUANT and BINARY are random, so their p-values are spread between 0 and 1: of 30 tests, one or two fall below 0.05 by chance (here one for QUANT, none for BINARY). That is why a study is judged against the number of its tests – 0.05/30 here, and 5×10⁻⁸ for the million tests of a genome-wide study.

sort -g -k16,16 assoc.FIRTH.glm.logistic.hybrid | cut -f 3,10,13,16,17 | head -4

For FIRTH, v1 stands out. Every carrier of its ALT allele is a case – complete separation, for which ordinary logistic regression finds no odds ratio: it grows without end. PLINK then uses Firth's regression (FIRTH? Y), and still marks the result UNFINISHED: the odds ratio is in effect infinite. ERRCODE says when a result is not to be trusted as it stands. No other variant passes 0.05/30 = 0.0017. The p-values as a picture, in Python (a small "Manhattan plot"):

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

os.chdir(os.path.expanduser("~/gwas"))
fig, ax = plt.subplots(figsize=(7, 3))
for name, file in [("QUANT", "assoc.QUANT.glm.linear"), ("BINARY", "assoc.BINARY.glm.logistic.hybrid"),
                   ("FIRTH", "assoc.FIRTH.glm.logistic.hybrid")]:
    t = pd.read_csv(file, sep="\t")
    ax.plot(t["POS"] / 1000, -np.log10(t["P"]), "o", markersize=4, label=name)
ax.axhline(-np.log10(0.05 / 30), linestyle="--", color="grey", linewidth=1)   # 0.05 for 30 tests
ax.set_xlabel("position (kb)")
ax.set_ylabel("−log10 p")
ax.legend()
plt.show()

Scores

head -3 /data/examples/plink2/synthetic_scores.tsv
plink2 --pfile chr21 --score /data/examples/plink2/synthetic_scores.tsv 1 2 3 header-read cols=+scoresums --out score
head -4 score.sscore

A score adds up, for each person, the copies of chosen alleles, each multiplied by a weight. --score FILE 1 2 3 takes the variant, the allele and the weight from columns 1, 2 and 3 of the file (header-read: its first line is a header). score.sscore has, for each person, the alleles counted (ALLELE_CT, two for each of the 20 variants), how many of them were the named allele, and the score – the average over the alleles (WEIGHT_AVG) and, with cols=+scoresums, the sum (WEIGHT_SUM). A polygenic risk score is made this way, with the effect sizes of a published GWAS as the weights. The weights of synthetic_scores.tsv are made up: these scores mean nothing.

More, and what is not here

--export vcf (vcf bgz: compressed) writes a VCF file again, --make-bed PLINK 1's files, and --export A a table of 0, 1 or 2 copies for each person and variant, to read in Python or R.
--make-king-table estimates the kinship of each pair of people (KING), and --king-cutoff removes one of each related pair. It needs thousands of common variants: on 600 variants, most of them rare, its numbers are noise. /data/examples/plink2/keep48.samples.txt lists 48 of the people, for --keep.
plink2 --help --glm explains an option; plink2 alone lists the main ones. The manual is at cog-genomics.org/plink/2.0.
Not in this build: --pca (principal components, which correct an association study for ancestry) and --pheno-svd need the LAPACK library, which this WebAssembly build of PLINK 2 does not have – PLINK says so and stops. --script and --rerun (options read from a file) are not available either; plink2 $(cat FILE) reads options from a file.

PLINK 2 here is the development version of 12 October 2026 (2.0.0-b.1), compiled to WebAssembly from its source by the genomics-wasm project; it gives the same files as PLINK 2 built for Linux from the same source (see Programs and versions). It reads its files on threads of their own, which needs a cross-origin isolated page – the workbench makes its page so as it starts (see Limits). It computes on one thread and in 640 MiB of working memory, the least that PLINK 2 accepts: the terminal adds --threads 1 --memory 640 require to each command unless you give them – --threads 2 and --memory 1024 at most. In the Galaxy tab, the tool PLINK 2 runs the same program on a VCF of the history.

The 1000 Genomes data: The 1000 Genomes Project Consortium (2015) A global reference for human genetic variation. Nature 526: 68–74. doi:10.1038/nature15393

References

The papers to cite when you use this program in your work – with the version you used (Programs and versions):

PLINK 2: PLINK's citation instructions (cog-genomics.org/plink/2.0/general_usage#cite) ask for the program – its version and web site – and the paper of PLINK 2, and also the papers of the methods you used (the Hardy–Weinberg exact test, KING, Firth regression …), which they list.; Chang CC, Chow CC, Tellier LCAM, Vattikuti S, Purcell SM, Lee JJ (2015) Second-generation PLINK: rising to the challenge of larger and richer datasets. GigaScience 4: 7. doi:10.1186/s13742-015-0047-8; Purcell S, Chang C, Demaille B. PLINK 2.0. www.cog-genomics.org/plink/2.0

Peaks: ChIP-seq and ATAC-seq (MACS3)

In ChIP-seq, the DNA that a protein holds – a transcription factor, a modified histone – is pulled out of the cell with an antibody against the protein and sequenced, so its reads pile up where the protein sits on the genome. MACS3 finds those piles, the peaks: places with more reads than chance – and the control experiment – would put there. Its command hmmratac does the same for ATAC-seq, whose reads mark the open chromatin.

The examples work in a folder of your own:

The reads

/data/examples/chipseq/ holds MACS's own test data: the reads of a ChIP-seq of CTCF – a protein that binds one DNA motif at tens of thousands of places in the genome and folds the chromatin into loops – and of its control, the same chromatin sequenced without the antibody (the "input"). There are about 50,000 of each, all on human chromosome 22, at their positions on hg19 (GRCh37).

zcat /data/examples/chipseq/CTCF_SE_ChIP_chr22_50k.bed.gz | head -3
zcat /data/examples/chipseq/CTCF_SE_ChIP_chr22_50k.bed.gz | wc -l
zcat /data/examples/chipseq/CTCF_SE_CTRL_chr22_50k.bed.gz | wc -l

Each line is a read where it mapped: chromosome, start and end (BED counts from 0) and its strand. The reads are 101 bases long, but the DNA fragments they were read from are longer: a read is one end of its fragment, on the strand it was read from. MACS3 reads BAM files as well (-f BAM) – sorted or not, it needs no index – and the pairs of a paired-end run (-f BAMPE, -f BEDPE; see below).

Calling the peaks

macs3 callpeak -t /data/examples/chipseq/CTCF_SE_ChIP_chr22_50k.bed.gz -c /data/examples/chipseq/CTCF_SE_CTRL_chr22_50k.bed.gz \
    -f BED -g 52000000 -n ctcf -B --outdir ctcf

-t is the ChIP (the "treatment"), -c the control, -f the format of both. -g is the effective genome size, the part of the genome that reads can map to: hs for human, mm for mouse, ce and dm for worm and fly, or a number – here the reads come from chromosome 22 alone, about 51 million bases, and MACS's own tests give 52,000,000. -n names the files, --outdir is their folder, and -B writes the pile-ups as bedGraph files as well. MACS3 is a Python program: each command starts the Python that comes with it, with NumPy, SciPy and the rest, which takes a few seconds – the first time a little longer, as about 45 MB are loaded.

Its messages say what it does, step by step ([RSS unavailable]: in the browser it cannot measure the memory it uses).

#1 it reads the tags – the reads – and keeps at most one at the same position on the same strand (--keep-dup 1): more are taken for copies made by PCR. 49,622 ChIP reads become 48,047.
#2 it builds its model: around a binding site, the reads of the + strand pile up on its left and those of the − strand on its right, half a fragment away on each side. From 469 pairs of such piles MACS works out the fragment length d: 228 bases. Each read is then extended to d bases, in its direction, and the fragments are piled up.
#3 it calls the peaks: at each position the ChIP's pile-up is compared with the number of fragments expected there, the local λ – the most of the control's pile-up over d, 1,000 and 10,000 bases around, and of the genome's average, scaled to the same number of reads – with a Poisson test. Its p-values become q-values (Benjamini–Hochberg), and the regions with q below 0.05 (-q 0.05) are the peaks. The control is what makes the test fair: regions that every sample has more reads in – open chromatin, repeats, copies of the genome – have more control reads too.

The peaks

ls ctcf
head -3 ctcf/ctcf_peaks.narrowPeak
wc -l ctcf/ctcf_peaks.narrowPeak
sort -k9,9gr ctcf/ctcf_peaks.narrowPeak | head -3

730 peaks. A line of ctcf_peaks.narrowPeak (ENCODE's narrowPeak format: BED with four more columns) is a peak: chromosome, start, end, name; a score (−10·log₁₀q, as a whole number of at most 1000); the strand (none: .); the fold enrichment at the summit – the pile-up over the local λ –, −log₁₀p and −log₁₀q; and where the summit, the highest point of the pile-up, lies, counted from the start. The strongest peak (the last command sorts by the q-value) has 58 times as many fragments at its summit as expected, and a q-value of 10⁻¹⁸⁸. The other files: ctcf_peaks.xls, the same peaks as a table for a spreadsheet, with the command and the settings at the top (a text file, despite its name; open it in the Files tab); ctcf_summits.bed, the summit of each peak; ctcf_treat_pileup.bdg and ctcf_control_lambda.bdg (with -B), the pile-up of the ChIP and the local λ along the chromosome; and ctcf_model.r, an R script that draws the model.

These are MACS's own results: its test of the same command keeps the peaks it expects in run_callpeak_narrow0_peaks.narrowPeak. Apart from the names, which come from -n, the two files are the same, byte for byte:

diff <(cut -f 1-3,5- ctcf/ctcf_peaks.narrowPeak) <(cut -f 1-3,5- /data/examples/chipseq/run_callpeak_narrow0_peaks.narrowPeak) && echo the same

On the genome

The fold enrichment along the chromosome is the ChIP's pile-up divided by the local λ: macs3 bdgcmp works it out from the two bedGraph files. Then the peaks and the enrichment go to the Genome browser, on hg19 with its genes (a genome hosted by the IGV team: it needs the internet).

macs3 bdgcmp -t ctcf/ctcf_treat_pileup.bdg -c ctcf/ctcf_control_lambda.bdg -m FE -o ctcf/ctcf_FE.bdg
igv --genome hg19 ctcf/ctcf_peaks.narrowPeak ctcf/ctcf_FE.bdg --locus chr22:39,700,000-39,800,000

Nine peaks in these 100,000 bases. -m chooses what bdgcmp computes: FE, logFE, logLR, ppois (−log₁₀ of the Poisson p-value at each position), qpois, subtract …; macs3 bdgpeakcall calls peaks from such a track with a cutoff of your own.

The model is drawn by R: the script writes ctcf_model.pdf in the folder it runs in.

cd ctcf && Rscript ctcf_model.r && cd ..
open ctcf/ctcf_model.pdf

Its first page shows the two piles of tags, the + strand (red) left and the − strand (blue) right of the binding sites; its second the correlation of the two strands for each shift between them, highest at 228.

The peaks in Python

import os
import pandas as pd
import matplotlib.pyplot as plt

cols = ["chrom", "start", "end", "name", "score", "strand", "fold_enrichment", "log10_p", "log10_q", "summit"]
peaks = pd.read_csv(os.path.expanduser("~/peaks/ctcf/ctcf_peaks.narrowPeak"), sep="\t", names=cols)
peaks["length"] = peaks["end"] - peaks["start"]
print(len(peaks), "peaks; their length:", peaks["length"].min(), "to", peaks["length"].max(), "– median", peaks["length"].median())
print(peaks.nlargest(3, "log10_q")[["chrom", "start", "end", "fold_enrichment", "log10_q"]])
plt.scatter(peaks["fold_enrichment"], peaks["log10_q"], s=8)
plt.xscale("log")
plt.xlabel("fold enrichment at the summit"); plt.ylabel("−log10 q")
plt.show()

No peak is shorter than d, 228 bases (--min-length, which is d unless given). The stronger the enrichment, the smaller the q-value – but not in step: the q-value depends on the numbers of reads too, and a few fragments over a very low λ make a large fold enrichment with a modest q.

Paired-end reads

With pairs MACS3 needs no model: each pair shows the whole fragment, from one read to its mate. The same experiment, read from both ends:

zcat /data/examples/chipseq/CTCF_PE_ChIP_chr22_50k.bedpe.gz | head -2
macs3 callpeak -t /data/examples/chipseq/CTCF_PE_ChIP_chr22_50k.bedpe.gz -c /data/examples/chipseq/CTCF_PE_CTRL_chr22_50k.bedpe.gz \
    -f BEDPE -g 52000000 -n ctcf_pe --outdir ctcf_pe
wc -l ctcf_pe/ctcf_pe_peaks.narrowPeak
bedtools intersect -u -a ctcf/ctcf_peaks.narrowPeak -b ctcf_pe/ctcf_pe_peaks.narrowPeak | wc -l

MACS's BEDPE is a fragment per line: chromosome, start and end (not the BEDPE of bedtools). The fragments are 253.5 bases long on average – the model's 228 was close. 715 peaks, and 693 of the 730 above overlap one of them. From a BAM file of paired reads, -f BAMPE does the same.

ATAC-seq: open chromatin, with hmmratac

In ATAC-seq the enzyme Tn5 cuts the DNA where it is not wrapped around nucleosomes and adds the adapters of the sequencer there; the fragments between two cuts are read from both ends. A fragment shorter than about 150 bases comes from open DNA; longer ones span one, two or three nucleosomes. /data/examples/atacseq/ has the 5,771 fragments that fall on chromosome I of baker's yeast (230,218 bases), of MACS's own test of HMMRATAC – a yeast ATAC-seq (SRA run SRR1822137).

zcat /data/examples/atacseq/yeast_SRR1822137_chrI.bedpe.gz | awk '{ n[$3 - $2 < 150 ? "shorter than 150" : "150 or longer"]++ } END { for (k in n) print k, n[k] }'
macs3 hmmratac -i /data/examples/atacseq/yeast_SRR1822137_chrI.bedpe.gz -f BEDPE --randomSeed 12345 -n atac --outdir atac
head -3 atac/atac_accessible_regions.narrowPeak
wc -l atac/atac_accessible_regions.narrowPeak

HMMRATAC learns the lengths of four kinds of fragment – free of nucleosomes, and spanning one, two or three – from the data (its messages: means of 50, 174, 357 and 549 bases); it splits the signal into those four, trains a hidden Markov model with three states – open, nucleosome, background – on regions with strong signal, and decodes the chromosome with it. The accessible regions are the open stretches and their flanking nucleosomes: 56 here (atac_accessible_regions.narrowPeak, with the summit). atac_model.json is the model, which --model uses again on other data. It picks fragments at random to learn their lengths: --randomSeed makes a run repeatable. With so few fragments the model depends on that draw – try another seed – while a real experiment has millions. The positions are those of the yeast genome sacCer3:

igv --genome sacCer3 atac/atac_accessible_regions.narrowPeak --locus chrI:120,000-150,000

macs3 has fourteen commands; macs3 COMMAND -h lists the options of each. Besides callpeak, hmmratac and bdgcmp: callpeak --broad joins nearby peaks into broad regions, for histone marks such as H3K27me3 (NAME_peaks.broadPeak and .gappedPeak); bdgpeakcall and bdgbroadcall call peaks from a bedGraph track; bdgdiff finds the differences between two conditions; cmbreps combines replicates; bdgopt changes a track; filterdup removes duplicate reads; predictd works out only the fragment length; pileup piles up the reads; randsample takes a sample of them; refinepeak refines the summits; callvar calls variants in the peaks – from reads of a ChIP and its control. In the Galaxy tab, the tool MACS3 callpeak (section Peak calling) runs macs3 callpeak on BAM or BED datasets.

MACS3 3.0.5 is compiled to WebAssembly from its released source – its C and Cython parts, with those of hmmlearn – by the genomics-wasm project, and gives the same peaks, tracks and tables as MACS3 on Linux, byte for byte – its models the same but for the last decimals (see Programs and versions). It runs as one process: callvar -m N works through its N parts one after another, with the same result. It keeps all the reads in memory: a chromosome, or the genome of yeast, is fine; the tens of millions of reads of a whole human genome need more memory than a browser tab has (see Limits).

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

MACS3: MACS3's README asks for the 2026 paper when MACS 2 or 3 is used, and for the 2008 paper when MACS 1 is; the 2008 paper describes the shifting model and the local lambda that MACS3 still uses. For macs3 hmmratac, the paper of HMMRATAC too.; Doherty P, Hu Q, Zhuang Z, Zhang H, Penikalapati SC, Liu S, et al. (2026) MACS3: A Peak-calling Platform for Bulk and Single-cell Regulatory Genomics. Genomics, Proteomics & Bioinformatics qzag097. doi:10.1093/gpbjnl/qzag097; Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. (2008) Model-based Analysis of ChIP-Seq (MACS). Genome Biology 9(9): R137. doi:10.1186/gb-2008-9-9-r137
HMMRATAC (macs3 hmmratac): Tarbell ED, Liu T (2019) HMMRATAC: a Hidden Markov ModeleR for ATAC-seq. Nucleic Acids Research 47(16): e91. doi:10.1093/nar/gkz533
igv.js (the genome browser): Robinson JT, Thorvaldsdottir H, Turner D, Mesirov JP (2023) igv.js: an embeddable JavaScript implementation of the Integrative Genomics Viewer (IGV). Bioinformatics 39(1): btac830. doi:10.1093/bioinformatics/btac830

Intervals and annotation

BED and GFF files describe regions: genes, exons, peaks. bedtools compares and combines them; gffread converts annotations and writes the sequences of transcripts and proteins.

The examples work in a folder of your own:

The files

head -3 /data/examples/annotation/cyp2c_exons.bed
grep -v "^#" /data/examples/annotation/cyp2c_slices.gff3 | cut -f 1-5,9 | head -4

The BED file lists the exons of eight transcripts (named TRANSCRIPT_exonN): those of CYP2C19 and CYP2C9, and of two pseudogenes that lie in their introns, MTND4P19 and MTND4P20. The GFF3 file describes the same genes in full: genes, transcripts, exons and the coding parts (CDS).

BED counts from 0 and its end is not included (a region 100–200 is bases 101 to 200); GFF and VCF count from 1. bedtools takes care of that when it compares them.

bedtools intersect: what overlaps what

# the variants that fall in an exon (of any of the transcripts)
bedtools intersect -header -u -a /data/examples/variants/NA12878.calls.vcf.gz -b /data/examples/annotation/cyp2c_exons.bed | grep -v "^##" | cut -f 1-5
# and those that do not
bedtools intersect -v -a /data/examples/variants/NA12878.calls.vcf.gz -b /data/examples/annotation/cyp2c_exons.bed | wc -l

-u: each feature of A that overlaps B, once; -v: those that do not; -wa -wb: both features of each overlap; -header keeps the header of A (needed for a valid VCF).

Coverage, merging, sequences

minimap2 -ax sr /data/examples/reads/reference.fa /data/examples/reads/NA12878_R1.fastq /data/examples/reads/NA12878_R2.fastq 2> /dev/null | samtools sort -o mapped.bam -
# depth along the reference, as a BedGraph (positions with the same depth form one line)
bedtools genomecov -ibam mapped.bam -bg > coverage.bedgraph
head -3 coverage.bedgraph
# the exons of all eight transcripts, joined where they overlap or touch
sort -k1,1 -k2,2n /data/examples/annotation/cyp2c_exons.bed | bedtools merge -i - | head -3
# the DNA of the exons
bedtools getfasta -fi /data/examples/reads/reference.fa -bed /data/examples/annotation/cyp2c_exons.bed -name | head -2

Also: bedtools coverage (how much of each region is covered), bedtools closest, bedtools complement, bedtools slop (widen regions; it needs a genome file of sequence lengths: cut -f 1,2 reference.fa.fai). bedtools alone lists them all.

gffread: annotations and their sequences

# GFF3 to GTF
gffread /data/examples/annotation/cyp2c_slices.gff3 -T -o cyp2c.gtf
head -3 cyp2c.gtf
# the spliced transcripts (-w) and their proteins (-y)
gffread /data/examples/annotation/cyp2c_slices.gff3 -g /data/examples/reads/reference.fa -w transcripts.fa -y proteins.fa
grep ">" proteins.fa

Other options: -x (coding sequences only), -E (report problems in the file), --table. The proteins can go straight to an alignment (Alignments and trees) or to BLAST (BLAST).

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

bedtools: Quinlan AR, Hall IM (2010) BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26(6): 841–842. doi:10.1093/bioinformatics/btq033
gffread: Pertea G, Pertea M (2020) GFF Utilities: GffRead and GffCompare. F1000Research 9: 304. doi:10.12688/f1000research.23297.2

Alignments and trees

Pairwise alignment (needleman_wunsch, smith_waterman), multiple alignment (kalign, muscle), trees (FastTree), and long DNA alignments (lastz).

The examples work in a folder of your own:

Two sequences

needleman_wunsch --printscores GATTACA GCATGCT
smith_waterman --printseq ACACACTA AGCACACA

needleman_wunsch aligns two sequences from end to end (global); smith_waterman finds the best-matching parts (local) and lists every good hit. The scores they use unless told otherwise differ: needleman_wunsch --match 1 --mismatch -2 --gapopen -4 --gapextend -1, smith_waterman --match 2 --mismatch -2 --gapopen -2 --gapextend -1. For proteins --scoring BLOSUM62 (also PAM30, PAM70, BLOSUM80). With --file pairs.fasta they align the sequences of a file two at a time.

needleman_wunsch --scoring BLOSUM62 --gapopen -10 --gapextend -1 --printscores \
  MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHF MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESF

lcs SEQUENCE (also from seq-align) lists the substrings that occur twice in a sequence, longest first.

Multiple alignment

kalign -i /data/examples/sequences/myoglobins.fasta -o myoglobins.afa
muscle -in /data/examples/sequences/myoglobins.fasta -out myoglobins_muscle.afa
open myoglobins.afa

Both write aligned FASTA: every sequence has the same length, with - for gaps. The Files tab shows an alignment in colour, column by column. Clustal format: kalign -f clu, muscle -clw. Give kalign its output file with -o: this version also prints its banner on the standard output, so kalign … > file would put the banner into the file.

A tree

fasttree myoglobins.afa > myoglobins.nwk
cat myoglobins.nwk
open myoglobins.nwk

FastTree builds an approximately-maximum-likelihood tree (JTT+CAT for proteins; -nt for DNA, with -gtr for the GTR model). The numbers before the branch lengths are local support values (0–1). The Files tab draws a .nwk file as a tree; in R, ape::read.tree() reads it, in Python Bio.Phylo.

FastTree -nt -gtr /data/examples/sequences/tRNA-Thr.afa > trna.nwk 2> /dev/null && cat trna.nwk

Long DNA: lastz

lastz aligns long DNA sequences – a region against a chromosome, two related genomes. [multiple] after the target file means it has more than one sequence; quote it, because [ ] mean something to the shell.

samtools faidx /data/examples/reads/reference.fa human_CYP2C19:1000-3000 > piece.fa
lastz '/data/examples/reads/reference.fa[multiple]' piece.fa --format=general | head -4

Formats: --format=general (a table), maf, sam, axt; --notransition --step=20 make it faster on big sequences.

Alignments and trees are also tools in Galaxy (Kalign, MUSCLE, FastTree), and BLAST searches run at EBI from the Databases tab.

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

seq-align (needleman_wunsch, smith_waterman): seq-align has no paper; cite the algorithms it carries out.; Turner I. seq-align: Smith-Waterman and Needleman-Wunsch alignment in C. github.com/noporpoise/seq-align; Needleman SB, Wunsch CD (1970) A general method applicable to the search for similarities in the amino acid sequence of two proteins. Journal of Molecular Biology 48(3): 443–453. doi:10.1016/0022-2836(70)90057-4; Smith TF, Waterman MS (1981) Identification of common molecular subsequences. Journal of Molecular Biology 147(1): 195–197. doi:10.1016/0022-2836(81)90087-5
Kalign: Lassmann T (2020) Kalign 3: multiple sequence alignment of large datasets. Bioinformatics 36(6): 1928–1929. doi:10.1093/bioinformatics/btz795
MUSCLE: Edgar RC (2004) MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Research 32(5): 1792–1797. doi:10.1093/nar/gkh340
FastTree: Price MN, Dehal PS, Arkin AP (2010) FastTree 2 – Approximately Maximum-Likelihood Trees for Large Alignments. PLoS ONE 5(3): e9490. doi:10.1371/journal.pone.0009490
LASTZ: Harris RS (2007) Improved pairwise alignment of genomic DNA. PhD thesis, The Pennsylvania State University. github.com/lastz/lastz

Text tools

Most bioinformatics files are text. The GNU tools – grep, sed, awk, sort, cut, uniq, wc – answer many questions in one line, without a program of their own.

The examples work in a folder of your own:

Look and count

head -3 /data/examples/tables/yeast_genes.csv | cut -c 1-120
wc -l /data/examples/tables/yeast_genes.csv
grep -c ">" /data/examples/sequences/myoglobins.fasta

Columns: cut, sort, uniq

# essential genes (column 3 of the yeast table: 1 = essential)
cut -d, -f3 /data/examples/tables/yeast_genes.csv | tail -n +2 | sort | uniq -c
# the variant positions of a VCF, as a list
bcftools view -H /data/examples/variants/NA12878.calls.vcf.gz | cut -f 1,2 | head -5

Careful with cut -d, on CSV files: a field with a comma inside quotes is cut there too. In the yeast table the location ("Golgi, vesicles & vacuole") and the description have such commas, so every column after the location – the chromosome is column 16 – would be counted wrong. gawk can read the quoted fields (FPAT describes what a field looks like); pandas and R read CSV files properly too:

# genes per chromosome (column 16), most first
awk -v FPAT='[^,]*|"[^"]*"' 'NR > 1 { print $16 }' /data/examples/tables/yeast_genes.csv | sort | uniq -c | sort -rn | head -5

grep: lines that match

grep -i "kinase" /data/examples/tables/yeast_genes.csv | wc -l
grep -v "^#" /data/examples/annotation/cyp2c_slices.gff3 | cut -f 3 | sort | uniq -c

-i ignore case, -v lines that do not match, -w whole words, -c count, -E extended regular expressions, -o only the matching part, -P Perl patterns.

awk: tables and arithmetic

# mean quality of the variant calls
bcftools view -H /data/examples/variants/NA12878.calls.vcf.gz | awk '{ s += $6; n++ } END { printf "%d calls, mean QUAL %.1f\n", n, s / n }'
# the gene lines of a GFF file: column 3 is "gene" (grep -w gene would find the mRNAs too – their
# column 9 says Parent=gene:…)
awk -F'\t' '$3 == "gene"' /data/examples/annotation/cyp2c_slices.gff3 | cut -f 1,4,5,9 | cut -c 1-100
# GC content of the reads (the sequence is every 4th line, from the 2nd; gsub counts what it replaces)
awk 'NR % 4 == 2 { n += length($0); gc += gsub(/[GC]/, "") } END { printf "%d reads, GC %.1f%%\n", NR / 4, 100 * gc / n }' /data/examples/reads/NA12878_R1.fastq

sed: change text

# shorter sequence names
sed 's/^>MyG_/>/' /data/examples/sequences/myoglobins.fasta | grep ">"
# lines 5 to 8 of a file
sed -n '5,8p' /data/examples/reads/NA12878_R1.fastq

sed -i changes a file in place (on your own copy – the examples are read-only).

JSON and the web: curl, jq

curl -s "https://rest.ensembl.org/lookup/symbol/homo_sapiens/CYP2C19?content-type=application/json" | jq '{id, start, end, strand, biotype}'
curl -s -o P02185.fasta "https://rest.uniprot.org/uniprotkb/P02185.fasta" && head -2 P02185.fasta

curl and wget reach web services that allow web pages to use them (Ensembl, UniProt, NCBI E-utilities, EBI, the PDB, …). Others fail with a message about the server; download such files with your browser and upload them in the Files tab.

Compressed files and archives

gzip -c /data/examples/sequences/pBR322.fasta > pBR322.fasta.gz
zcat pBR322.fasta.gz | head -2
tar -czf examples.tar.gz -C /data/examples sequences
tar -tzf examples.tar.gz
zip -r sequences.zip /data/examples/sequences > /dev/null && unzip -l sequences.zip

gzip, gunzip and zcat are carried out by bgzip; zip, unzip and tar (with z for .tar.gz) work on your files.

Scripts and pipelines

Commands that work once belong in a script – so that they run again tomorrow, on new data, and show exactly what was done.

The examples work in a folder of your own:

A bash script

Write the commands into a file – in the Files tab, or with a here-document as below – and run it with bash:

cat > count_reads.sh << 'EOF'
#!/usr/bin/env bash
set -euo pipefail          # stop at the first error
for f in "$@"; do
  n=$(( $(wc -l < "$f") / 4 ))
  echo -e "$(basename "$f")\t$n"
done
EOF
bash count_reads.sh /data/examples/reads/*.fastq

"$@" are the file names given to the script; set -euo pipefail makes it stop at the first command that fails, instead of carrying on with half a result. After chmod +x count_reads.sh it also starts as ./count_reads.sh.

A whole analysis in one script

cat > call_variants.sh << 'EOF'
#!/usr/bin/env bash
set -euo pipefail
ref=/data/examples/reads/reference.fa
r1=/data/examples/reads/NA12878_R1.fastq
r2=/data/examples/reads/NA12878_R2.fastq
fastp -i "$r1" -I "$r2" -o trimmed_1.fastq -O trimmed_2.fastq -h fastp.html -j fastp.json 2> fastp.log
minimap2 -ax sr -R '@RG\tID:NA12878\tSM:NA12878' "$ref" trimmed_1.fastq trimmed_2.fastq 2> minimap2.log | samtools sort -o mapped.bam -
samtools index mapped.bam
bcftools mpileup -f "$ref" mapped.bam 2> /dev/null | bcftools call -mv -Oz -o calls.vcf.gz
bcftools filter -i 'QUAL>=20 && INFO/DP>=5' -Oz -o filtered.vcf.gz calls.vcf.gz
echo "variants: $(bcftools view -H filtered.vcf.gz | wc -l)"
EOF
bash call_variants.sh

Python and R scripts

cat > gc.py << 'EOF'
import sys
from Bio import SeqIO
from Bio.SeqUtils import gc_fraction
for rec in SeqIO.parse(sys.argv[1], "fasta"):
    print(rec.id, f"{gc_fraction(rec.seq):.3f}", sep="\t")
EOF
python gc.py /data/examples/reads/reference.fa
Rscript -e 'x <- read.csv("/data/examples/tables/yeast_genes.csv"); print(table(x$essential))'

Both use the same Python and R as the notebooks. A script's plots go to files: plt.savefig("plot.png") in Python; in R, png("plot.png") … dev.off() (otherwise Rplots.pdf).

Snakemake

Snakemake runs a pipeline from rules: each rule says which files it makes from which, and Snakemake runs only what is missing or out of date. The terminal has a browser edition of Snakemake 9 – snakemake -n shows the plan, snakemake --cores 1 runs it, snakemake --dag | dot -Tsvg -o dag.svg draws it:

mkdir -p ~/smk && cd ~/smk
cat > Snakefile << 'EOF'
rule all:
    input: "results/flagstat.txt"

rule map:
    input:
        ref="/data/examples/reads/reference.fa",
        r1="/data/examples/reads/NA12878_R1.fastq",
        r2="/data/examples/reads/NA12878_R2.fastq"
    output: "results/mapped.bam"
    shell: "minimap2 -ax sr {input.ref} {input.r1} {input.r2} | samtools sort -o {output} -"

rule flagstat:
    input: "results/mapped.bam"
    output: "results/flagstat.txt"
    shell: "samtools flagstat {input} > {output}"
EOF
snakemake -n
snakemake --cores 1
cat results/flagstat.txt

Run snakemake --cores 1 again: nothing to do, because every result is newer than its inputs. Snakemake keeps a record of how each file was made: snakemake --summary.

Pipelines written for a computer or a cluster often give each rule a conda environment (conda: "envs/mapping.yaml") and are run with snakemake --use-conda. They run here as they are. Every program is installed already, so no environments are made: --use-conda only checks that each environment file is there and can be read.

Graphviz

dot draws graphs written in the DOT language – workflows, networks, decision trees:

echo 'digraph { rankdir=LR; reads -> fastp -> minimap2 -> bcftools -> variants }' | dot -Tsvg -o steps.svg
open steps.svg

Keep a record

Save session above the terminal turns the commands you typed into a notebook or a script (see Saving a session) – a good start for a script of your own.

References

The papers to cite when you use this program in your work – with the version you used (Programs and versions):

Snakemake: Mölder F, Jablonski KP, Letcher B, Hall MB, van Dyken PC, Tomkins-Tinch CH, et al. (2025) Sustainable data analysis with Snakemake. F1000Research 10: 33. doi:10.12688/f1000research.29032.3

Your files

Your files are in your home folder, ~ (/home/user). The notebooks, the terminal, Galaxy and the Files tab all use the same files.

Where your files are kept

There is no server: your files are kept in this browser, on this computer (in its storage for web sites, IndexedDB). They are there when you come back – in the same browser, on the same computer, at the same web address.

Download what you want to keep

Another browser, another computer or a private window does not have your files.
Clearing the browser's data for this site (cookies and site data) deletes them – some computers do that by themselves when you log out (shared university PCs often do).
Browsers may delete site data when the disk is nearly full.

My files → Download all my files (.zip) (top right) keeps everything; download single notebooks with Download in the Notebook tab.

The dot at the top right, Saved in this browser, shows that the latest changes are stored. My files also says how much space you use; browsers give a site a share of the free disk space – usually several gigabytes.

The Files tab

Open a file by clicking it. Text files open in an editor (with colours for Python, R, shell, YAML …): edit, then Save or Ctrl+S. Several files open as tabs.
Other files are shown in a way that fits them: tables (CSV, TSV) as tables, FASTA files with their sequences, alignments in colour, trees (.nwk) drawn, pictures, PDF files (with zoom and page numbers; their text can be selected and copied), HTML reports (without their scripts). Open as text shows the text itself.
Notebooks (.ipynb) open in the Notebook tab, structures (.pdb, .cif) in the 3D viewer, networks (.sif, .graphml, .xgmml, .gml, .cyjs, .cx2) in the Networks tab.
The ⋯ of a file or folder: open in another tab (Notebook, 3D viewer, Genome browser, Networks, Galaxy), download, rename, duplicate, copy the path, open the folder in the terminal, unzip, delete. A table (CSV, TSV) opens in Networks as a table of edges, or as values for the nodes of the network that is open.
The buttons at the top: upload files from your computer, a new file, a new folder, show hidden files (names that start with a dot).

Upload and download

Upload: drag files from your computer onto the Files tab, or use the upload button or My files → Add files from my computer…. Folders can be uploaded too. A .zip can be unpacked with its ⋯ menu (or unzip in the terminal).
Download: the ⋯ of a file; a folder downloads as a .zip. In the terminal: download FILE.
Very large files (several gigabytes) do not fit: work on a part of them (head, samtools view -s, seqtk sample) on your own computer first.

The example data

/data/examples comes with the workbench and cannot be changed. The Files tab shows it under Examples (read-only); its README.md says where each file comes from. To change a file, copy it first – its ⋯ menu has Copy to my files, or in the terminal:

cp -r /data/examples/reads ~/example-reads
ls -l ~/example-reads

Folder	What is in it
`reads/`	Illumina exome reads of the reference sample NA12878 (3,519 pairs, 76 bases) from two regions of chromosome 10, and their reference (`reference.fa`: the CYP2C19 and CYP2C9 regions of GRCh37)
`variants/`	the variants called from those reads (`NA12878.calls.vcf.gz`)
`annotation/`	the genes and exons of the two regions (GFF3, BED)
`sequences/`	pBR322 (DNA), myoglobins of eight animals (protein), aligned tRNAs
`tables/`	5,481 yeast genes in 19 columns – names, essentiality, protein length and abundance, location, chromosome … (`yeast_genes.csv`; the README describes each column)
`structures/`	sperm whale myoglobin, PDB 1MBN

Start afresh

My files → Start afresh deletes all your files, Galaxy's history and the workbench's settings in this browser – after you have downloaded what you want to keep. The settings of the AI assistant and of the help pane are kept.

Galaxy

Galaxy runs analysis tools from forms, keeps every result in a history with a record of how it was made, and joins tools into workflows. This one works like usegalaxy.org and its sister servers – and runs in your browser, with the same programs as the terminal.

The three columns

Tools (left): Upload Data, and the tools in sections. Search them with the box above the list.
The middle: the form of a tool, a dataset, the details of a job, or the workflow editor.
History (right): your datasets, newest at the top, numbered. Grey: waiting; yellow: running; green: ready; red: failed. Click a dataset to open it.

Hide the help pane (Ctrl+\) to give Galaxy room.

Get data

Upload Data has three sources:

From your files – the files in your home folder and the example data; tick them and press Start. The Files tab sends files here too (⋯ → Upload to Galaxy).
From your computer – choose files or drop them.
Paste text – a few sequences, a list of positions …

Galaxy decides the datatype from the name and content (fastqsanger for reads, fasta, bam, vcf, bed, tabular, hmm3 for profile HMMs, bus for the records of kallisto bus …); tools offer only the datasets of the types they read. If it guessed wrong, choose the type in the upload form. Compressed text files (.fastq.gz, .fa.gz, .hmm.gz …) are uncompressed on the way in.

Run a tool

Click a tool. Its form lists its inputs – choose datasets from your history – and its settings, with Galaxy's own help under the form and, at its end, the papers to cite when you use the tool (References), as on a Galaxy server.
Run tool. The new datasets appear at the top of the history; wait for them to turn green.
Open a dataset: shows it (tables as tables, HTML reports, trees drawn), downloads it, copy to my files puts it in ~/galaxy for the terminal and notebooks. BAM, VCF and BED datasets have to show them in the Genome browser; networks (the datatypes sif and xgmml) have to open them in the Networks tab.

A first analysis with the example reads: upload NA12878_R1.fastq, NA12878_R2.fastq and reference.fa; run fastp on the two read files; Map with minimap2 with the reference and the two trimmed outputs; bcftools mpileup on the BAM, then bcftools call, then bcftools filter. Each step takes seconds.

The tools

Section	Tools
Quality control	fastp, Cutadapt, seqtk_comp
Mapping	Map with minimap2, Bowtie2, Map with BWA-MEM, HISAT2 (RNA-seq reads, across introns), Map with BWA-MEM2 – each single- or paired-end
SAM/BAM	Samtools flagstat, stats, idxstats, depth, view (filter); mosdepth
RNA-seq	featureCounts, StringTie, Kallisto quant, Salmon quant
Single-cell	bustools sort, correct, count, inspect, text (the BUS records of `kallisto bus` to a matrix of cells and genes)
Peak calling	MACS3 callpeak (the peaks of ChIP-seq, with or without a control)
Variant calling	bcftools mpileup, call, filter, view, query, stats; FreeBayes; Insert indel qualities and Call variants (LoFreq); Strelka Somatic (tumour and normal)
Genotypes (GWAS)	PLINK 2: quality control, filters, LD pruning, kinship, association tests and scores on a VCF of many samples
Metagenomics	Kraken2 (with the example database), Bracken
Annotation	Prodigal Gene Predictor (genes of bacteria and archaea); hmmbuild, hmmsearch, hmmscan, phmmer, nhmmer, hmmalign (HMMER: profile HMMs)
Intervals	bedtools Genome Coverage, Intersect intervals
FASTA/FASTQ	seqtk_seq (FASTQ to FASTA), seqtk_sample
Alignments and trees	Kalign, MUSCLE, FastTree
Text manipulation	Select lines, Cut columns, Sort, Line/Word/Character count, Concatenate

The tool IDs and versions are those of the Galaxy Tool Shed, so a workflow made here opens on usegalaxy.org (see Workflows) – except Kalign and MUSCLE, the five bustools tools, MACS3 callpeak and PLINK 2: the Tool Shed has no tool for these programs or versions (its plink tool runs PLINK 1.9, its MACS2 callpeak MACS2 – the form of MACS3 callpeak here is that one's), so they are the workbench’s own, and a server shows those steps as missing. Kallisto quant and Salmon quant carry the IDs of the newest tools of the Tool Shed, which run kallisto 0.48.0 and Salmon 1.10.1; here they run the terminal's kallisto 0.52.0 and Salmon 2.8.0, so a server's results can differ a little. Call variants with LoFreq runs lofreq call on one core where a server runs lofreq call-parallel. Kraken2 carries the ID of the Tool Shed's newest Kraken2 tool, which runs Kraken2 2.17.1; here it runs the terminal's 2.17.2. Its database is the example one: a server offers the databases its administrators installed – as for hmmscan, whose built-in database here holds the three profiles of the HMMER tutorial, where a server has Pfam. The HMMER tools run on one thread, without --cpu; Galaxy's tools pass the three filter thresholds with --max, which HMMER refuses, so here --max goes alone. Prodigal's form has no SQN format, which Prodigal 2.6.3 does not write. featureCounts, Salmon quant and PLINK 2 work on several threads: they need the page to be cross-origin isolated (see Limits).

How was this made? Dataset details

The of a dataset shows Galaxy's record of it: the tool and its version, the versions of the programs (dependencies), the command line that ran, every parameter, and the program's messages. That record is why Galaxy analyses can be repeated – and the first place to look when a dataset is red. opens the tool again with the same inputs and settings.

The history

Click the history's name to rename it.
× deletes a dataset; the history's ⋯ shows deleted datasets (to undelete them) and purges them, which frees their space.
⋯ → Open the history as a notebook (or download it as .ipynb): every upload and job as a %%bash cell with its command, inputs and parameters. Run the notebook and it makes the same datasets again, as files data1.fastq, data2.bam …
⋯ → Start a new history clears it (your workflows stay).

The history and its datasets are kept in this browser with your files (see where your files are kept).

Different from a Galaxy server

The tools are a selection, run with the programs of the terminal, one job at a time.
There is one history and no user account; nothing is shared with others.
Datasets have no genome build ("database: ?") unless a tool passes one on.

References

The papers to cite when you use this program in your work – with the version you used (Programs and versions):

Galaxy: The Galaxy Community (2024) The Galaxy platform for accessible, reproducible, and collaborative data analyses: 2024 update. Nucleic Acids Research 52(W1): W83–W94. doi:10.1093/nar/gkae410

Galaxy workflows

A workflow is a saved graph of tools: each step's inputs are connected to the outputs of earlier steps. It runs again with one click – on the same data or on new data.

Extract a workflow from your history

After a few tools have run, Extract workflow (at the top of the history) turns the jobs into a workflow: the datasets you uploaded become its inputs, and each job a step. It opens in the editor.

Build one in the editor

Workflows → Create a new workflow. Give it a name (the box at the top).
+ Input dataset adds an input; click it to set its label and datatype (for example fastqsanger for reads).
Click tools on the left: each becomes a box. Its inputs are on the left side, its outputs on the right.
Drag from an output dot to an input dot to connect them. Only matching datatypes connect: while you drag, the inputs that fit light up.
Click a box to change its settings (right column). With nothing selected, the right column lists what is still missing.

Click a connection to remove it; × on a box removes the step. Tidy layout arranges the boxes. Changes are saved as you go.

Run it

Run asks for a dataset of your history for each input, lists the steps, and runs them in order. Invocations shows each run and the state of its jobs; the outputs appear in your history, with the full record of each job ().

Take it to a Galaxy server

.ga downloads the workflow in Galaxy's own format. On usegalaxy.org (or usegalaxy.eu, .org.au …): Workflows → Import, choose the file. The steps name the Tool Shed tools and versions they use; where the server has another version of a tool, it says so when the workflow opens. Steps with Kalign or MUSCLE come up as missing tools (they have no Tool Shed counterpart): replace them with the server's own alignment tool. to my files keeps the .ga file in ~/galaxy/workflows, and Import a .ga file opens a workflow made elsewhere – if its tools are among the ones here.

Workflows and notebooks

A workflow records the steps; a history records one run of them. To keep the commands of a run as code, use the history's ⋯ → Open the history as a notebook (see Galaxy).

Databases

The Databases tab asks public bioinformatics databases your questions, straight from the browser: proteins and structures, genes and genomes, variants, interactions, pathways, papers and small molecules. Every answer can be kept – saved among your files, or sent to a notebook as Python or R code that asks the same question again.

The cards

Group	Card	Ask it for
Proteins	UniProt	protein entries: function, names, location, disease and variants, domains, sequence
	Protein Data Bank	experimental 3D structures
	AlphaFold DB	the predicted structure of a protein, with its confidence
	InterPro	the families, domains and sites of a protein
Genes and genomes	Ensembl	genes and transcripts, sequences, the effect of a variant (VEP), orthologues, the genes in a region
	NCBI	Gene, Nucleotide, Protein, PubMed, dbSNP, ClinVar, Taxonomy, Assembly; sequences as FASTA or GenBank
	UCSC Genome Browser	the DNA of any region of a genome
	European Nucleotide Archive	sequencing runs and their FASTQ files; any sequence by accession
Sequence searches	BLAST	similar sequences in UniProt, the PDB, AlphaFold DB or ENA
Interactions and networks	STRING	protein–protein interaction networks; the functions a list of genes shares
	IntAct	molecular interactions shown by experiments, with their papers
	NDEx	public networks: pathways, interactomes, the networks of papers
Function and pathways	Reactome	the pathways of a protein, with diagrams
	Gene Ontology	the GO terms of a protein; terms of GO and other ontologies
	g:Profiler	the GO terms and pathways a list of genes shares (enrichment)
Literature and chemistry	Europe PMC	papers, with links to free full texts
Literature and chemistry	PubChem	drugs, metabolites and ligands: properties, drawings, 3D models

Ask a question

Click a card. Cards that answer several kinds of question have tabs at the top – Ensembl has Gene, Sequence, Variant effect (VEP), Orthologues and Genes in a region.
Fill in the form, or click one of the examples after Try: – they fill in the form and search.
Search. The answer appears under the form: a title and a short summary, then a table (click a column heading to sort by it), a picture or a sequence.

← All databases goes back to the cards. Help opens the card's page in this pane; About the API opens the service's own documentation.

Keep the answer

The buttons above an answer:

Save to my files (a sequence, a report) and Save the table (.tsv) – it asks where, in ~/data unless you change it; structures go to ~/structures. From there the terminal, the notebooks and Galaxy can use them.
Send to notebook → Python (requests) or R: adds a cell with code that makes the same query to the notebook that is open. Run it, and the answer is in Python (a dictionary, or a pandas table) or in R – ready for the next step of an analysis.
Raw data (JSON): what the service sent back, unchanged – every field, also those the table leaves out.
Copy the URL: the query as a web address. Paste it into a browser, into curl, or into code.
Some answers have more: Show in the 3D viewer (AlphaFold, PubChem), Show in the Genome browser (Ensembl genes), Open in Networks (STRING, IntAct, NDEx), a FASTA file for each transcript (Ensembl), Open at … the service's own page.

Your queries – as a notebook

Under the cards, Your queries lists the questions you asked in this visit, newest first – those of the UniProt and PDB cards too. Save as a Python notebook (or … as an R notebook) writes them all into one notebook, ~/notebooks/database-queries-DATE.ipynb: for each query a heading, its URL and a code cell that makes it again. It is a record of where your data came from; download it from the Notebook tab like any other notebook.

The list is kept until you reload the page

Save it as a notebook before you leave.

Asking from code

Each card uses a REST API: a web address goes in, data comes out – JSON for most, FASTA or text for sequences. The same works in a notebook:

import requests

r = requests.get("https://rest.ensembl.org/lookup/symbol/homo_sapiens/CYP2C19",
                 headers={"Content-Type": "application/json"})
r.raise_for_status()          # stop here if the service said no
gene = r.json()
gene["id"], gene["seq_region_name"], gene["start"], gene["end"]

library(jsonlite)
gene <- fromJSON("https://rest.ensembl.org/lookup/symbol/homo_sapiens/CYP2C19?content-type=application/json")
c(gene$id, gene$seq_region_name, gene$start, gene$end)

In the browser, R can only make GET requests (a URL); for services that need a POST request – BLAST, VEP with positions, the PDB search – use Python. More in the recipe Databases from code.

Good practice

Databases change. New releases add, correct and remove entries; a query can give another answer next month. Note the date – the notebook of your queries does – and, for positions on a genome, the assembly: GRCh38 (hg38) and GRCh37 (hg19) number the same base differently.
Annotations are predictions. VEP calls rs4244285 synonymous – yet it is CYP2C19*2, a loss-of-function allele that changes the splicing of the gene (see Genes, genomes and variants). Look at the clinical significance and the literature before you conclude.
The services are shared by scientists everywhere. Ask for what you need, not thousands of queries in a quick loop: Ensembl allows about 15 requests a second, NCBI 3. A request that gets no answer at all is sent again (twice) by the tab; for code, see the recipe.
When a query fails, the message says why: Not found (404) – check the identifier; did not accept the query (400) – check what you typed; slow down (429) – wait a minute; not answering (5xx) or no answer that the page could read – the service is busy or down: try again later.

Public services

Only what you type into a form is sent to the service (for BLAST, also your e-mail address). These are public services that keep logs: never type or paste personal or patient data into them.

UniProt

UniProtKB, the UniProt Knowledgebase, is the reference collection of protein sequences and of what is known about each protein. The UniProt card shows its entries laid out like the pages of uniprot.org, from live data.

Reviewed and unreviewed

UniProtKB has two parts. Swiss-Prot (★ Reviewed) entries are read and written by expert curators, from the literature; there are about 575,000 of them. TrEMBL (Unreviewed) entries – about 150 million – are annotated automatically, mostly from genome sequences (numbers of release 2026_03). When a reviewed entry exists, use it.

Each entry has an accession (P04637), which never changes and is the one to cite, and an entry name (P53_HUMAN). Its annotation score (1 to 5) says how much is known about it, not how correct it is.

Search

Type a protein name, a gene name, an accession or any words, and Search. Tick Reviewed (Swiss-Prot) only and Human (Homo sapiens) only to narrow the list. An accession on its own opens the entry. The results come 25 at a time (Show 25 more); click an accession to open it.

UniProt's query language works in the box:

Query	Finds
`gene:CYP2C19 AND organism_id:9606`	the human entries of a gene
`reviewed:true AND length:[100 TO 120]`	reviewed proteins of 100 to 120 amino acids
`ec:3.4.21.4`	proteins with that enzyme activity (EC number)
`keyword:KW-0675 AND organism_id:9606`	human proteins with the keyword Receptor
`cc_disease:Alzheimer AND reviewed:true`	reviewed entries that mention the disease

An entry

The top of an entry gives the protein's name, gene, organism, status, annotation score and length, with buttons: AlphaFold model in 3D, Save FASTA (to your files) and Send to notebook (Python or R code that fetches the entry). The sections are listed on the left:

Section	What it says
Function	what the protein does, catalytic activity, GO terms
Names & Taxonomy	the names of the protein and gene, the organism and its lineage
Subcellular Location	where in the cell, and the membrane topology
Disease & Variants	diseases, natural variants, mutagenesis experiments
PTM / Processing	signal peptides, chains, modified residues, disulfide bonds, glycosylation
Expression	tissues and developmental stages
Interaction	subunits and binary interactions (IntAct); the protein's STRING network or its IntAct interactions in the Networks tab
Structure	the PDB entries (which parts of the protein each covers) and the AlphaFold model
Family & Domains	domains, regions, motifs, repeats, and the family
Sequence & Isoforms	the sequence – find a motif, go to a residue, copy or save it as FASTA – and how the other isoforms differ

Many statements carry an evidence tag: experimental evidence, or inferred by similarity to another protein (ECO codes). Positions are numbered along the canonical sequence of the entry.

Features in 3D

Domains, sites, variants and other features are drawn along the sequence. Click one to see where it is in the 3D viewer: the positions are found on the structure that is loaded – a PDB entry of this protein, or its AlphaFold model – by aligning its sequence with the entry's. If no structure of the protein is loaded, or the residues are not in it (experimental structures often contain only part of a protein), it offers to load the AlphaFold model, which covers the whole sequence.

In Structure, the coverage map shows which residues each PDB entry contains; click a bar, or 3D in the table, to load that structure, and Details to open it in the Protein Data Bank card.

From code

An entry as JSON is https://rest.uniprot.org/uniprotkb/P04637.json, its sequence …/P04637.fasta, and a search https://rest.uniprot.org/uniprotkb/search?query=…. Send to notebook writes the code; or:

import requests
import pandas as pd

r = requests.get("https://rest.uniprot.org/uniprotkb/search",
                 params={"query": "gene:CYP2C19 AND reviewed:true", "fields": "accession,id,organism_name,length", "format": "tsv"})
from io import StringIO
pd.read_csv(StringIO(r.text), sep="\t")

Protein Data Bank

The Protein Data Bank (PDB) holds the 3D structures of proteins, DNA, RNA and their complexes that have been determined by experiment – over 260,000 of them. Each has a four-character PDB ID, such as 1MBN. The card searches the RCSB PDB, shows an entry, and loads it into the 3D viewer.

Search

Type a PDB ID, a protein name, a UniProt accession or any words. Auto decides which it is; the menu next to it can say so instead (Protein name, UniProt accession, Full text). Narrow the results with Human only, Contains DNA and the method (X-ray, cryo-EM, NMR); sort by Best match, Newest first or Best resolution. The results come ten to a page, each with Load in 3D and Details.

An entry

Method and resolution. X-ray crystallography and cryo-EM give a resolution in ångström (Å): the smaller the number, the finer the detail. Roughly: at 1.2–1.5 Å and better single atoms can be told apart; around 2 Å the side chains are clear; beyond about 3.5 Å the path of the chain shows, and the side chains are largely modelled. NMR structures have no resolution; they are a set of models.
The paper that describes the structure.
Macromolecules: each kind of molecule, its chains (A, B …), type, length, organism, and the UniProt entry of the protein – click it to open it in the UniProt card.
Small molecules: ligands, ions and sugars, each with its three-letter code (HEM is haem, HOH water) and a drawing.
Biological assemblies: the file holds the asymmetric unit – the chains the experiment saw together – which is not always the molecule as it exists in the cell. The assembly says which it is (a monomer, a dimer …).
Files: Look inside the .pdb file shows its records with the columns coloured; Save .pdb to my files or Save .cif keeps the file in ~/structures; Send to notebook adds Python or R code that fetches the entry and its file.

The file

A structure is a text file with one line per atom. In the old .pdb format the columns are fixed:

ATOM      1  N   VAL A   1      -2.900  17.600  15.500  1.00  0.00           N
ATOM      2  CA  VAL A   1      -3.600  16.400  15.300  1.00  0.00           C

After ATOM and the atom's number come the atom name (CA is the alpha carbon), the residue (VAL), the chain (A), the residue number, the x, y and z coordinates in ångström, the occupancy and the B-factor (0 in this old file), and the element.

ATOM lines are the protein or nucleic acid, HETATM lines the ligands and water. The B-factor says how much an atom moves or how uncertain its position is (in AlphaFold models this column holds the pLDDT confidence instead). Very large structures exist only in the newer .cif (mmCIF) format. Both are read by the 3D viewer and by Biopython (Bio.PDB) in a Python notebook.

grep -c '^ATOM' /data/examples/structures/1mbn.pdb
grep '^HETATM' /data/examples/structures/1mbn.pdb | awk '{print $4}' | sort | uniq -c

AlphaFold DB and InterPro

Two cards for what a protein's sequence says about it: AlphaFold DB has a predicted 3D structure for almost every protein in UniProt; InterPro finds its families, domains and important sites.

AlphaFold DB

AlphaFold predicts a protein's structure from its sequence. AlphaFold DB (EMBL-EBI and Google DeepMind) holds the predictions for over 200 million UniProt proteins. Type a UniProt accession and Search:

The table names the model (AF-P69905-F1), the protein and its organism, and the mean pLDDT.
The picture is the predicted aligned error (PAE).
Show in the 3D viewer loads the model, coloured by confidence; Save the model (.pdb) to my files keeps it in ~/structures.

How far to trust a model

pLDDT	Colour	Meaning
above 90	dark blue	very high confidence: side chains placed well too
70 – 90	light blue	confident: the backbone is right
50 – 70	yellow	low confidence
below 50	orange	very low: often a region with no fixed structure (disordered), not a real shape

pLDDT is a confidence for each residue, about its own surroundings. It says nothing about how the domains lie relative to each other: that is what the PAE plot shows. Its colour in column x (the scored residue, along the top) and row y (the aligned residue, down the side) is the error AlphaFold expects in the position of residue x when the model is lined up on residue y: dark is a small error. Dark squares along the diagonal are domains placed confidently; light areas between them mean their arrangement is uncertain.

In the 3D viewer the model's panel has the pLDDT profile (click it to go to a residue), the PAE plot (drag a box over it: the two stretches of residues are selected, with the mean expected error between them), Hide pLDDT < 50, and Superpose onto… an experimental structure. A model is a prediction: where an experimental structure exists, compare them.

InterPro

InterPro (EMBL-EBI) brings together the protein signatures of several databases – Pfam, PROSITE, SMART, CDD and others – and groups those that describe the same thing into InterPro entries. Type a UniProt accession, choose whose signatures to show, and Search: they are drawn along the sequence and listed with their positions.

Type	Is
Family	a group of proteins that share an evolutionary origin, usually along their whole length
Domain	a distinct unit – of structure or function – that also occurs in other proteins
Homologous superfamily	proteins related by structure, even where their sequences are not similar
Repeat	a short unit that occurs several times in a protein
Site	a few important residues: an active site, a binding site, a conserved site, a modification site

The UniProt entry of a protein (its Family & Domains section) has the same information in a shorter form; InterPro shows the signatures of each database separately. Pfam's signatures are profile HMMs: to search with profiles yourself – a few families, or your own – use HMMER in the terminal (Genes and protein domains).

Genes, genomes and variants

Four cards for genes, genomes and sequencing data: Ensembl, NCBI, the UCSC Genome Browser and the European Nucleotide Archive.

Ensembl

Ensembl (EMBL-EBI) annotates the genomes of vertebrates and model organisms: genes, transcripts, proteins, variants, and how genes of different species are related. The card asks its REST API, which serves the latest release – for human, on the GRCh38 assembly.

Gene: A species and a gene symbol (CYP2C19) or Ensembl gene ID (ENSG…). The answer gives the gene's ID, type (biotype), position and strand, and a table of its transcripts – which is the canonical one, how many exons, how long the protein is – with buttons that save each transcript's cDNA, CDS or protein sequence. Save the gene sequence (FASTA); for human genes, Show in the Genome browser (GRCh38).
Sequence: The sequence of an Ensembl gene, transcript or protein ID. genomic: the whole region, introns included (with extra bases on each side if you ask); cDNA: the exons of a transcript, untranslated regions included; CDS: from the start codon to the stop codon; protein.
Variant effect (VEP): The Variant Effect Predictor: an rs ID (rs4244285), an HGVS description (ENST00000646891.2:c.1799T>A), or a chromosome, position, reference and alternative base on GRCh38 (10 94781859 G A). The answer gives the most severe consequence; the known variants at that position, with their clinical significance; allele frequencies (1000 Genomes, gnomAD); and for each transcript the consequence, its impact (HIGH, MODERATE, LOW, MODIFIER), the amino-acid change and the SIFT and PolyPhen predictions. A consequence is worked out from the position and the transcripts alone: rs4244285 comes out as synonymous (Pro227Pro), yet it is CYP2C19*2, the commonest loss-of-function allele of the gene – the base change makes a new splice site, and the mRNA is spliced wrongly. Read the clinical significance and the papers too.
Orthologues: The genes in another species – or in all of them – that descend from one gene of their last common ancestor: orthologues are separated by a speciation event (paralogues by a duplication within a genome).
Genes in a region: Every gene between two positions of a chromosome (at most 5 Mb).

GRCh38 and GRCh37

The same base has different positions on the two human assemblies: rs4244285 is at chr10:94,781,859 on GRCh38 and chr10:96,541,616 on GRCh37 (hg19). The example data of the workbench, and many older papers, use GRCh37. In code, Ensembl's GRCh37 server is https://grch37.rest.ensembl.org; the UCSC card has both.

NCBI

NCBI's E-utilities search the Entrez databases. Search: choose the database – Gene, Nucleotide, Protein, PubMed, dbSNP, ClinVar, Taxonomy or Assembly – and type search terms. Entrez's field tags narrow a search:

Search	In	Finds
`CYP2C19[sym] AND human[orgn]`	Gene	the human gene with that symbol
`myoglobin[title] AND refseq[filter]`	Protein	RefSeq myoglobin sequences
`BRCA2[gene] AND pathogenic[clinsig]`	ClinVar	pathogenic variants of BRCA2
`clopidogrel CYP2C19 AND 2020:2026[dp]`	PubMed	papers published 2020–2026

The table shows the records found, at most 50; Nucleotide and Protein records each have a FASTA button. Download sequences fetches sequences by accession, as FASTA, GenBank (nucleotide) or GenPept (protein) – the GenBank format has the annotation too: genes, coding sequences, the source organism, the paper.

Accessions say what a record is: NM_ a RefSeq mRNA, NP_ a RefSeq protein, NC_ a RefSeq chromosome or genome, J01749 a GenBank submission. The number after the dot is the version: NM_000769.4 is the fourth version of that record.

In a Python notebook, Biopython's Bio.Entrez asks the same E-utilities (give it your e-mail address: Entrez.email = "…").

UCSC Genome Browser

The DNA of any region of a genome assembly: human hg38 (GRCh38) and hg19 (GRCh37), mouse, zebrafish, fly, worm and yeast. Write the region as chr10:96541580-96541650, counting from 1 as genome browsers do; at most 1 Mb. Bases in lower case are repeats (the sequence is soft-masked); upper and lower case are the same bases. For hg38 and hg19, Show this region in the Genome browser.

The UCSC API itself counts from 0 (start=96541579 for position 96,541,580); the code that Send to notebook writes takes care of it.

European Nucleotide Archive

ENA (EMBL-EBI) – with its partners NCBI's SRA and DDBJ – keeps the raw data of published sequencing experiments. Sequencing runs takes a study (PRJEB37886), sample or run (SRR098401) accession, or a query in ENA's language (tax_eq(9606) AND library_strategy="AMPLICON"), and lists the runs: organism, instrument, library strategy and layout, number of reads, and the FASTQ files. Click a file to download it into ~/reads.

Mind the size

A human genome or exome run is many gigabytes (the two FASTQ files of SRR098401 are 7.5 and 7.7 GB): far more than fits in a browser. Small runs – of viruses, bacteria, amplicons – are a few megabytes to a few hundred. The recipe Databases from code shows how to take just the first reads of a large file.

Sequence by accession gives any ENA/GenBank sequence as FASTA (MN908947.3 is the first SARS-CoV-2 genome).

BLAST

BLAST finds the sequences in a database that are similar to yours – the first thing to do with an unknown sequence. The BLAST card sends the search to EMBL-EBI, where it runs for seconds to a few minutes, and shows the hits.

A search

Program – which depends on what your sequence is and what you search:

Program Your sequence Searches

blastp protein proteins

blastn DNA DNA

blastx DNA, translated in all six reading frames proteins

tblastn protein DNA, translated in all six frames

Protein searches find relatives much further back in evolution than DNA searches: when the sequence codes for a protein, compare proteins (blastp, or blastx from DNA).
Database. For proteins: UniProtKB/Swiss-Prot (reviewed – a good first choice), the PDB (proteins with a known structure), AlphaFold DB, human UniProtKB, or all of UniProtKB (slow). For DNA: the coding sequences in ENA (human, mammals, mouse, all), or virus sequences.
Sequence: FASTA, or just the letters.
E-value cut-off and the number of Hits to keep.
Your e-mail address: EBI asks for one with each job, and writes to it only if something goes wrong. Remember it in this browser keeps it here, on this computer only.
Run. The job's state is shown while it waits and runs.

The examples (Try:) fill in the form – a myoglobin protein, a piece of human DNA; press Run after typing your e-mail address.

Reading the hits

Column	Means
E-value	how many hits this good would turn up by chance in a database this size. The smaller, the better: below 1e-10 the similarity is almost certainly real (the sequences share an ancestor); around 0.01–1 it may be chance.
Bits	the bit score: the quality of the alignment, independent of the database's size. Higher is better.
Identity %	the share of identical residues in the aligned part
Aligned	the length of the aligned part – a hit can cover only a piece of your sequence (one domain, say)

The full BLAST report (under the table) has the alignments themselves. Save the full report keeps it; Save the table (.tsv) keeps the hits; Open the results at EBI shows them on EBI's site, where they stay for 7 days. Click a hit to open its entry (UniProt, PDB or ENA).

From code

Send to notebook → Python writes the whole job – submit, wait, fetch the hits into a pandas table. Put your e-mail address where it says YOUR_EMAIL before you run it: addresses are not written into notebooks, which you may share. R in the browser cannot submit a BLAST job (it needs a POST request); use a Python notebook.

Without BLAST

For two sequences you already have, align them in the terminal: needleman_wunsch (whole sequences), smith_waterman (the best matching parts), lastz (long DNA) – see Alignments and trees. There is no local BLAST database in the workbench.

Interactions and networks

Which proteins work together? STRING collects known and predicted associations between proteins and tests lists of genes for shared functions; IntAct holds the interactions that experiments have shown, curated from papers; NDEx shares whole networks – pathways, interactomes, the networks of published papers. Each answer can be opened in the Networks tab, to style with your own data and analyse.

STRING

STRING scores how likely two proteins are to work together, from seven kinds of evidence (channels): experiments, curated databases, genes switched on together (co-expression), genes that sit together in genomes, are fused, or are found in the same species, and papers that mention them together (text mining). The combined score, from 0 to 1, puts them together.

Interaction network: type gene or protein names (or UniProt or Ensembl IDs), one per line, choose the species, and Search.

Minimum confidence: high (0.7) keeps the associations STRING is most sure of; medium (0.4) is STRING's default.
Add neighbours: adds the proteins most strongly connected to yours, to show the neighbourhood.
The picture is STRING's own drawing; the table has a row for each pair, with the combined score and the scores of its main channels. Open in Networks makes the same network in the Networks tab, with the edges as wide as their scores.

An association is not always a physical contact: two enzymes of one pathway are associated without ever touching. Look at where a score comes from before you call it an interaction – or, in the Networks tab's From STRING…, ask for physical interactions only.

Functional enrichment: type a list of genes – those an experiment changed, say. STRING says which GO terms, pathways (KEGG, Reactome, WikiPathways) and other annotations are more common among them than among all the genes of the species, with false discovery rates (FDR). Show keeps one category. Their network in Networks opens the network of the list.

IntAct

IntAct is a database of molecular interactions, curated by hand from published experiments: which molecules bind, how that was shown (yeast two-hybrid, co-immunoprecipitation, a crystal structure …), and in which paper. Type a UniProt accession (P04637) or a protein name; several, separated by spaces, find the interactions of each.

IntAct stores each experiment as a separate piece of evidence; the table joins them, one row for each pair of molecules, with the number of pieces of evidence, the methods and the PubMed IDs of the papers.
The MI score (0 to 1) is IntAct's confidence in an interaction: higher when more experiments, of more kinds and from more papers, show it. Above 0.6 is usually called high confidence, above 0.45 medium.
Open in Networks: the network of the molecules you asked for and their partners, with an edge for each pair, as wide as its MI score. Isoforms of a protein (P04637-2) are counted as the protein itself.

NDEx

The Network Data Exchange is a public library of networks: curated pathways (NCI-PID, SIGNOR, WikiPathways), interactomes, and networks that authors share with their papers. Search it by words – a pathway, a gene, a disease, an author. The table lists each network's size and its Rights: each network has its own licence, so read it before you reuse one. Open in Networks opens a network (of up to 50,000 edges) with its layout, and all the columns of its nodes and edges.

Which one?

	STRING	IntAct	NDEx
What it holds	Associations, physical or functional, scored from many kinds of evidence – including predictions and text mining	Physical interactions shown by experiments, each with its paper	Whole networks, made by others
Coverage	Almost every protein of thousands of species	Fewer proteins, studied ones first	Whatever has been shared
Use it for	A first look at a list of genes; their shared functions	Evidence you can cite for a single interaction	A published pathway or network to compare with

Well-studied proteins have more known partners because they have been studied more: a hub in a network made from these databases can be a popular protein rather than a central one.

From code

Send to notebook gives the code of each query. In Python, a STRING network becomes a networkx graph:

import requests, pandas as pd, networkx as nx

params = {"identifiers": "\r".join(["TP53", "MDM2", "CDKN1A", "ATM", "CHEK2"]), "species": 9606,
          "required_score": 400, "caller_identity": "bioinformatics_workbench"}
edges = pd.DataFrame(requests.post("https://string-db.org/api/json/network", data=params).json())
G = nx.from_pandas_edgelist(edges, "preferredName_A", "preferredName_B", edge_attr="score")
sorted(G.degree, key=lambda x: -x[1])

Public services

Only what you type is sent – gene names, accessions, search words. Never type or paste personal or patient data. STRING's and IntAct's data are free to reuse under CC BY 4.0: cite them (below) when you publish a network or a figure made from them.

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

STRING: Szklarczyk D, Nastou K, Koutrouli M, Kirsch R, Mehryary F, Hachilif R, et al. (2025) The STRING database in 2025: protein networks with directionality of regulation. Nucleic Acids Research 53(D1): D730–D737. doi:10.1093/nar/gkae1113
IntAct: del Toro N, Shrivastava A, Ragueneau E, Meldal B, Combe C, Barrera E, et al. (2022) The IntAct database: efficient access to fine-grained molecular interaction data. Nucleic Acids Research 50(D1): D648–D653. doi:10.1093/nar/gkab1006
NDEx: Pratt D, Chen J, Pillich R, Rynkov V, Gary A, Demchak B, et al. (2017) NDEx 2.0: A Clearinghouse for Research on Cancer Pathways. Cancer Research 77(21): e58–e61. doi:10.1158/0008-5472.CAN-17-0606

Pathways, GO and enrichment

What does a protein do? Reactome names the pathways it takes part in, the Gene Ontology card lists the standard terms for its function, the processes it is part of and where it is in the cell, and g:Profiler finds what a whole list of genes has in common. Which proteins it works with – STRING, IntAct and NDEx – is on the page Interactions and networks.

Reactome

Reactome is a curated database of human biological pathways – reactions, the molecules in them, and how they connect – with predictions for other species. Type a UniProt accession: the answer lists the pathways the protein is in, with Diagram for those that have one; Disease marks pathways of a disease.

Gene Ontology

The Gene Ontology (GO) is a dictionary of terms that describe what genes do, so that all databases say it the same way. It has three parts (aspects):

Aspect	Describes	Example
Molecular function	what the protein does itself	GO:0003700 DNA-binding transcription factor activity
Biological process	the larger process it takes part in	GO:0006915 apoptotic process
Cellular component	where in the cell it is	GO:0005634 nucleus

GO terms of a protein (from QuickGO): type a UniProt accession, and choose an aspect or all three. Each term has an evidence code that says how it was assigned: EXP, IDA, IMP, IPI … from experiments; IBA, ISS from related genes; IEA automatically, without a curator; TAS, NAS from statements in papers. The qualifier says how the protein relates to the term (enables, involved in, located in, part of, or NOT).

Look up a term (from the Ontology Lookup Service) searches the terms of GO and other ontologies: the Human Phenotype Ontology (symptoms and signs), Mondo (diseases), Uberon (anatomy), the Cell Ontology, ChEBI (chemicals), the Sequence Ontology (the names of variant consequences used by VEP).

g:Profiler

Lists of genes – the genes changed in an experiment – are compared against GO and the pathway databases with an enrichment (over-representation) analysis: which terms occur more often in the list than among all the genes of the species? g:Profiler's g:GOSt answers for GO, KEGG, Reactome and WikiPathways at once. Type the genes (names or Ensembl IDs), choose the organism and the sources, and Search.

The p-values are adjusted for testing thousands of terms with g:Profiler's own method, g:SCS, made for the way GO terms overlap; terms below 0.05 are listed.
Genes is the number of your genes with the term, over the number of genes in the genome with it: a term with 3 of your genes out of 20 in the genome says more than 3 out of 2,000.
g:Profiler says which of your names it did not know – check their spelling, or use Ensembl IDs.

The test compares your list with the whole genome. If your genes could only come from a smaller set – the genes detected in your experiment – the whole genome makes many terms look enriched; on g:Profiler's own site you can give that set as the background. STRING's enrichment (on Interactions and networks) and the Networks tab's Analyse → Functional enrichment do the same test, and can use the nodes of a network as the background.

In a notebook, the recipe Databases from code asks UniProt about many proteins at once and reads the answer as a table: add go_p, go_f or go_c (biological process, molecular function, cellular component) to the fields of its query, and the table has their GO terms.

References

The papers to cite when you use this program in your work – with the version you used (Programs and versions):

g:Profiler: Kolberg L, Raudvere U, Kuzmin I, Adler P, Vilo J, Peterson H (2023) g:Profiler—interoperable web service for functional enrichment analysis and gene identifier mapping (2023 update). Nucleic Acids Research 51(W1): W207–W212. doi:10.1093/nar/gkad347

Papers and small molecules

Europe PMC searches the life-science literature; PubChem describes small molecules – drugs, metabolites, ligands – and gives their 3D models.

Europe PMC

Europe PMC holds the abstracts of PubMed and more – preprints, patents, guidelines – and millions of full texts. Type search terms; tick Only open access to keep the papers whose full text is open – free to read, with a licence that allows reuse (many more papers are free to read without such a licence). The table lists the most relevant first, with the journal, year, how often each has been cited and whether it is open access; click a title to open the paper.

Its query syntax narrows a search:

Search	Finds
`CYP2C19 clopidogrel`	papers that mention both
`"variant effect predictor"`	the exact phrase
`AlphaFold AND PUB_TYPE:review`	reviews only
`TITLE:"long-read sequencing" AND PUB_YEAR:[2024 TO 2026]`	the phrase in the title, in recent years
`AUTH:"Jumper J" AND AlphaFold`	papers by an author (initials after the surname) on a topic – an author's name alone also finds others with that name

A search finds papers; reading them is still your job. Check what kind of paper it is (a trial, a review, a case report), how large the study was, and whether others have since confirmed it.

PubChem

PubChem (NCBI) is the largest open collection of chemical information. Type the name of a compound or its PubChem CID (compound ID): the answer has its 2D drawing and its properties – formula, molecular weight, IUPAC name, SMILES, XLogP (how fat-soluble it is), and how many hydrogen-bond donors and acceptors it has.

Show in the 3D viewer loads PubChem's 3D model (a computed shape, as an .sdf file).
Save the 3D model (.sdf) keeps it in ~/structures.

The ligands in a PDB structure are named by three-letter codes (HEM); its Protein Data Bank entry lists them, with their names.

The 3D viewer

The 3D viewer shows the structures of proteins, nucleic acids and small molecules – from the Protein Data Bank, from AlphaFold DB or from your files. Style and colour them, find binding sites, measure distances, and save pictures for your reports.

Load a structure

PDB ID (top left): type 1MBN and Load – fetched from the RCSB PDB.
AlphaFold: a UniProt accession and Load model – the predicted structure, coloured by its confidence (see AlphaFold DB).
From your files: click a .pdb, .cif or .sdf file in the Files tab, or type open FILE.pdb in the terminal. Open file takes one from your computer.
From the Databases tab: Load in 3D in the Protein Data Bank card, Show in the 3D viewer in AlphaFold DB and PubChem, 3D buttons in a UniProt entry.

The workspace on the left lists what is loaded – each structure with its chains and ligands. Click a name to select it, double-click to centre on it, and use the eye to show or hide it. Several structures can be loaded together.

Move it

Rotate	drag with the left mouse button (one finger on a touch screen)
Move sideways	drag with the right button, or `Ctrl` + left button (two fingers)
Zoom	the scroll wheel or a pinch, or `Shift` + left button
Select	click an atom; `Shift`-click adds to the selection; `Esc` clears it
Centre on a residue	double-click it

The buttons at the top right: reset the view, centre on the selection, spin, change the background, full screen, and these controls ().

Styles and colours

The tabs above the picture – Display, Colour, Select, Labels & key, Measure, Views, Save image – hold the tools. Styles and colours apply to the selection, or to everything when nothing is selected: the bar above the picture always says which.

Style	Shows
Cartoon	the fold of the backbone: helices as spirals, strands as arrows, loops as tubes
Tube	a smooth line through the backbone
Sticks, Ball & stick	every bond – for ligands and side chains
Spheres	atoms at their real size: the space a molecule fills
Surface, See-through	the molecular surface – shape, pockets and grooves; see-through shows what is inside
Lines	thin bonds

Colour schemes: by chain, by molecule, rainbow (blue at the N-terminus to red at the C-terminus), secondary structure, hydrophobicity (orange hydrophobic, blue hydrophilic), residue type (acidic red, basic blue, polar green, hydrophobic sandy), element, pLDDT (AlphaFold models), B-factor (experimental structures), or one colour of your choice.

Select

Click atoms, names in the workspace, or letters in the sequence bar under the picture (it shows the sequence of one chain; choose which). The Select tab also has Box select (drag a rectangle), Select nearby (residues within some ångström of the selection), and a box for a selection query such as A:64 (residue 64 of chain A), A:60-70, HEM.

Binding sites and contacts

Select a ligand – click HEM in the workspace – and use Binding site: the residues around it are shown as sticks, with lines for the contacts they make: hydrogen bonds, salt bridges, π-stacking, cation–π, halogen bonds and metal coordination. Hydrophobic contacts are many, and are left out: the command contacts HEM hydrophobic adds them.

Measure, label, compare

Measure: click two atoms for the distance between them, in ångström.
Labels & key: labels on residues, a title on the picture, and a colour key (legend).
Views: Save view remembers the picture – orientation, styles, colours and labels – like a slide; click it to go back to it.
Superpose: lay one structure on another (an AlphaFold model on an experimental structure); the RMSD says how far apart they are. In the model's panel: Superpose onto…, or the command superpose.

The command line

Under the picture, Command line opens a command line: load 1MBN, color rainbow A, show sticks HEM; type help for the list. Most buttons write the command they stand for there too – a record of how you made a picture. All the commands are in the command reference.

Save your work

Save image: a PNG picture at 1×–4× the screen size, on the current or a transparent background, with or without the labels and key. Download PNG saves it to your computer; PNG to my files keeps it in ~/figures, for a notebook or a report; Copy image puts it on the clipboard for Word or PowerPoint.
Coordinates: a structure as it is shown – superposed, if it was – as a .pdb file, downloaded or kept in ~/structures.
Views → Save session downloads everything – structures, styles, labels, views – as one file; Open session continues from it.

The structures, styles and views are not kept when the page is reloaded: save a session to continue later.

References

The papers to cite when you use this program in your work – with the version you used (Programs and versions):

NGL Viewer (the 3D viewer): Rose AS, Bradley AR, Valasatava Y, Duarte JM, Prlić A, Rose PW (2018) NGL viewer: web-based molecular graphics for large complexes. Bioinformatics 34(21): 3755–3758. doi:10.1093/bioinformatics/bty419; Rose AS, Hildebrand PW (2015) NGL Viewer: a web application for molecular visualization. Nucleic Acids Research 43(W1): W576–W579. doi:10.1093/nar/gkv402

Viewer commands

The command line under the 3D picture takes these commands. Several can go on one line, separated by ; – a script that rebuilds a picture exactly. ↑ brings back earlier commands.

What a command acts on

Many commands end with what they act on; without it they act on the selection, or on everything when nothing is selected.

Write	For
`A` · `chain A`	chain A
`A:64` · `A:60-70` · `A:His64`	a residue, a range of residues
`HEM`	a ligand, by its name in the workspace
`ligands` · `water` · `ions` · `sugars` · `polymer` · `all`	kinds of molecule
`A or HEM`	both
`sele`	the selection
`3FXI A`	chain A of one of several loaded structures

The commands

Command	Does
`load 1MBN`	load an entry of the PDB (`load 1MBN again`: a second copy)
`load af P69905`	load the AlphaFold model of a UniProt accession
`load https://…/file.pdb`	load a structure from a web address (files of yours: `open FILE` in the terminal)
`remove 1MBN` · `clear`	remove one structure · everything, labels too
`show cartoon A`	add a style: `cartoon`, `tube`, `sticks`, `ballstick`, `spheres`, `surface`, `glass` (see-through), `lines`
`hide sticks A` · `hide all water`	remove one style, or all of them
`only sticks HEM`	just this one style
`undisplay B` · `display B` · `isolate A` · `showall`	hide and show parts, like the eyes in the workspace
`color rainbow A`	a colour scheme: `chain`, `molecule`, `rainbow`, `ss`, `hydrophobicity`, `restype`, `element`, `cpk`, `plddt`, `bfactor`, `default`
`color orange HEM` · `color #2a9d8f A:60-70`	one colour, by name or hex code
`select A:64` · `select add HEM` · `select invert` · `select none`	the selection
`select near HEM 4.5`	the residues within 4.5 Å
`site HEM 4.5`	the binding site: nearby residues as sticks, with their contacts
`contacts HEM` · `contacts HEM hydrophobic` · `contacts off`	lines for hydrogen bonds, salt bridges and other contacts
`center HEM` · `reset` · `spin on`	the camera
`background black` · `opacity 0.4`	the background; how see-through `glass` is
`label A:64` · `label3d "Distal His" A:64` · `label2d "Myoglobin" 0.05 0.05` · `labels clear`	labels on residues (`label A:64` selects A:64 first – the selection is replaced); text that moves with the molecule; a title fixed on the picture (x, y from 0 to 1)
`key on` · `key auto` · `key add red "Haem"` · `key title "Colours"` · `key off`	the colour key
`distance A:64.NE2 A:93.NE2` · `distance clear`	the distance between two atoms (chain:residue.atom)
`superpose AF-P02185-F1:A onto 1MBN:A`	lay one chain on another; reports the RMSD
`view save "Pocket"` · `view 1` · `view Pocket`	save a view; go back to it
`export png 2` (or `2x`)	download the picture at 2× the screen size
`camera` · `orient …`	print the camera's orientation as an `orient` command, to put it back exactly
`wait 1000` · `echo text` · `tab colour`	for scripts: pause (ms), show a message, open a tab of the viewer
`help`	the list of commands

Scripts

Each button runs a script in the viewer – copy the commands into the command line to change them:

load 1MBN; hide all; show cartoon; color ss; only sticks HEM; color element HEM; site HEM 4.5; center HEM

load 1MBN; color rainbow; show sticks A:64 or A:93; label A:64 or A:93; distance A:64.NE2 A:93.NE2

The AlphaFold model of sperm whale myoglobin, laid on the experimental structure:

load 1MBN; load af P02185; superpose AF-P02185-F1:A onto 1MBN:A

The genome browser

The genome browser (igv.js, the web edition of IGV) draws data along a genome: mapped reads, variants, genes, coverage. It is where you check, with your own eyes, what a program has claimed – a variant call, a deletion, a badly mapped region.

The genome

Choose it at the top left:

Example: CYP2C19 / CYP2C9 slices – the reference of the example reads: two pieces of human chromosome 10 (GRCh37), called human_CYP2C19 and human_CYP2C9, with their genes. Use it for anything made from the example data.
Hosted by the IGV team – whole genomes, with their genes: human hg38 and hg19, mouse, zebrafish, fly, worm, yeast. They need the internet.
A FASTA file of yours – any FASTA file in your home folder, once it has an index: samtools faidx genome.fa.

Another genome replaces the one shown, and the tracks you added to it go – they belong to that genome. If the new genome cannot be loaded (no internet, a server that does not answer), the one you had comes back with its tracks, and a message says why.

The chromosome names must match

A track is drawn only on a genome whose chromosome (sequence) names are the ones in the file. Reads mapped to reference.fa are on human_CYP2C19, and show only on the example genome; a VCF file that says chr10 needs hg19 or hg38 – the one it was made on.

Add your files

Add my files… lists the files among yours that the browser can show: BAM (with its .bai index – samtools index makes it), VCF (.vcf, or .vcf.gz with its .tbi), BED, GFF and GTF, bedGraph (also MACS's .bdg), the peaks of MACS3 (.narrowPeak, .broadPeak, .gappedPeak), wig, bigWig and bigBed.
In the terminal: igv mapped.bam calls.vcf.gz, igv --locus human_CYP2C19:11,566-11,666 to go somewhere, and igv --genome hg19 peaks.narrowPeak to change the genome first – a hosted one (hg38, hg19, mm39, mm10, danRer11, dm6, ce11, sacCer3) or a FASTA file of yours (see Peaks).
From the Files tab (⋯ → Show in the Genome browser; a FASTA file: Use as the genome), from a Galaxy dataset (), and from the Ensembl and UCSC cards of the Databases tab.
From my computer… shows files from your computer without copying them into your files – choose a BAM file together with its index, or a FASTA file with its .fai to use it as the genome.

Map the example reads and look at them:

mkdir -p ~/igv && cd ~/igv
minimap2 -ax sr /data/examples/reads/reference.fa /data/examples/reads/NA12878_R1.fastq /data/examples/reads/NA12878_R2.fastq 2> /dev/null | samtools sort -o mapped.bam -
samtools index mapped.bam
igv mapped.bam /data/examples/variants/NA12878.calls.vcf.gz --locus human_CYP2C19:11,566-11,666

Find your way

Type a position – human_CYP2C19:11,566-11,666, chr10:96,541,616 (on hg19) – or, on the hosted genomes, a gene name (CYP2C19) in the box at the top, and press Enter.
Drag the picture sideways to move; + and − (or a double-click) to zoom. Reads appear when you are zoomed in far enough.
Click a read, a variant or a gene for its details. The ⚙ of a track changes how it is drawn: colours, height, expanded or squished.

Reading the reads

The coverage track on top shows how many reads cover each base. Where many reads disagree with the reference, the bar is coloured by the bases they have.
Each bar is a read, coloured by its strand: pink forward, blue reverse. A coloured letter in a read is a base that differs from the reference; faint letters are bases of low quality.
A purple I is an insertion, a black line a deletion. Faint reads have mapping quality 0 (MAPQ 0): they fit equally well somewhere else in the genome.
A real heterozygous SNV shows the other base in about half the reads, on both strands; a sequencing error shows up in one or two reads, or on one strand only.

Variants

In a VCF track, each variant is a bar (its colour shows the frequency of the alternative allele); under it, each sample's genotype: dark blue heterozygous, cyan homozygous for the alternative allele, grey homozygous for the reference.

Save a picture

Save image downloads the current view as a PNG file.

References

The papers to cite when you use this program in your work – with the version you used (Programs and versions):

igv.js (the genome browser): Robinson JT, Thorvaldsdottir H, Turner D, Mesirov JP (2023) igv.js: an embeddable JavaScript implementation of the Integrative Genomics Viewer (IGV). Bioinformatics 39(1): btac830. doi:10.1093/bioinformatics/btac830

The Networks tab

Networks of interacting genes and proteins – or of anything else with links between things: open them, colour and size them with your own data, look for their hubs and clusters, and save pictures and files of them. The tab works like the desktop program Cytoscape, and reads and writes its files; it is built on Cytoscape.js.

The parts of the tab

The toolbar: the network you are looking at (click it to switch to another, rename, duplicate or close it), then Open, Data, Layout, Select, Analyse and Save, a box to find nodes, undo and redo.
The network, with the number of nodes and edges at the top left, the legend at the bottom left (− folds it away) and, at the bottom right, zoom, fit and a button that hides the panels to give the network the whole tab.
The side panel: Style (colours, sizes, shapes, labels), Select (filters), Results (of the analyses) and Details (of what is selected).
The tables underneath: one row for each node, or each edge, with all their columns. Drag the bar above them to make them taller.

Open a network

From	How
A file	Open → From my files… or From my computer…. The tab reads SIF (`.sif`), GraphML (`.graphml` – networkx, igraph, Gephi, yEd), XGMML (`.xgmml`, Cytoscape's XML), GML (`.gml`), Cytoscape JSON (`.cyjs`, also Cytoscape.js's own JSON), CX2 (`.cx2`, NDEx) – and tables of edges (`.csv`, `.tsv`): it shows the first rows and asks which columns hold the two ends of each edge.
STRING	Open → From STRING…: type genes or proteins; STRING finds them and returns their network – known and predicted associations, with a confidence score (0–1) and the evidence behind it. The edge width shows the score.
A table of genes	Open → From a table of genes…: the results of DESeq2, edgeR or limma. Choose the column of gene names or IDs and the filters (adjusted p-value, fold change): the genes that pass go to STRING, and every column of the table becomes a column of the nodes. See the recipe Expression data on an interaction network.
IntAct	Open → From IntAct…: the experimentally shown interactions of proteins (UniProt accessions or names), one edge for each pair, with IntAct's MI score and the number of experiments.
NDEx	Open → From NDEx…: search the public networks of NDEx – pathways (NCI-PID, WikiPathways, SIGNOR), interactomes, the networks of papers – and open one, with its own layout. Each network has its own licence, shown in the list.
The other tabs	Files (⋯ → Open in Networks), the terminal (`network FILE`), a Galaxy dataset of the formats `sif` or `xgmml`, and the STRING, IntAct and NDEx cards of the Databases tab (Open in Networks). From a Python notebook: `!network graph.graphml`.

network --example karate
network --string TP53 MDM2 CDKN1A ATM CHEK2
network /data/examples/networks/karate_club.graphml

What goes to the services

STRING, IntAct, NDEx and g:Profiler get the names you send them – gene names, accessions, search words – and nothing else: your tables stay in the browser. Never send personal or patient data. STRING's and IntAct's data are CC BY 4.0: cite them (below) when you publish a network made from them.

Find your way

Drag the background to move; the mouse wheel (or a pinch) zooms; F or the fit button shows the whole network.
Click a node or an edge to select it; Shift+click adds to the selection; Shift+drag draws a box around several. Drag a node to move it. A double-click zooms in on a node and its neighbours.
Right-click a node: its neighbours, its interactors from STRING, its UniProt entry, its AlphaFold structure in the 3D viewer, its gene in the Genome browser (human, hg38); right-click the background to add a node there.
Find nodes (the box in the toolbar) marks the nodes whose names contain what you type – several names, separated by spaces, at once; Enter selects them and zooms to them.
Delete deletes what is selected; Ctrl+Z undoes (deletions, layouts, imports, analyses), Ctrl+Shift+Z redoes; Esc selects nothing. (The keys work once you have clicked the network.)

The tables

The node table has a row for each node and a column for each attribute: its name, the columns of the file it came from, those of tables you imported, and those the analyses add. The edge table does the same for edges, and Network describes the network as a whole.

Click a column's name to sort by it; click a row to select its node (or edge) and go to it. only selected shows the rows of what is selected; the box filters the rows by any text.
Double-click a cell to change it. The ⋯ of a column: colour or size the nodes by it, a histogram, select by its values, rename it, delete it.
Save writes the table as a CSV file into your files.

Your data on the network

Data → Import a table of values onto the nodes… joins a table to the nodes: each row whose name (or ID) matches a node adds its columns to that node. The tab finds the pair of columns that match best – gene names with gene names, Ensembl IDs with Ensembl IDs – and says how many nodes found a row. Networks from STRING keep what you typed in the column query, so Ensembl IDs from DESeq2 match there.

Style

Every look of a node or an edge – fill colour, size, shape, label, border, opacity; edge colour, width, line, arrow – is the same for all, or mapped from a column, in the Style panel:

Continuous (a column of numbers): a range of colours or sizes. For fold changes, the blue – grey – red colours (blue down, grey unchanged, red up) with 0 in the middle; for p-values, −log₁₀ so that the smallest p-values give the largest nodes. 95% of them sets the ends of the scale so that a few extreme values do not wash the others out (those get the colour of the end, ≥ in the legend).
Discrete (any column): a colour, shape or size for each value – clusters, kinds of molecule, the two clubs of the karate club. The first eight values get eight colours chosen to be told apart by colour-blind readers too; more values are grey unless you choose a colour for them.
Label: any column, or none.
Chart: a ring of slices on each node, one for each column ticked, coloured by its value – the fold changes of several experiments on one network.

The ready-made styles set several at once: Plain, Confidence (edges wider for higher STRING or IntAct scores), Expression (fill by the fold change, size by the adjusted p-value), Clusters (a colour for each cluster) and Degree (size and colour by the number of edges). The legend follows the mappings; it is drawn in saved pictures too.

Layouts

Layout	For
Force-directed (fCoSE)	The usual first choice: nodes repel each other, edges pull them together, so connected groups come together. Separate pieces are packed side by side.
… strong edges shorter	The same, with edges of high score (STRING, IntAct) shorter.
Circle, Grid	Every node visible and in order – small networks, or to start editing.
Concentric	The most connected nodes in the middle.
Hierarchical (dagre)	Directed networks – signalling, regulation, pathways – in layers from top to bottom.
Tree	A breadth-first tree from the selected node (or the most connected one).
Groups	A circle for each value of a column: one per cluster after Analyse → Clusters.
Spring-embedded (CoSE)	Cytoscape.js's older force-directed layout.

A layout of the selected nodes only rearranges a part and leaves the rest. Layouts can be undone.

Select, edit, combine

Select: all, none, the inverse, the neighbours of what is selected, the edges between the selected nodes, nodes without edges, the shortest path between two selected nodes – and, in the Select panel, by the values of columns (padj < 0.01, |log2FoldChange| ≥ 2, a name that contains a word …), with several conditions.
What is selected becomes a new network (with the edges between its nodes), can be hidden, or deleted.
Data: add a node, an edge between two selected nodes, the interactors of the selected nodes from STRING, STRING's interactions between the nodes of a network that has none; merge two networks (union, intersection, difference – nodes matched by name or another column); treat the edges as directed (arrows).

Analyse

Analyse computes the statistics of the network (degree, betweenness, closeness, clustering coefficient …), ranks its hubs, finds clusters (MCODE, MCL, Louvain), connected components and k-cores, the shortest path between two nodes, and the functions and pathways its genes share. Each result appears in the Results panel and adds columns to the node table. What each one means, and how it was checked: Network analysis.

Save and share

Save →	What you get
Picture (PNG)	The whole network as drawn, at 1×, 2× or 4×, with the legend at the right – to download or to keep among your files.
Picture (SVG)	The same as a vector picture, for papers and posters (Inkscape, Illustrator, PowerPoint open it).
Cytoscape JSON (`.cyjs`)	The network with its positions and every column; it opens again here (with its style) and in Cytoscape (File → Import → Network from File).
GraphML, SIF	For networkx, igraph, Gephi, yEd and Cytoscape; SIF keeps only the nodes, edges and kinds of interaction.
Node table, edge table	CSV files of the tables.
To a notebook	Adds code to the open notebook that reads the network: for Python, the network is saved as GraphML in `~/networks` and read with networkx; for R, its node and edge tables are saved as CSV files there and made into a network with igraph.

The networks that are open are kept in this browser, in ~/.networks/session.json: they are there again when you come back. To keep a network for good, or take it elsewhere, save it as a file.

import networkx as nx

# a network made in Python, opened in the Networks tab
G = nx.les_miserables_graph()
nx.write_graphml(G, "lesmis.graphml")
!network lesmis.graphml

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

Cytoscape.js (the Networks tab): Franz M, Lopes CT, Huck G, Dong Y, Sumer O, Bader GD (2016) Cytoscape.js: a graph theory library for visualisation and analysis. Bioinformatics 32(2): 309–311. doi:10.1093/bioinformatics/btv557; Franz M, Lopes CT, Fong D, Kucera M, Cheung M, Siper MC, et al. (2023) Cytoscape.js 2023 update: a graph theory library for visualization and analysis. Bioinformatics 39(1): btad031. doi:10.1093/bioinformatics/btad031
fCoSE (the force-directed layout): Balci H, Dogrusoz U (2022) fCoSE: A Fast Compound Graph Layout Algorithm with Constraint Support. IEEE Transactions on Visualization and Computer Graphics 28(12): 4582–4593. doi:10.1109/TVCG.2021.3095303; Dogrusoz U, Giral E, Cetintas A, Civril A, Demir E (2009) A layout algorithm for undirected compound graphs. Information Sciences 179(7): 980–994. doi:10.1016/j.ins.2008.11.017
Cytoscape (the desktop program whose tools the Networks tab follows): Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. (2003) Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Research 13(11): 2498–2504. doi:10.1101/gr.1239303
STRING: Szklarczyk D, Nastou K, Koutrouli M, Kirsch R, Mehryary F, Hachilif R, et al. (2025) The STRING database in 2025: protein networks with directionality of regulation. Nucleic Acids Research 53(D1): D730–D737. doi:10.1093/nar/gkae1113
IntAct: del Toro N, Shrivastava A, Ragueneau E, Meldal B, Combe C, Barrera E, et al. (2022) The IntAct database: efficient access to fine-grained molecular interaction data. Nucleic Acids Research 50(D1): D648–D653. doi:10.1093/nar/gkab1006
NDEx: Pratt D, Chen J, Pillich R, Rynkov V, Gary A, Demchak B, et al. (2017) NDEx 2.0: A Clearinghouse for Research on Cancer Pathways. Cancer Research 77(21): e58–e61. doi:10.1158/0008-5472.CAN-17-0606

Network analysis

The Analyse menu of the Networks tab measures the network and its nodes, ranks its hubs, divides it into clusters and asks what its genes have in common. Everything but the enrichment is computed in your browser, in the background – Stop cancels it; the enrichment sends the names of the genes, and nothing else, to STRING or g:Profiler.

Before you start

Every analysis treats the network as undirected: an edge A → B joins A and B. Several edges between the same two nodes count as one, and self-loops are left out (the statistics count them).
Each analysis adds its results to the Results panel and columns to the node table: sort by them, colour or size the nodes by them, filter on them, save them. Ctrl+Z takes them away again.
The statistics and the hub rankings follow every shortest path between every two nodes: a second or two for 2,000 nodes, about 10 seconds for 5,000 nodes and 25,000 edges, half a minute for 10,000 nodes. Clusters take seconds.

Network statistics

Analyse → Network statistics computes what Cytoscape's NetworkAnalyzer computes for an undirected network, with the same column names, so that results can be compared with papers that used it. For each node:

Column	What it is
Degree	The number of edges at the node.
AverageShortestPathLength	The mean distance, in edges, to the nodes it can reach.
ClosenessCentrality	1 / AverageShortestPathLength: high for nodes close to the rest of their part of the network.
BetweennessCentrality	The share of the shortest paths between other nodes that pass through it, from 0 to 1 (counted within its connected component). High: a bridge between parts of the network.
Stress	The number of shortest paths that pass through it.
ClusteringCoefficient	Of the pairs of its neighbours, the share that are neighbours of each other: 1 when they all are.
NeighborhoodConnectivity	The mean degree of its neighbours.
TopologicalCoefficient	How much it shares neighbours with other nodes.
Eccentricity	The largest distance from it to a node it can reach.
Radiality	(Δ + 1 − AverageShortestPathLength) / Δ, where Δ is the diameter of its component: near 1 for central nodes.

The edge table gets EdgeBetweenness: the number of shortest paths that run along the edge. The Results panel has the network as a whole – number of nodes and edges, average number of neighbours, diameter and radius, characteristic path length (the mean distance between two nodes), clustering coefficient, density, heterogeneity, centralization, connected components – and three charts.

Degree distribution, on logarithmic axes, with a power law fitted to it (y = a·x^b, as NetworkAnalyzer fits it). Many networks have a few hubs and many nodes with few edges, and their points fall near a line here. A line on these axes is a hint, not a proof, that a network is "scale-free": a careful test needs more than a straight line, and finds that few real networks pass it (Broido & Clauset 2019).
Shortest path lengths: how many pairs of nodes are 1, 2, 3 … edges apart.
Clustering coefficient by degree: falling with the degree when hubs join groups that are not joined to each other.

Hubs

Analyse → Hubs ranks the nodes by the methods of cytoHubba, a Cytoscape app for finding the important nodes of a network. Tick the methods; the Results panel lists the top ones of each (equal scores, equal rank), and every score goes into the node table as hub_MCC, hub_Degree …

Method	A node scores high when …
MCC (maximal clique centrality)	it is in many cliques (groups in which every node is joined to every other), and large ones: the sum of (\|C\| − 1)! over its maximal cliques C. A node whose neighbours are not joined to each other scores its degree.
DMNC	the largest group of its neighbours that are joined to each other is dense: E / N^1.7 for its N nodes and E edges.
MNC	that group is large: N.
Degree	it has many neighbours.
EPC (edge percolated component)	it stays connected to many nodes when edges are removed at random: 1,000 copies of the network, each edge kept with probability ½; the sum, over the copies, of the nodes it is connected to, divided by the number of nodes.
BottleNeck	many shortest-path trees – one from each node – have more than a quarter of their nodes under it.
EcCentricity	the farthest node from it is near.
Closeness	it is near all the others: the sum of 1/distance.
Radiality	it is near all the others, compared with the diameter.
Betweenness	many shortest paths between other nodes pass through it.
Stress	the same, each path counted.
ClusteringCoefficient	its neighbours are neighbours of each other.

MCC, DMNC, MNC, Degree and the clustering coefficient look at a node and its neighbours; the others at the whole network. In the paper that introduced it, MCC found the essential proteins of yeast best; degree is simpler and not far behind; betweenness and BottleNeck find the nodes that link groups. A node that several methods put near the top is a better candidate than one that only one method likes. The buttons under the list select the top nodes, make a network of them and their neighbours (cytoHubba's "expanded subnetwork"), or colour them by rank.

Clusters and communities

Analyse →	Finds	Every node in a cluster?	Column
Clusters: MCODE	Dense regions – in protein networks often complexes	No: most nodes are in none	`MCODE_Cluster`, `MCODE_Score`
Clusters: MCL	Groups that random walks stay inside	Yes, in exactly one	`MCL_Cluster`
Communities: Louvain	A division with many edges inside groups and few between them	Yes, in exactly one	`Louvain_Community`
Connected components	Parts that no path joins	Yes	`Component`
k-cores	The core number of each node: the largest k for which it belongs to a k-core – a part of the network in which every node has at least k neighbours within the part	–	`CoreNumber`

The clusters are numbered from the best (MCODE) or the largest (the others); the nodes get the colour of their cluster, and each cluster in the Results panel can be selected, made into a new network, or tested for enrichment. Lay out by cluster puts each cluster in a circle of its own.

MCODE

MCODE (Bader & Hogue 2003) first scores every node: it finds the most densely connected core of the node's neighbourhood (the highest k-core) and multiplies its density by k. It then grows a cluster from the best-scoring node, taking in the neighbours whose scores are within 20% of the seed's (Node score cutoff 0.2), and so on outwards; then from the best node not yet in a cluster, and so on. A cluster is kept when it has a 2-core (k-core); the haircut removes nodes joined to it by a single edge; fluff (off) would add dense neighbours. The score of a cluster is its density times its number of nodes. The settings are those of the Cytoscape app clusterMaker2, whose implementation the tab follows step by step.

MCL

The Markov cluster algorithm (van Dongen; Enright et al. 2002) follows random walks on the network: they stay longer within groups that are well joined. In rounds, it spreads the walks one step further (expansion) and strengthens the strong paths and weakens the weak ones (inflation: each value to the power r, here 2), until nothing changes; each group of nodes that the walks no longer leave is a cluster. Larger inflation, more and smaller clusters (1.4 to 6 are usual). With Edge weights – STRING's score, say – strong edges hold clusters together. It runs as the Python package markov_clustering runs it, with a self-loop on each node as strong as its strongest edge (as Cytoscape's clusterMaker2 adds them); values below 0.001 are dropped, and it stops after 100 rounds (the panel says if it stopped before it converged).

Louvain

The Louvain method (Blondel et al. 2008) looks for the division with the highest modularity: the share of edges inside groups, minus the share expected if the same nodes, with the same numbers of edges, were joined at random (Newman 2006). It moves single nodes to the neighbouring group that raises it most, then treats each group as one node and starts again, until it stops rising. A modularity above about 0.3 means a clear division. Resolution above 1 gives more, smaller communities; below 1, fewer and larger. It runs as networkx's louvain_communities does, but takes the nodes in the order of the table instead of shuffling them: the same network always gives the same communities here (networkx's can differ a little from run to run).

Paths

Select two nodes, then Analyse → Shortest path between the two selected nodes (also in the Select menu): the path with the fewest edges is selected, and listed in the Results panel. When there are several equally short paths, it is one of them. Nodes in different components have no path between them.

Functional enrichment

Analyse → Functional enrichment asks which GO terms (processes, functions, cellular components), pathways (KEGG, Reactome, WikiPathways) and other annotations are more common among a set of genes than among all the genes of the species – an over-representation test, corrected for testing thousands of terms. The genes can be all the nodes, the selected ones, or a cluster (the enrichment link of each cluster).

STRING adds its own annotations (publications, tissues, protein domains, STRING's clusters), reports false discovery rates (Benjamini–Hochberg), and can use the nodes of the network as the background – tick Background when the network is a selection (say, the genes detected in your experiment) rather than a random set.
g:Profiler reports p-values adjusted with its own method, g:SCS, made for the way GO terms overlap.

Click a term to select its genes on the network; Save the table writes the terms as CSV. Two cautions: genes that are connected in STRING tend to share annotations, because shared annotations are part of what connects them; and an enriched term says the genes are annotated with it more often than by chance, not that the term explains your experiment.

Try it on the karate club

Wayne Zachary recorded the friendships among the 34 members of a university karate club (Zachary 1977). A quarrel between the instructor (member 1, "Mr. Hi") and the club's administrator (34, "Officer") split the club in two; the column club says which side each member joined, and the example colours the nodes by it. It is small enough to check the analyses by eye:

Analyse →	What you should see
Network statistics	78 edges; diameter 5, radius 3; characteristic path length 2.408; clustering coefficient 0.571; density 0.139.
Hubs	Degree: 34 (17 friends), 1 (16), 33 (12), 3 (10), 2 (9). Betweenness: 1, 34, 33, 3, 32. MCC: 1, then 2 and 3 (equal), 4, 34.
MCODE	Three clusters: 1, 2, 3, 4, 8, 14 (score 5.6); 15, 16, 19, 21, 23, 33, 34 (3.67); 5, 6, 7, 11 (2.67).
MCL (inflation 2)	Two clusters, of 19 and 15 members: the two clubs, but for members 3 and 9.
Louvain	Four communities, modularity 0.419: two for each side of the quarrel.
k-cores	The highest core is 4: members 1, 2, 3, 4, 8, 9, 14, 31, 33 and 34.
Shortest path, 17 and 26	17 → 6 → 1 → 32 → 26: four steps.

How the analyses were checked

The analyses were written for the workbench (file assets/js/net-algo.js) and compared with independent programs on five networks – the karate club, the characters of Les Misérables, the Florentine families, a random network and the STRING network of the recipe – 215 comparisons, all equal:

networkx 3.4.2: degrees, distances, closeness, betweenness and edge betweenness, clustering coefficients, eccentricity, neighbour degrees, components, core numbers, maximal cliques and modularity;
markov_clustering (Python): MCL, with and without weights;
line-by-line Python versions of clusterMaker2's MCODE and of networkx's Louvain (without the shuffle), and of the definitions in the cytoHubba paper for its other methods.

Where the published definitions leave a choice, the tab makes these:

Betweenness and stress count each pair of nodes once (not once from each end); NetworkAnalyzer's BetweennessCentrality is divided by (N − 1)(N − 2)/2, for N the nodes of the node's component.
EPC draws its 1,000 random networks from a fixed seed, so that the same network always gives the same scores.
BottleNeck takes the tree of a breadth-first search from each node, each node under the node that reached it first, and counts a node (itself included) under the nodes above it.
Radiality (cytoHubba): the sum of (Δ + 1 − distance) to the other nodes of the component, divided by their number, as radiality was first defined (Valente & Foreman 1998), times the component's share of all the nodes.

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

NetworkAnalyzer (network statistics): Assenov Y, Ramírez F, Schelhorn SE, Lengauer T, Albrecht M (2008) Computing topological parameters of biological networks. Bioinformatics 24(2): 282–284. doi:10.1093/bioinformatics/btm554; Doncheva NT, Assenov Y, Domingues FS, Albrecht M (2012) Topological analysis and interactive visualization of biological networks and protein structures. Nature Protocols 7(4): 670–685. doi:10.1038/nprot.2012.004
Shortest paths and betweenness: Brandes U (2001) A faster algorithm for betweenness centrality. The Journal of Mathematical Sociology 25(2): 163–177. doi:10.1080/0022250X.2001.9990249
cytoHubba (hubs): Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY (2014) cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Systems Biology 8(S4): S11. doi:10.1186/1752-0509-8-S4-S11
Radiality, as first defined: Valente TW, Foreman RK (1998) Integration and radiality: Measuring the extent of an individual's connectedness and reachability in a network. Social Networks 20(1): 89–105. doi:10.1016/S0378-8733(97)00007-5
MCODE (clusters): The second paper for clusterMaker2, whose implementation of MCODE the Networks tab follows step by step.; Bader GD, Hogue CWV (2003) An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics 4(1): 2. doi:10.1186/1471-2105-4-2; Utriainen M, Morris JH (2023) clusterMaker2: a major update to clusterMaker, a multi-algorithm clustering app for Cytoscape. BMC Bioinformatics 24(1): 134. doi:10.1186/s12859-023-05225-z
MCL (Markov clustering): Enright AJ, Van Dongen S, Ouzounis CA (2002) An efficient algorithm for large-scale detection of protein families. Nucleic Acids Research 30(7): 1575–1584. doi:10.1093/nar/30.7.1575; van Dongen S (2008) Graph Clustering Via a Discrete Uncoupling Process. SIAM Journal on Matrix Analysis and Applications 30(1): 121–141. doi:10.1137/040608635
Louvain communities: The second paper for modularity, which the Louvain method makes as large as it can.; Blondel VD, Guillaume JL, Lambiotte R, Lefebvre E (2008) Fast unfolding of communities in large networks. Journal of Statistical Mechanics: Theory and Experiment 2008(10): P10008. doi:10.1088/1742-5468/2008/10/P10008; Newman MEJ (2006) Modularity and community structure in networks. Proceedings of the National Academy of Sciences 103(23): 8577–8582. doi:10.1073/pnas.0601602103
NetworkX: Hagberg AA, Schult DA, Swart PJ (2008) Exploring Network Structure, Dynamics, and Function using NetworkX. Proceedings of the Python in Science Conference, pp. 11–15. doi:10.25080/TCWV9851
Degree distributions (are networks scale-free?): Broido AD, Clauset A (2019) Scale-free networks are rare. Nature Communications 10(1): 1017. doi:10.1038/s41467-019-08746-5
STRING: Szklarczyk D, Nastou K, Koutrouli M, Kirsch R, Mehryary F, Hachilif R, et al. (2025) The STRING database in 2025: protein networks with directionality of regulation. Nucleic Acids Research 53(D1): D730–D737. doi:10.1093/nar/gkae1113
g:Profiler: Kolberg L, Raudvere U, Kuzmin I, Adler P, Vilo J, Peterson H (2023) g:Profiler—interoperable web service for functional enrichment analysis and gene identifier mapping (2023 update). Nucleic Acids Research 51(W1): W207–W212. doi:10.1093/nar/gkad347
The karate club network (the example): Zachary WW (1977) An Information Flow Model for Conflict and Fission in Small Groups. Journal of Anthropological Research 33(4): 452–473. doi:10.1086/jar.33.4.3629752

The IDE tab: RStudio and Spyder

Most of your work in the workbench will be done in notebooks, but some jobs are easier in an integrated development environment, or IDE. The IDE tab has two of them: RStudio for R and Spyder for Python. Both are imitations of the real programs, made for teaching, and they work in the same way as far as possible, so what you learn here carries over to RStudio and Spyder on your own computer.

A notebook or an IDE?

A notebook is the best place to explore your data. Each piece of code sits next to its results, and you can write notes in between, so the notebook becomes a record of what you did and what you found. An IDE is built for writing programs rather than records, and it is worth switching to one when:

your code has grown beyond a few cells, and you would rather keep it in a script that you can run again from the top whenever your data change;
you are writing a report in R Markdown or Quarto, or an interactive Shiny app (both in RStudio);
your code doesn’t do what you expected, and you want to go through it one line at a time and watch your variables change (Spyder’s debugger);
you want to see all of your variables at a glance, and look inside data frames and arrays as tables;
you will be using RStudio or Spyder on your own computer later, and want to get to know them first.

The two ways of working go well together. You can try something out in a notebook and then copy the code that works into a script, or write functions in a script and load them into a notebook with source("functions.R") in R or %run functions.py in Python.

Each IDE has its own R or Python

RStudio and Spyder each run their own copy of R or Python, separate from the one that the notebooks use, just as RStudio and Spyder start R and Python processes of their own on a computer. A variable that you create in a notebook therefore doesn’t exist in Spyder, and the other way round: if you need it in both places, run the code that creates it in both. What they do share are your files. A script in Spyder can read a table that a notebook wrote, and a notebook can open a plot that RStudio saved.

The first time you open an IDE, your browser has to download R or Python, which can take a minute or so (R alone is about 30 MB). After that it keeps a copy, so the IDE opens much more quickly on your next visit. An IDE that you have opened keeps running while you use the other tabs, and you can have both open at once: the buttons at the top of the IDE tab switch between them.

Opening your files in an IDE

In the Files tab, an R script or an R Markdown or Quarto document has an Open in RStudio button above it, and a Python script has Open in Spyder. You will find the same commands in the menu of a file (its ⋯ button, or a right-click).
In the Terminal, rstudio analysis.R and spyder script.py open files in the IDE, and rstudio or spyder on its own simply opens it.
Inside an IDE, use its own File ▸ Open (Ctrl+O) or its Files pane.

As in the real programs, a file that you edit in an IDE is saved only when you save it (Ctrl+S). RStudio shows the name of a file with unsaved changes in red, and Spyder adds an asterisk to it. If you try to leave or reload the page while a file has unsaved changes, your browser asks whether you really want to.

The AI assistant and the IDEs

When code fails in RStudio or Spyder, a ✦ Ask the AI assistant about this error button appears under the error. It sends the code and the error message to the assistant, together with the script you are working on, the end of the console, and the names and types of your variables (but not their values). If you go from the IDE tab to the assistant and ask a question of your own, the same information goes with it, unless you untick your code in the IDE under the message box. Code in the assistant’s answers can then go straight into the script that is open in RStudio or Spyder, with Put into RStudio or Put into Spyder. See The AI assistant for how to connect a model.

What is different from the real programs

Everything in the IDE tab runs inside your web browser, and that brings a few limits that RStudio and Spyder on a computer don’t have:

R and Python can’t start other programs, so system() in R and subprocess in Python don’t work. Instead, a line that starts with ! in Spyder’s console runs a command in the workbench’s terminal, and RStudio has a Terminal tab of its own next to the Console.
Spyder’s debugger and input(), and stopping code while it is running in either IDE, need the page to be cross-origin isolated. The workbench arranges this itself in most browsers (see when the debugger is not available).
Reports can be knitted to web pages (HTML), but not to PDF or Word, because the programs needed for those (pandoc and LaTeX) don’t run in a browser.
There are no projects, no version control (Git) and no plug-ins or add-ins.

RStudio is made by Posit, and Spyder by the Spyder project. The IDEs in this tab are not made or endorsed by either of them.

RStudioR scripts, R Markdown and Quarto reports, Shiny apps, and RStudio’s panes
SpyderPython scripts and cells, the IPython console, the Variable Explorer and the debugger

RStudio

RStudio is the IDE that most people use for R. The one in the IDE tab is an imitation, made for teaching, that works in the same way: you write scripts and reports in an editor, run them in the Console, and see your objects, plots and help pages in the panes around them.

The four panes

Pane	What it is for
Source (top left)	the scripts and documents you are editing, each in a tab of its own, and data frames opened with `View()`
Console (bottom left)	R itself: type a command and press `Enter`. Next to it, Terminal is a shell with the same programs as the workbench’s Terminal tab, and Background Jobs shows what R prints while it knits a report.
Environment (top right)	the objects in R’s memory, and Import Dataset for reading in a data file; the History tab lists the commands you have run
Files, Plots, Packages, Help and Viewer (bottom right)	your files, the plots you have drawn, the installed packages, R’s help pages, and web pages such as knitted reports and Shiny apps

You can drag the borders between the panes to change their sizes. The menus at the top have the same commands as in RStudio, and Help ▸ Keyboard Shortcuts Help lists the keys.

Scripts

A script is a text file of R commands, saved with the extension .R. To start one, choose File ▸ New File ▸ R Script, or click New script in RStudio under this example, which reads the workbench’s table of yeast genes and draws a plot:

genes <- read.csv("/data/examples/tables/yeast_genes.csv")
str(genes)

# how many genes in each part of the cell are essential (1) or not (0)?
table(genes$location, genes$essential)

boxplot(protein_length ~ essential, data = genes, log = "y",
        names = c("not essential", "essential"),
        ylab = "Protein length (amino acids)")

Put the cursor on the first line and press Ctrl+Enter (⌘+Enter on a Mac). The line is run in the Console, and the cursor moves on to the next line, so you can keep pressing the same keys to work through the script. A command that spans several lines, such as the boxplot() above, is run as a whole, and if you select some lines first, just those lines are run. Source (Ctrl+Shift+S) runs the whole script.

Save your script with Ctrl+S; the first time, RStudio asks for a name. Because the script holds every step of your analysis, you can clear R’s memory at any time with Session ▸ Restart R, and get all your objects and plots back by running the script again.

Packages

RStudio works like R on your own computer: a package has to be installed once, with install.packages(), before library() can load it. (The notebooks are different: they install a package by themselves the first time that library() asks for it.) The packages come ready-built for webR, the version of R that runs in a browser. Most CRAN packages are there, including ggplot2 and dplyr, and so are many Bioconductor packages, such as DESeq2, edgeR and Biostrings.

install.packages("ggplot2")
library(ggplot2)

genes <- read.csv("/data/examples/tables/yeast_genes.csv")
ggplot(genes, aes(gc_content, codon_bias)) +
  geom_point(alpha = 0.3) +
  labs(x = "GC content (%)", y = "Codon bias")

Installing a package takes a little while, because the package and the other packages it needs have to be downloaded. You only have to do it once: the workbench remembers what you have installed, and fetches it again whenever library() needs it, on later visits too. The Packages pane lists what is installed, and Tools ▸ Install Packages… does the same as install.packages().

Looking at your data

The Environment pane lists your objects. Click the arrow next to a data frame to see its columns, or click the table icon at the end of its row (or type View(genes)) to open it as a table in the Source pane, where you can sort and filter it. Import Dataset, at the top of the pane, reads a CSV or Excel file through a dialog, and shows you the code that it runs, which you can copy into your script. Plots appear in the Plots pane, where the arrows take you back to earlier plots and Export saves one as a picture. To read R’s documentation, type ?mean or help(mean), or search in the Help pane.

Reports: R Markdown and Quarto

An R Markdown document (.Rmd) contains your text, written in Markdown, and your R code in chunks. You can run a single chunk with the green ▶ at its top right, and its output appears underneath it. Knit (Ctrl+Shift+K) runs the whole document from the beginning, in a new R session, and turns it into a web page with your text, code and results. The page opens in the Viewer and is saved next to the document as an .html file. Quarto documents (.qmd) are written in much the same way, and have a Render button instead of Knit.

---
title: "Essential genes in yeast"
output: html_document
---

Are the proteins of essential genes longer than the others?

```{r}
genes <- read.csv("/data/examples/tables/yeast_genes.csv")
tapply(genes$protein_length, genes$essential, median)
```

The median lengths are given for non-essential (0) and essential (1) genes.

Because knitting starts with an empty R session, the document has to create or read in everything it uses. An object that you made in the Console, but not in the document, won’t be there, and knitting stops with an error such as object 'genes' not found. That is what makes a knitted report trustworthy: it shows that the code, exactly as written, produces the results. An R Notebook (File ▸ New File ▸ R Notebook) is different. Its Preview saves a web page of the document as it is on screen, with the outputs of the chunks you have run, without running anything again.

Shiny apps

A Shiny app is an interactive web page written in R, in which sliders, menus and buttons change what R calculates and draws. File ▸ New File ▸ Shiny Web App… creates a small example app in a folder of its own, and Run App starts it in the Viewer pane. The first time, RStudio offers to install the shiny package, which takes a minute or two. While an app is running, the Console is busy, just as it would be on your computer, so press Stop (or Esc) when you have finished with it.

Shiny apps rely on the page’s service worker, a small helper that the browser runs for the workbench. If your browser doesn’t allow service workers, as some private windows don’t, RStudio tells you that the app can’t be shown.

Where R looks for files

RStudio’s R starts in your home folder, ~, so read.csv("data/counts.csv") reads the file ~/data/counts.csv. getwd() shows which folder R is in, and Session ▸ Set Working Directory changes it, for example to the folder of the script you are editing. When you knit a document, its code runs in the document’s own folder. The example data are in /data/examples, and you can read them from anywhere by giving their full paths.

Keys

`Ctrl`+`Enter`	run the current line or the selected lines
`Ctrl`+`Shift`+`S`	source the whole script
`Ctrl`+`Shift`+`Enter`	run the current chunk (in a script: source it, showing each command)
`Ctrl`+`Alt`+`I`	insert a new chunk
`Ctrl`+`Shift`+`K`	knit, render or preview the document
`Alt`+`-`	type `<-`
`Ctrl`+`Shift`+`M`	type the pipe, `\|>`
`Tab`	complete a name
`Esc` in the Console	stop the command that is running
`Ctrl`+`Shift`+`F10`	restart R

A few of RStudio’s usual shortcuts, such as Ctrl+W and Ctrl+Shift+N, belong to the browser, so they don’t work here. The menus have the same commands.

Spyder

Spyder is an IDE for Python that was designed with scientists in mind. The one in the IDE tab is an imitation, made for teaching, that looks and works like Spyder 6: you write scripts in the Editor, run them in the IPython Console, and look at your variables, plots and documentation in the panes beside them.

The panes

Pane	What it is for
Editor (left)	your scripts, each in a tab of its own, with syntax colouring, completion and warnings about likely mistakes
IPython Console (bottom right)	Python itself: type a command and press `Enter`. This is the real IPython, so its own commands, such as `%timeit` and `%pip`, work too. The History tab next to it lists the commands you have run.
Variable Explorer (top right)	the variables in the console’s memory, with their type, size and value; double-click one to look inside it. The tabs next to it are Help (documentation), Plots (your figures), Files and the Debugger.

The View menu brings any pane to the front, and View ▸ Reset window layout puts the panes back where they started. Spyder opens in its dark theme; you can switch to a light one in Tools ▸ Preferences, which also has settings for the Editor, for running files and for the console.

Running code

Try this script: New file in Spyder, under it, opens it in the Editor.

import pandas as pd
import matplotlib.pyplot as plt

# %% Read the table of yeast genes
genes = pd.read_csv("/data/examples/tables/yeast_genes.csv")
print(genes.shape)

# %% The fraction of essential genes in each part of the cell
by_place = genes.groupby("location")["essential"].mean().sort_values()
print(by_place)

# %% Plot it
by_place.plot.barh(figsize=(6, 4))
plt.xlabel("fraction of genes that are essential")
plt.tight_layout()
plt.show()

There are three ways of running a script, and you will soon find yourself using all of them:

`F5`	Run file saves the file and runs all of it, in the file’s own folder. Spyder types `%runfile` into the console, followed by the file’s name, so you can see what it is doing.
`Ctrl`+`Enter`	Run cell runs the cell that the cursor is in. A cell is a part of a script that begins with a line starting with `# %%`, and the Editor shades the one you are in. `Shift`+`Enter` runs the cell and moves on to the next one.
`F9`	Run selection runs the lines you have selected, or the line that the cursor is on, as if you had typed them into the console.

Whichever you choose, the code runs in the console’s Python, so the variables it creates are still there afterwards. You can look at them in the Variable Explorer, or carry on working with them in the console.

The IPython Console

The console is where your code runs. Type a command and press Enter. If the command isn’t finished yet, such as the first line of a for loop, you get a continuation prompt so that you can type the rest. ↑ and ↓ bring back earlier commands, and Tab completes names. A few of IPython’s own commands are worth knowing: name? shows the documentation of an object, %timeit times a line of code, %pip install NAME installs a package, and a line that begins with ! runs a command in the workbench’s terminal, as in !ls -l.

import numpy as np
x = np.random.default_rng(1).normal(size=10_000)
print(x.mean(), x.std())
%timeit np.sort(x)

NumPy, pandas and matplotlib are loaded while Spyder starts. The many other packages that come with Pyodide, such as SciPy, scikit-learn, statsmodels and Biopython (import Bio), are loaded the first time a script or a command imports them. Pure-Python packages from PyPI can be installed with %pip install NAME; they last until you reload the page.

Consoles ▸ Restart kernel (Ctrl+.) starts Python afresh, and Remove all variables clears the variables without restarting. Ctrl+C in the console stops the code that is running.

The Variable Explorer

The Variable Explorer lists the variables in the console’s Python, with their type, their size and a summary of their value, and colours the rows by type. Double-click a variable to open it: a data frame or an array as a table that you can sort and scroll through, a list or a dictionary with its contents, and a number or a piece of text in a box where you can change it. The right-click menu has more, such as plotting an array or drawing it as a histogram. Modules, functions and names that begin with an underscore are hidden; the filter button at the top of the pane changes what is shown.

Plots, help and completion

Figures that your code draws with matplotlib appear in the Plots pane, and in the console as well unless you turn that off in the pane’s menu. The arrows take you back to earlier figures, and the buttons at the top save a figure as a PNG file, copy it, zoom in or remove it.

To read the documentation of a function, put the cursor on its name, in the Editor or the console, and press Ctrl+I: the Help pane shows it. You can also type a name into the Object box at the top of that pane. As you type, Tab or Ctrl+Space offers to complete a name, and Ctrl+G (or Ctrl and a click) takes you to where a function was defined.

The Editor also checks your code as you write it, with pyflakes. A red dot in the margin marks an error, such as a bracket that is never closed, and a yellow triangle marks a warning, such as a name that is used but never defined, or a module that is imported but never used. Hover over the mark to read the message.

The debugger

When a script doesn’t do what you expect, the debugger lets you run it one line at a time and watch your variables as they change. Click just to the left of a line number, or press F12, to put a breakpoint on that line: a red dot appears. Then choose Debug file (Ctrl+F5). The script runs until it reaches the breakpoint, and stops there with the line it is about to run marked in the Editor. The Variable Explorer now shows the variables at that point, and the console’s prompt changes to IPdb [1]:, where you can type Python to look at anything you like.

`Ctrl`+`F10`	Step: run the current line, and stop at the next one
`Ctrl`+`F11`	Step into: go inside the function that the current line calls
`Ctrl`+`Shift`+`F11`	Step return: carry on until the current function returns
`Ctrl`+`F12`	Continue: run on to the next breakpoint
`Ctrl`+`Shift`+`F12`	Stop debugging

The same commands are on the toolbar, in blue, and in the Debug menu, and the Debugger pane shows how you got to the current line: which function called which. If a script has already stopped with an error, typing %debug in the console takes you to the line where it failed, with the variables as they were at that moment. To have that happen by itself whenever a script fails, tick Directly enter debugging when errors appear in Tools ▸ Preferences.

When the debugger is not available

The debugger, input() and stopping code with Ctrl+C all need Python to pause and wait for you while it is in the middle of running something. A web page can only do that when it is cross-origin isolated. The workbench arranges this itself in most browsers, but not in every situation (in some private windows, for example). If the page is not isolated, the console says so when it starts. Everything else still works, but input() stops with an error instead of waiting for an answer, the debugger can’t be started, and the only way to stop code that runs for too long is Restart kernel. Reloading the page, or opening it in an ordinary window, usually helps; What is different in the browser explains more.

Files and the working directory

The console starts in your home folder, ~. The working directory is shown at the right-hand end of the toolbar, where you can also change it, and the Files pane shows what is in it. When you run a whole file with F5, the working directory changes to the file’s folder, so a script can read the files next to it by their names alone; Run ▸ Configuration per file changes that. In the Files pane, right-click a CSV file and choose Import data to read it into a pandas data frame.

What is different from Spyder on a computer

The console’s Python runs inside your browser (Pyodide), and is separate from the Python of the notebooks.
Python can’t start threads or other processes, so threading, multiprocessing and subprocess don’t work; use !command for the programs of the terminal.
Figures are drawn in the Plots pane and the console only: there are no separate interactive windows, and desktop toolkits such as tkinter and Qt are not available.
There are no projects, no plug-ins and no version control.

The AI assistant

The AI assistant answers questions about your code, an error message, a program or a method – and puts the code from its answers into a notebook or the terminal with one click. Its Agent tab goes further: it carries out a whole task in the terminal (see The AI agent).

Before you start

Never paste personal or patient data into an AI tool, here or anywhere else.
Check what it gives you. AI models make mistakes, invent options and references, and sound sure when they do. Run the code and look at the result.
Follow your course's rules on using AI in assessed work, and say where you used it.

Connect a model

The assistant uses an AI service of your choice, with a key of your own (or one you have been given). Your questions go from your browser straight to that service. Click Not connected or (AI settings), and choose the Service:

Google Gemini – make a free key with a Google account at aistudio.google.com/apikey (you must be 18 or over). On the free tier Google may use what you send to improve its products, and people may read it – except for users in the UK, the European Economic Area and Switzerland, whose free use Google treats like paid use (it keeps what is sent for a limited time, to detect misuse, and does not use it to improve its products). Each model allows a few requests a minute. When one model is busy, the page asks another, and the answer says which one replied.
Anthropic (Claude) – a key from the Claude Console, platform.claude.com (paid).
An OpenAI-compatible service – its address (base URL) and model: a service of your university, or a model running on your own computer, if it allows requests from web pages.

Paste the key, press Test, then Save. The model's name then shows at the top of the tab.

Where the key is kept

Only in this browser tab, in memory: reloading the page forgets it. Tick Keep the key until I close this browser tab to keep it across reloads. It is never written into your files, into downloads, or into the browser's long-term storage. Forget the key removes it. A key is like a password: do not share it, and do not paste it into notebooks or files.

Ask

Type your question and press Enter (Shift+Enter for a new line); Stop ends an answer. Under the box, choose what goes with your question:

the open notebook's code – its code cells, not their output, not your data;
your file names – names and sizes, not what is in the files;
your last commands – when you come from the Terminal, your last three commands and the end of what they printed;
the open network – when you come from the Networks tab: the network's name, size and columns, its best-connected nodes, the names of the selected ones and the latest results (not the values in its tables);
your code in the IDE – when you come from the IDE tab: the script open in RStudio or Spyder (or the lines you have selected in it), the end of its console, and the names and types of its variables (not their values). A question from the IDE tab sends this instead of the open notebook’s code.

It also knows which programs the terminal has, and the versions of Python and R. From the other tabs: Ask AI on a notebook cell, or Ask the AI to explain this error under a cell that failed – this sends the cell's code with its error message (up to 2,500 characters), or with the beginning of its output (up to 1,200 characters); ✦ Ask the AI assistant about this error after a command that failed in the terminal – the command, its folder and the end of what it printed (up to 3,000 characters) – and under an error in RStudio or Spyder, which sends the code that failed (for a script that was run as a whole, the script), the error message and the IDE’s working directory, with your code in the IDE if that is ticked. The open notebook's code and your file names go with them if they are ticked under the box.

Without a key, common errors still get a short explanation from the page itself, and Copy a prompt for another AI tool copies your question, with its context, to paste into an AI tool you are allowed to use.

Code in the answers

Each block of code in an answer has buttons:

Python or R: Add to notebook, or Add and run – into the notebook that is open. Once you have opened Spyder or RStudio in the IDE tab, Put into Spyder or Put into RStudio puts the code into the script that is open there, at the cursor; when you came to the assistant from that IDE, this is the first button.
A whole R Markdown or Quarto document: Open in RStudio, as a new document that you can knit.
Shell commands: Run in the terminal (commands with rm, rmdir, mv, shred, truncate or find -delete, and > into a path that starts with ~/ or /home, ask first – other commands can overwrite files too, such as > results.txt: read them before you run them), Type in the terminal (a single command, put at the prompt without running it: check it, then press Enter), or Add to notebook as a %%bash cell.
A file with a name – a Snakefile, a script – Save as that name; Copy for any block.

⋯ (More) has Download this conversation (Markdown); New starts a new conversation.

Use it well

Say what you want and what you have: the file and its format, the program, the error message as it is.
Ask for explanations, not only answers: explain each option of this command, why did this fail, what would go wrong if …. You learn more, and you notice its mistakes.
Check against the source: man samtools or samtools view --help in the terminal, help(function) in Python, ?function in R, and these help pages.
Look at the results, not only at the code: the number of reads, the first lines of a file, a plot. Code that runs is not code that is right.

The AI agent

Give the agent a whole task – map the example reads, call the variants and tell me how many there are. It carries it out in the Terminal, one command at a time, reads what each command prints before it decides the next, and reports. Every run is kept as a notebook.

Give it a task

Connect a model first (see The AI assistant), then open the Agent tab.
Type a task – or click one under Try a task – and Send.
Watch: each step shows the command and what it printed, here and in the Terminal. Stop ends the task at any time.
Read its report – then check it (below).

The agent works in a folder of its own, ~/ai-agent/run-1, run-2 …; programs put the files they make there. A further task goes on in the same folder, and the agent is told what was done before – until New starts a new folder.

What it may do

Run the commands of the terminal: the bioinformatics programs, the text tools, python and Rscript scripts, Snakemake. Not cd – every step starts in its folder, and it uses paths instead.
Read your files. What it changes outside its folder – a file it deletes or overwrites by mistake – is put back at the end of each step. What you change yourself while it works is yours, and stays.
Use the internet (curl, wget) only if you allow it: AI settings → Let the agent use the internet. Without it, curl and wget are refused and the agent is told not to use the network from Python or R – but a Python or R script could still reach the web, so read what it runs.

Limits: at most 20 steps for a task; a command that runs for more than 150 seconds is stopped; the model sees the last 3,000 characters of what each command prints. While the agent works, the terminal is its own: wait, or press Stop.

What is sent to the AI service

The task, the commands and everything they print – including what the agent reads from your files. Do not let it work on files with personal or patient data.

The record

Each run is written into its folder as a notebook, agent-record.ipynb: the task, each command as a %%bash cell with its output, and the report. Under the agent's card: Open the record (in the Notebook tab), .ipynb (download it), Folder (show it in Files). Run the notebook's cells again to see that the commands really make what the report says. The cells run in the notebook's own folder: to keep the agent's files as they are, copy the whole folder and run the copy of the record.

Check its work

Read the commands. Are the options what you would choose? The right reference, the right files, the right genome assembly?
Look at the results yourself: open the files, count, put the reads in the genome browser. The report can be wrong even when the commands are right – and a model may describe output it has not seen.
Could you do it without the agent? If not, ask the Chat tab to explain each step, and do it yourself in a notebook.

Recipes

A recipe is a whole analysis on the example data, step by step: what to do, why, and the code that does it. Each one opens as a notebook with one button, so that you can run it cell by cell, change it and keep it.

Reads to variants, in a notebookQC, mapping, calling, a table of the variants, VEP, the genome browser
RNA-seq: reads to differential expressionHISAT2, spliced reads in the genome browser, featureCounts, DESeq2 and edgeR
Expression data on an interaction networkDESeq2's results on a STRING network: colour the genes by fold change, find the hubs and clusters, and ask what the clusters do
From a protein to its structureUniProt, Biopython, a PDB file, an AlphaFold model, superposition
A table, in Python and in RRead, summarise, plot and test the yeast gene table
Databases from coderequests, batches and loops, retries, Biopython's Entrez, part of a big file, R
Keep and share your workEvery tab as a notebook; downloads; running a notebook elsewhere; versions

Reads to variants, in a notebook

A whole analysis in one notebook: check and trim the example reads, map them, call the variants, look at them in a table and a plot, ask Ensembl what they do – and keep every step, so the notebook is the record of the work.

The button writes the steps below into a new Python notebook in ~/notebooks; run its cells from the top. Or add the cells one by one with Add to notebook.

1. A folder for the work

The results go into a folder of their own. %cd makes it the folder of the notebook's code and of its %%bash cells.

import os
os.makedirs(os.path.expanduser("~/recipe-variants"), exist_ok=True)
%cd ~/recipe-variants

2. Check and trim the reads

fastp removes adapters, drops reads of poor quality, and reports on the reads before and after (it trims low-quality ends too when asked: --cut_tail). Its report, fastp.html, opens in the Files tab.

%%bash
fastp -i /data/examples/reads/NA12878_R1.fastq -I /data/examples/reads/NA12878_R2.fastq \
      -o trimmed_1.fastq -O trimmed_2.fastq -h fastp.html -j fastp.json 2> fastp.log
grep -A 4 "after filtering" fastp.log | head -10

3. Map them to the reference

The reads come from two regions of chromosome 10; reference.fa holds those two regions. -R names the sample in the BAM file, so that the VCF file names it too.

%%bash
minimap2 -ax sr -R '@RG\tID:NA12878\tSM:NA12878' /data/examples/reads/reference.fa \
         trimmed_1.fastq trimmed_2.fastq 2> minimap2.log | samtools sort -o mapped.bam -
samtools index mapped.bam
samtools flagstat mapped.bam

4. Call the variants

%%bash
bcftools mpileup -f /data/examples/reads/reference.fa -a AD,DP mapped.bam 2> /dev/null \
  | bcftools call -mv -Oz -o calls.vcf.gz
bcftools filter -i 'QUAL>=20 && INFO/DP>=5' -Oz -o filtered.vcf.gz calls.vcf.gz
bcftools index -t filtered.vcf.gz
echo "variants called: $(bcftools view -H calls.vcf.gz | wc -l), after filtering: $(bcftools view -H filtered.vcf.gz | wc -l)"

5. The variants as a table

bcftools query writes the fields you name as a table; pandas reads it.

%%bash
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t%INFO/DP[\t%GT]\n' filtered.vcf.gz > variants.tsv
head -3 variants.tsv

import pandas as pd

v = pd.read_csv("variants.tsv", sep="\t", names=["chrom", "pos", "ref", "alt", "qual", "depth", "genotype"])
v["type"] = ["SNP" if len(r) == 1 and len(a) == 1 else "indel" for r, a in zip(v.ref, v.alt)]
v.groupby(["chrom", "type", "genotype"]).size()

6. Quality and depth

Calls with few reads have low quality: the filter of step 4 removed the weakest.

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(6, 3.5))
for gt, d in v.groupby("genotype"):
    ax.scatter(d.depth, d.qual, label=gt, s=20)
ax.set_xlabel("reads covering the position (DP)")
ax.set_ylabel("call quality (QUAL)")
ax.legend(title="genotype")
plt.tight_layout()
plt.savefig("quality_depth.png", dpi=150)
plt.show()

7. What do the variants do?

Ensembl's Variant Effect Predictor (VEP) says what each variant changes. The positions in the VCF file are on the two slices of the reference; on chromosome 10 of GRCh37 the slices start at 96,530,001 and 96,699,001, so the positions are moved before they are sent – to Ensembl's GRCh37 server.

import time
import requests

def post_json(url, data, tries=4):
    """POST data as JSON; when the service is busy or does not answer, wait and try again."""
    for attempt in range(tries):
        try:
            r = requests.post(url, json=data, headers={"Content-Type": "application/json", "Accept": "application/json"})
            if r.status_code not in (429, 500, 502, 503, 504):
                r.raise_for_status()
                return r.json()
        except requests.exceptions.ConnectionError:
            pass
        time.sleep(2 ** attempt)
    raise RuntimeError(f"no answer from {url} – try again in a few minutes")

start = {"human_CYP2C19": 96_530_000, "human_CYP2C9": 96_699_000}
lines = [f"10 {p + start[c]} . {r} {a} . . ." for c, p, r, a in zip(v.chrom, v.pos, v.ref, v.alt)]
answer = post_json("https://grch37.rest.ensembl.org/vep/homo_sapiens/region", {"variants": lines})
vep = pd.DataFrame([{"position (GRCh37)": x["input"].split()[1], "consequence": x["most_severe_consequence"],
                     "known as": ", ".join(c["id"] for c in x.get("colocated_variants", []) if c["id"].startswith("rs"))}
                    for x in answer])
vep

Position 96,541,616 is rs4244285, the CYP2C19*2 allele. VEP calls it a synonymous variant – the changed codon still codes for proline – yet it is the commonest reason for a CYP2C19 enzyme that does not work: the change creates a new splice site in exon 5, the mRNA is spliced wrongly, and no working enzyme is made. A prediction is not the whole story: read what is known about the variants that matter (Europe PMC, ClinVar in the Databases tab). NA12878 carries one copy (genotype 0/1); people with two copies turn the drug clopidogrel into its active form poorly.

8. Look at the evidence

In the workbench, igv shows the reads and the calls in the Genome browser (this line works only here, not in Jupyter):

!igv mapped.bam filtered.vcf.gz --locus human_CYP2C19:11,566-11,666

Keep it

Restart Python and run the notebook from the top: if every cell runs, the notebook is a complete record. Download → Jupyter notebook (.ipynb) keeps it; the files it made are in ~/recipe-variants. The same steps as Galaxy tools: Galaxy; as a Snakemake pipeline: Scripts and pipelines.

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

fastp: Chen S, Zhou Y, Chen Y, Gu J (2018) fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34(17): i884–i890. doi:10.1093/bioinformatics/bty560
minimap2: Li H (2018) Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34(18): 3094–3100. doi:10.1093/bioinformatics/bty191; Li H (2021) New strategies to improve minimap2 alignment accuracy. Bioinformatics 37(23): 4572–4574. doi:10.1093/bioinformatics/btab705
SAMtools: Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. (2021) Twelve years of SAMtools and BCFtools. GigaScience 10(2): giab008. doi:10.1093/gigascience/giab008; Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. (2009) The Sequence Alignment/Map format and SAMtools. Bioinformatics 25(16): 2078–2079. doi:10.1093/bioinformatics/btp352
BCFtools: Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. (2021) Twelve years of SAMtools and BCFtools. GigaScience 10(2): giab008. doi:10.1093/gigascience/giab008
Ensembl VEP: McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GRS, Thormann A, et al. (2016) The Ensembl Variant Effect Predictor. Genome Biology 17: 122. doi:10.1186/s13059-016-0974-4
igv.js (the genome browser): Robinson JT, Thorvaldsdottir H, Turner D, Mesirov JP (2023) igv.js: an embeddable JavaScript implementation of the Integrative Genomics Viewer (IGV). Bioinformatics 39(1): btac830. doi:10.1093/bioinformatics/btac830

RNA-seq: reads to differential expression

Two halves of an RNA-seq analysis. First, map real RNA-seq reads with HISAT2, see reads split across an intron, and count the reads of each gene with featureCounts. Then, from the gene counts of a published experiment, find the genes that a drug switches on or off – with DESeq2, and with edgeR for comparison.

The button writes the steps below into a new R notebook in ~/notebooks; the terminal's steps become %%bash cells. Run the cells from the top.

1. A folder for the work

setwd() makes it the folder of the notebook's R code and of its %%bash cells.

dir.create("~/recipe-rnaseq", showWarnings = FALSE)
setwd("~/recipe-rnaseq")

2. Map RNA-seq reads with HISAT2

256 pairs of reads from yeast RNA (study GSE110004), and chromosome I of the yeast genome as the reference. hisat2-build indexes the reference once; hisat2 maps the reads and reports how many it placed.

hisat2-build /data/examples/rnaseq/yeast_chrI.fa yeast_chrI > /dev/null 2>&1
hisat2 -x yeast_chrI -1 /data/examples/rnaseq/SRR6357070_1.fastq -2 /data/examples/rnaseq/SRR6357070_2.fastq \
       2> hisat2.log | samtools sort -o rnaseq.bam -
samtools index rnaseq.bam
cat hisat2.log

3. Reads across an intron

A read from a spliced RNA can start in one exon and end in the next. HISAT2 splits it: an N in the CIGAR string is the intron that the read skips. 80M366N21M means 80 bases on the genome, 366 skipped, then 21 more. How many reads skip an intron, and how long are the introns they skip?

samtools view rnaseq.bam | awk '$6 ~ /N/' | wc -l
samtools view rnaseq.bam | awk '$6 ~ /N/ {match($6, /[0-9]+N/); print substr($6, RSTART, RLENGTH)}' | sort | uniq -c | sort -rn

Most of them skip the same 366 bases: the intron of EFB1 (YAL003W), a yeast gene at I:142,174–143,160, whose intron holds a small nucleolar RNA gene (snR18). – the reads are drawn with a thin line across the intron, under the gene's exons.

4. Count the reads of each gene

featureCounts counts the read pairs that fall in each gene of the annotation. -s says which strand the reads come from: these were made with a reverse-stranded kit (-s 2). Compare the summaries – with the wrong strand, -s 1, hardly a pair finds its gene.

featureCounts -p --countReadPairs -s 2 -a /data/examples/rnaseq/yeast_chrI.gtf -o counts.txt rnaseq.bam 2> featurecounts.log
awk '$2 > 0' counts.txt.summary
tail -n +3 counts.txt | sort -k7,7nr | head -5 | cut -f 1,7
featureCounts -p --countReadPairs -s 1 -a /data/examples/rnaseq/yeast_chrI.gtf -o counts_s1.txt rnaseq.bam 2> /dev/null
awk '$2 > 0' counts_s1.txt.summary

EFB1 (YAL003W) has the most: 120 of the 239 pairs that were counted. A table like this, with a column for each sample, is where differential expression starts. With one sample of 256 read pairs there is nothing to compare – the next steps take the counts of a published experiment. (RNA-seq: counts and abundances does the same with StringTie, kallisto and Salmon.)

5. The counts of a whole experiment

Differential expression starts from a table: how many reads each gene got in each sample. For whole genomes and tens of millions of reads such tables are made on a server (featureCounts, as above, or Salmon); this one is from the Bioconductor package airway: human airway smooth muscle cells, four cell lines (cell), each untreated and treated with the steroid dexamethasone (dex) – Himes et al., PLoS ONE 2014.

counts <- as.matrix(read.delim("/data/examples/rnaseq/airway_counts.tsv", row.names = 1, check.names = FALSE))
samples <- read.delim("/data/examples/rnaseq/airway_samples.tsv", row.names = 1)
samples <- samples[colnames(counts), ]
samples$cell <- factor(samples$cell)
samples$dex <- factor(samples$dex, levels = c("untrt", "trt"))
dim(counts)
samples[, c("cell", "dex")]
colSums(counts) / 1e6   # reads per sample, in millions

6. DESeq2

The design ~ cell + dex compares treated with untreated cells of the same cell line. Genes with almost no reads are left out first. The first library(DESeq2) downloads about 40 packages – give it a minute; DESeq() takes about as long.

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = counts, colData = samples, design = ~ cell + dex)
dds <- dds[rowSums(counts(dds)) >= 10, ]
dds <- DESeq(dds)
res <- results(dds, contrast = c("dex", "trt", "untrt"), alpha = 0.05)
summary(res)

DESeq() normalises for the depth of each sample (size factors), estimates how much each gene varies between replicates (dispersion), and tests every gene. padj is the p-value adjusted for testing thousands of genes: a cut-off of 0.05 keeps the false discoveries near 5%.

7. The genes that changed most

top <- as.data.frame(res[order(res$padj), ])
head(top[, c("baseMean", "log2FoldChange", "padj")], 10)
write.csv(top, "deseq2_dex_vs_untreated.csv")

A log2 fold change of 2 is four times as much RNA in treated cells; −2 a quarter. Third in the list is DUSP1 (ENSG00000120129), a well-known target of steroids such as dexamethasone. Search an ID in the Databases tab (Ensembl card) to see the gene's name.

8. Plots

The MA plot shows every gene: its mean count and its fold change, the significant ones in colour. The PCA plot of variance-stabilised counts shows which samples resemble each other – by treatment first, then by cell line.

plotMA(res, ylim = c(-5, 5))

vsd <- vst(dds, blind = FALSE)
plotPCA(vsd, intgroup = c("dex", "cell"))

plotCounts(dds, gene = "ENSG00000120129", intgroup = c("dex", "cell"))

9. The same with edgeR

edgeR models the counts in a similar way, and tests them with quasi-likelihood F-tests. How many of DESeq2's genes does it find too? (edgeR loads limma, whose plotMA hides DESeq2's: after this step, write DESeq2::plotMA(res).)

library(edgeR)
y <- DGEList(counts = counts, samples = samples)
keep <- filterByExpr(y, group = samples$dex)
y <- normLibSizes(y[keep, , keep.lib.sizes = FALSE])
design <- model.matrix(~ cell + dex, data = samples)
y <- estimateDisp(y, design)
fit <- glmQLFit(y, design)
qlf <- glmQLFTest(fit, coef = "dextrt")
summary(decideTests(qlf))

edger_genes <- rownames(topTags(qlf, n = Inf, p.value = 0.05))
deseq2_genes <- rownames(top)[which(top$padj < 0.05)]
c(DESeq2 = length(deseq2_genes), edgeR = length(edger_genes), both = length(intersect(deseq2_genes, edger_genes)))

10. On an interaction network

Which of these genes work together? The recipe Expression data on an interaction network puts the table that step 7 wrote onto the STRING network of the genes, coloured by their fold changes, and looks for its hubs and clusters.

Good to know

DESeq2 and edgeR need raw counts – whole numbers, not normalised values – and biological replicates. In the browser they run as on a computer: locfit, a package they need, was built for the workbench, and a check against R on Linux gave the same results. Download the notebook to run it in RStudio or Jupyter: BiocManager::install(c("DESeq2", "edgeR")) there.

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

HISAT2: Kim D, Paggi JM, Park C, Bennett C, Salzberg SL (2019) Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nature Biotechnology 37(8): 907–915. doi:10.1038/s41587-019-0201-4
SAMtools: Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. (2021) Twelve years of SAMtools and BCFtools. GigaScience 10(2): giab008. doi:10.1093/gigascience/giab008; Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. (2009) The Sequence Alignment/Map format and SAMtools. Bioinformatics 25(16): 2078–2079. doi:10.1093/bioinformatics/btp352
featureCounts: Liao Y, Smyth GK, Shi W (2014) featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30(7): 923–930. doi:10.1093/bioinformatics/btt656
DESeq2: Love MI, Huber W, Anders S (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15: 550. doi:10.1186/s13059-014-0550-8
edgeR: Chen Y, Chen L, Lun ATL, Baldoni PL, Smyth GK (2025) edgeR v4: powerful differential analysis of sequencing data with expanded functionality and improved support for small counts and larger datasets. Nucleic Acids Research 53(2): gkaf018. doi:10.1093/nar/gkaf018; Robinson MD, McCarthy DJ, Smyth GK (2010) edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26(1): 139–140. doi:10.1093/bioinformatics/btp616

Expression data on an interaction network

The genes that a drug switches on or off make a list; the proteins they code for work together in a network. This recipe takes the result of the RNA-seq recipe – the genes that the steroid dexamethasone changes in airway smooth muscle cells – onto their STRING network, colours them by how much they changed, and looks for the hubs and clusters of the network and for what the clusters do. Most of it happens in the Networks tab, with the buttons below; the last steps do the same in Python.

1. The results of DESeq2

A row for each gene of the experiment: its Ensembl ID, its mean count (baseMean), how much it changed in treated cells (log2FoldChange: 1 is twice as much RNA, −1 half) and the adjusted p-value (padj). It is the file deseq2_dex_vs_untreated.csv that the RNA-seq recipe writes; a copy is in /data/examples/networks/airway_deseq2.csv.

2. From the table to a network

(if you made your own table in the RNA-seq recipe, choose it there instead). The tab finds the column of IDs and DESeq2's columns padj and log2FoldChange by their names. Keep the genes with padj < 0.01 and |log2FoldChange| ≥ 1.5 – 409 genes –, leave human, medium confidence (0.4) and all associations, and Make the network.

The IDs go to STRING, which finds the proteins – 357 of the 409 genes. Most of the others make no protein (long non-coding RNAs, pseudogenes), or have IDs that Ensembl has since retired: the airway data use an older annotation.
STRING returns the associations among the 357 proteins with a score of at least 0.4: 512 of them.
Every column of the table is added to the nodes – the rows are matched on the IDs that were sent, which the nodes keep in the column query.
The nodes are coloured by log2FoldChange – red up, blue down, grey unchanged – and sized by padj, the smallest p-values largest. The ends of the colour scale hold 95% of the genes, so that a few extreme ones do not wash out the rest (≤ and ≥ in the legend).

The finished network (the button at the top) was made in the same way on 13 October 2026, from STRING 12.5, and opens without the internet. STRING changes from version to version: made today, the network can differ a little.

3. Read the network

About a hundred proteins have no association with the others at this confidence. , , and : the rest is one large component of 237 proteins and some small ones.

Where are the genes that every steroid experiment finds – the classic targets of the glucocorticoid receptor? (or type their names into Find nodes). All eight are red: switched on. ZBTB16 changed most – a log2 fold change of 7.35, over 150 times as much RNA. Click one for its details, and its row in the node table; right-click it for its UniProt entry or its AlphaFold structure.

4. Hubs

with MCC, Degree, Betweenness and Closeness ranks the proteins. The most connected are leptin (LEP, 26 partners), thrombospondin 1 (THBS1, 21), the chemokine CXCL12 (19) and FOXO1 (16); MCC puts CXCL12, THBS1, LEP, VCAM1 and TGFB2 first; betweenness LEP, THBS1, FOXO1, CAMK2A and GPX3 – the proteins that link the parts of the network. The hubs are a mix of genes switched on (LEP, THBS1, FOXO1) and off (CXCL12, VCAM1, TGFB2).

A hub of a network made from databases is partly a well-studied gene: leptin and the chemokines are the subject of thousands of papers, and STRING's evidence – experiments, curated databases, text mining – grows with them. A hub is a reason to look closer, not a finding on its own. (What each ranking measures: Network analysis.)

5. Clusters

finds the densely connected regions: fourteen with the usual settings, each in a colour of its own. to see them apart. The eight with the highest scores (↑ switched on, ↓ off):

Cluster	Proteins	What they are
1	COL4A3↑ COL4A4↑ COL11A1↑ ITGA10↑ ITGA11↓	Collagens and the integrins that bind them
2	ATF3↑ DUSP1↑ NR4A3↑ EGR2↓ EGR3↓	Immediate-early genes: switched on within minutes of a signal
3	15 proteins, all but two ↓: BDKRB1, BDKRB2, ADRA2A, CCKAR, MCHR1, ADCY1, ADCY4, ADCY8, CAMK2A …	Receptors that act through G proteins, the adenylyl cyclases that make cAMP, and what follows
4	PDK4↑ FABP4↑ FOXO1↑ SCD↓	The use of fats and sugars
5	NEXN↑ MYOZ2↑ PKP2↑ PRKAG2↑	Muscle and cell-junction proteins
6	MT1E↑ MT1M↑ MT1X↑ MT2A↑	Metallothioneins – among the first genes found to be switched on by glucocorticoids
7	HSPA2↑ DNAJB4↑ HSPA6↓ HSPA12B↓	Heat-shock proteins, which help other proteins fold
8	FZD5↑ FZD8↑ FRZB↑ WNT2↓	Wnt signalling

The other six are smaller or looser – among them the chemokines CCL8 and CXCL12 with the adhesion molecule VCAM1, all switched off (cluster 9), and the clock genes PER1, PER3 and DBP (cluster 14). Several clusters lean one way – the metallothioneins all up, the receptors of cluster 3 nearly all down: the steroid switches groups of proteins that work together on or off. Try Clusters: MCL or Communities: Louvain too: they put every protein in a group, so their groups are larger and looser.

6. What the clusters do

Each cluster in the Results panel has a link enrichment: its genes go to STRING, which lists the GO terms and pathways they share more than chance would. : ECM–receptor interaction and focal adhesion (KEGG), integrin cell surface interactions (Reactome), collagen trimers (GO). Click a term to select its genes on the network. Cluster 3 gives signalling by G-protein-coupled receptors, cAMP and calcium – and, among the WikiPathways, Albuterol and budesonide therapy for asthma: a bronchodilator and a steroid, for the disease these airway cells are studied for. Cluster 4 gives the PPAR pathway and the response to fatty acids; Analyse → Functional enrichment of the whole network, among others, the regulation of lipolysis, the extracellular matrix, signalling by G-protein-coupled receptors, and cytokine–cytokine receptor interaction.

Five genes are a small sample: an enriched term among them is a hint. And the test compares the genes with the whole genome, although only genes switched on in these cells could have been found; with the box Background, STRING compares them with the network instead.

7. Keep it

draws the network with its legend, at up to four times the size of the screen; Picture (SVG) for a poster or a paper. Save → Network: Cytoscape JSON keeps the network with its layout, its colours and every column, to open here again or in Cytoscape. The networks that are open are also kept in this browser for your next visit.

8. The same in Python

The network can be made and analysed in a notebook too, with STRING's API and networkx – useful to repeat it on other experiments, or to keep the analysis with the code of the rest of the project. The table, filtered as in step 2:

import pandas as pd, requests, time
import networkx as nx

de = pd.read_csv("/data/examples/networks/airway_deseq2.csv", index_col=0)
sig = de[(de.padj < 0.01) & (de.log2FoldChange.abs() >= 1.5)].sort_values("padj")
len(sig)

STRING's identifiers for the genes, then the network among them. STRING asks for a second's pause between requests, and for a name for the program that asks (caller_identity):

STRING = "https://string-db.org/api/json/"
me = "bioinformatics_workbench"
found = pd.DataFrame(requests.post(STRING + "get_string_ids", data={
    "identifiers": "\r".join(sig.index), "species": 9606, "limit": 1, "echo_query": 1, "caller_identity": me}).json())
time.sleep(1)
edges = pd.DataFrame(requests.post(STRING + "network", data={
    "identifiers": "\r".join(found.stringId), "species": 9606, "required_score": 400, "caller_identity": me}).json())
len(found), len(edges)

A networkx graph, with DESeq2's values on its nodes; its most connected and most central proteins:

G = nx.Graph()
for _, g in found.iterrows():
    G.add_node(g.preferredName, ensembl=g.queryItem,
               log2FoldChange=float(de.log2FoldChange[g.queryItem]), padj=float(de.padj[g.queryItem]))
G.add_weighted_edges_from(zip(edges.preferredName_A, edges.preferredName_B, edges.score), weight="score")
print(G.number_of_nodes(), "proteins,", G.number_of_edges(), "interactions,", nx.number_connected_components(G), "components")

degree = pd.Series(dict(G.degree())).sort_values(ascending=False)
betweenness = pd.Series(nx.betweenness_centrality(G, normalized=False)).sort_values(ascending=False)
pd.DataFrame({"degree": degree.head(8), "log2FoldChange": [G.nodes[n]["log2FoldChange"] for n in degree.head(8).index]})

The same numbers as the tab's – and betweenness puts LEP, THBS1 and FOXO1 first again (betweenness.head()). Louvain communities, and the network in the Networks tab:

communities = nx.community.louvain_communities(G, seed=1)
big = sorted((c for c in communities if len(c) > 1), key=len, reverse=True)
print(len(big), "communities, modularity", round(nx.community.modularity(G, communities), 3))
nx.write_graphml(G, "airway_string.graphml")

!network airway_string.graphml

networkx takes the nodes in a random order (fixed here with seed=1), the Networks tab in the order of its table: their communities are alike, but not the same. In the tab, colour the nodes again by log2FoldChange: the Style panel, Fill colour, continuous.

Good to know

STRING's associations are not all physical: two proteins of the same pathway are associated without touching. For a network of physical contacts only, choose physical interactions only in step 2, or IntAct (Interactions and networks). The gene names go to STRING; the table stays in the browser – never send personal or patient data to these services.

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

DESeq2: Love MI, Huber W, Anders S (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15: 550. doi:10.1186/s13059-014-0550-8
The airway experiment (the example): Himes BE, Jiang X, Wagner P, Hu R, Wang Q, Klanderman B, et al. (2014) RNA-Seq Transcriptome Profiling Identifies CRISPLD2 as a Glucocorticoid Responsive Gene that Modulates Cytokine Function in Airway Smooth Muscle Cells. PLoS ONE 9(6): e99625. doi:10.1371/journal.pone.0099625
STRING: Szklarczyk D, Nastou K, Koutrouli M, Kirsch R, Mehryary F, Hachilif R, et al. (2025) The STRING database in 2025: protein networks with directionality of regulation. Nucleic Acids Research 53(D1): D730–D737. doi:10.1093/nar/gkae1113
cytoHubba (hubs): Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY (2014) cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Systems Biology 8(S4): S11. doi:10.1186/1752-0509-8-S4-S11
MCODE (clusters): The second paper for clusterMaker2, whose implementation of MCODE the Networks tab follows step by step.; Bader GD, Hogue CWV (2003) An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics 4(1): 2. doi:10.1186/1471-2105-4-2; Utriainen M, Morris JH (2023) clusterMaker2: a major update to clusterMaker, a multi-algorithm clustering app for Cytoscape. BMC Bioinformatics 24(1): 134. doi:10.1186/s12859-023-05225-z
NetworkX: Hagberg AA, Schult DA, Swart PJ (2008) Exploring Network Structure, Dynamics, and Function using NetworkX. Proceedings of the Python in Science Conference, pp. 11–15. doi:10.25080/TCWV9851

From a protein to its structure

Sperm whale myoglobin, from its UniProt entry to its structure: the sequence and its properties, the haem pocket of the experimental structure 1MBN measured with Biopython, and the AlphaFold model compared with the experiment.

1. The UniProt entry

import requests

acc = "P02185"
entry = requests.get(f"https://rest.uniprot.org/uniprotkb/{acc}.json").json()
print(entry["uniProtkbId"], "–", entry["proteinDescription"]["recommendedName"]["fullName"]["value"])
print(entry["organism"]["scientificName"], "·", entry["sequence"]["length"], "amino acids")

# the binding sites UniProt lists: which residue binds what
for f in entry["features"]:
    if f["type"] == "Binding site":
        print(f["location"]["start"]["value"], f.get("ligand", {}).get("name", ""), f.get("description", ""))

His 65 binds the oxygen, His 94 the haem's iron. UniProt numbers the sequence from its first methionine; that methionine is removed in the mature protein, so the structure numbers every residue one lower (UniProt's His 65 is His 64 in 1MBN).

2. The sequence and its properties

from io import StringIO
from Bio import SeqIO
from Bio.SeqUtils.ProtParam import ProteinAnalysis

fasta = requests.get(f"https://rest.uniprot.org/uniprotkb/{acc}.fasta").text
record = SeqIO.read(StringIO(fasta), "fasta")
p = ProteinAnalysis(str(record.seq))
print(f"molecular weight {p.molecular_weight():.0f} Da, isoelectric point {p.isoelectric_point():.2f}, GRAVY {p.gravy():.3f}")
helix, turn, sheet = p.secondary_structure_fraction()
print(f"amino acids that favour helix: {helix:.0%}")

3. The haem pocket of 1MBN

The example file 1mbn.pdb is the first protein structure solved at high resolution (1960s). The haem's iron is held by His 93 below it – the proximal histidine – while His 64 above it lines the pocket where oxygen binds.

from Bio.PDB import PDBParser

s = PDBParser(QUIET=True).get_structure("1MBN", "/data/examples/structures/1mbn.pdb")
chain = s[0]["A"]
iron = next(a for a in s.get_atoms() if a.element == "FE")
for resno in (64, 93):
    res = chain[resno]
    print(f"{res.get_resname()} {resno}: NE2 is {iron - res['NE2']:.2f} Å from the iron")

# every residue with an atom within 5 Å of the haem
haem = [r for r in chain if r.get_resname() == "HEM"][0]
near = sorted({r.id[1] for r in chain for a in r for b in haem
               if r is not haem and r.id[0] == " " and a - b < 5})
print(len(near), "residues within 5 Å of the haem:", near)

4. The AlphaFold model

AlphaFold's file keeps its confidence, pLDDT, in the B-factor column.

import matplotlib.pyplot as plt

info = requests.get(f"https://alphafold.ebi.ac.uk/api/prediction/{acc}").json()[0]
with open(f"{info['entryId']}.pdb", "w") as fh:
    fh.write(requests.get(info["pdbUrl"]).text)
model = PDBParser(QUIET=True).get_structure("AF", f"{info['entryId']}.pdb")
plddt = [res["CA"].get_bfactor() for res in model[0]["A"]]

plt.figure(figsize=(6, 2.5))
plt.plot(range(1, len(plddt) + 1), plddt)
plt.axhline(70, ls="--", c="grey")
plt.xlabel("residue (UniProt numbering)")
plt.ylabel("pLDDT")
plt.tight_layout()
plt.show()

5. Model and experiment, superposed

The Cα atoms of the model (residue n + 1) and of 1MBN (residue n) are laid on each other as closely as they go; the RMSD (root-mean-square deviation: the square root of the mean of the squared distances) then says how far apart they are.

from Bio.PDB import Superimposer

exp_ca, model_ca = [], []
af_chain = model[0]["A"]
for res in chain:
    n = res.id[1]
    if res.id[0] == " " and "CA" in res and (n + 1) in af_chain:
        exp_ca.append(res["CA"])
        model_ca.append(af_chain[n + 1]["CA"])
sup = Superimposer()
sup.set_atoms(exp_ca, model_ca)
print(f"RMSD over {len(exp_ca)} Cα atoms: {sup.rms:.2f} Å")

An RMSD around 1 Å means the prediction matches the experiment closely. To see it, run in the 3D viewer: load 1MBN; load af P02185; superpose AF-P02185-F1:A onto 1MBN:A (Viewer commands). The viewer reports a slightly smaller RMSD over a few residues fewer: it pairs the residues by aligning the sequences, leaves out pairs that end up more than 2 Å apart (loose ends) and fits again, as structure programs often do – the code above keeps every pair.

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

Biopython: Cock PJA, Antao T, Chang JT, Chapman BA, Cox CJ, Dalke A, et al. (2009) Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics 25(11): 1422–1423. doi:10.1093/bioinformatics/btp163
NGL Viewer (the 3D viewer): Rose AS, Bradley AR, Valasatava Y, Duarte JM, Prlić A, Rose PW (2018) NGL viewer: web-based molecular graphics for large complexes. Bioinformatics 34(21): 3755–3758. doi:10.1093/bioinformatics/bty419; Rose AS, Hildebrand PW (2015) NGL Viewer: a web application for molecular visualization. Nucleic Acids Research 43(W1): W576–W579. doi:10.1093/nar/gkv402

A table, in Python and in R

Are essential genes different? The example table describes 5,481 yeast genes – whether each is essential, how long its protein is, how many partners it has, how abundant it is. The same analysis in both languages: read the table, summarise it, plot it, and test what the plots suggest.

Each step has the Python code first, then the R code; each button takes the code of one language.

1. Read the table

import pandas as pd

genes = pd.read_csv("/data/examples/tables/yeast_genes.csv")
print(genes.shape)
genes.head()

genes <- read.csv("/data/examples/tables/yeast_genes.csv")
dim(genes)
head(genes[, 1:8])

One row per gene. essential is 1 when the yeast cannot live without the gene; interaction_partners counts the proteins its protein is known to interact with; human_homolog is 1 when the gene has a human counterpart.

2. Summarise by group

genes.groupby("essential")[["protein_length", "interaction_partners", "protein_abundance"]].median()

aggregate(cbind(protein_length, interaction_partners, protein_abundance) ~ essential,
          data = genes, FUN = median)

Medians, not means: these columns are skewed – a few proteins have thousands of partners. (R's formula drops the genes with a missing value in any of the columns.)

3. Plot it

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(5, 3.5))
genes.boxplot(column="interaction_partners", by="essential", ax=ax)
ax.set_yscale("log")
ax.set_xticklabels(["not essential", "essential"])
ax.set_xlabel("")
ax.set_ylabel("interaction partners")
ax.set_title("")
fig.suptitle("")
plt.tight_layout()
plt.show()

library(ggplot2)
ggplot(genes, aes(x = factor(essential, labels = c("not essential", "essential")), y = interaction_partners + 1)) +
  geom_boxplot() +
  scale_y_log10() +
  labs(x = NULL, y = "interaction partners + 1") +
  theme_bw()

On a log scale, because the values span four orders of magnitude (the + 1 in R keeps genes with no partners on the plot).

4. Test it

The Mann–Whitney (Wilcoxon rank-sum) test asks whether one group's values tend to be larger, without assuming a normal distribution.

from scipy.stats import mannwhitneyu

ess = genes.loc[genes.essential == 1, "interaction_partners"]
non = genes.loc[genes.essential == 0, "interaction_partners"]
mannwhitneyu(ess, non)

wilcox.test(interaction_partners ~ essential, data = genes)

Are genes with a human counterpart more often essential? Two yes/no columns: a table of counts and a chi-squared test.

from scipy.stats import chi2_contingency

counts = pd.crosstab(genes.human_homolog, genes.essential)
print(counts)
print(f"essential: {counts.loc[1, 1] / counts.loc[1].sum():.0%} with a human homologue, {counts.loc[0, 1] / counts.loc[0].sum():.0%} without")
chi2_contingency(counts)

counts <- table(homolog = genes$human_homolog, essential = genes$essential)
counts
prop.table(counts, margin = 1)
chisq.test(counts)

5. Two measurements together

Highly expressed genes use the codons that are read fastest (high codon bias). Spearman's correlation works on ranks, so it does not mind skewed values.

from scipy.stats import spearmanr

d = genes.dropna(subset=["codon_bias", "protein_abundance"])
print(spearmanr(d.codon_bias, d.protein_abundance))
d.plot.scatter(x="codon_bias", y="protein_abundance", logy=True, s=4, alpha=0.4, figsize=(5, 3.5))
plt.show()

cor.test(genes$codon_bias, genes$protein_abundance, method = "spearman", exact = FALSE)
ggplot(genes, aes(codon_bias, protein_abundance)) +
  geom_point(size = 0.6, alpha = 0.4) +
  scale_y_log10() +
  theme_bw()

6. Save the results

summary = genes.groupby("location").agg(genes=("orf", "size"), essential=("essential", "mean")).sort_values("essential")
summary.to_csv("essential_by_location.csv")
summary

summary <- aggregate(essential ~ location, data = genes, FUN = mean)
summary <- summary[order(summary$essential), ]
write.csv(summary, "essential_by_location.csv", row.names = FALSE)
summary

The file is written next to the notebook, in ~/notebooks. A p-value says how surprising a difference would be by chance; how large the difference is – the medians, the percentages – is what to report with it.

Databases from code

The Databases tab asks one question at a time. Code asks a hundred – and keeps the record. This recipe shows the patterns: a request and its answer as a table, many at once, a polite loop, what to do when a service is busy, NCBI with Biopython, and the first reads of a file too big to download.

1. A request that tries again

Web services are shared and sometimes busy: now and then they answer 429 (too many requests) or 5xx (a problem of the server), or nothing at all. A small function waits and asks again – and stops with a clear error for anything else (a wrong identifier gives 400 or 404).

import time
import requests
import pandas as pd

def ask(method, url, tries=4, **kwargs):
    """Make a request; when the service is busy or does not answer, wait and try again."""
    for attempt in range(1, tries + 1):
        try:
            r = requests.request(method, url, **kwargs)
            if r.status_code not in (429, 500, 502, 503, 504):
                r.raise_for_status()          # 400, 404 …: a mistake in the request – stop
                return r
            problem = f"HTTP {r.status_code}"
        except requests.exceptions.ConnectionError:
            problem = "no answer"
        if attempt == tries:
            raise RuntimeError(f"{problem} from {url} after {tries} tries – try again later")
        print(f"{problem} – trying again in {2 ** attempt} s")
        time.sleep(2 ** attempt)

def get_json(url, params=None):
    return ask("GET", url, params=params, headers={"Accept": "application/json"}).json()

gene = get_json("https://rest.ensembl.org/lookup/symbol/homo_sapiens/CYP2C19")
gene["id"], gene["seq_region_name"], gene["start"], gene["end"]

2. Many at once

Many services take a list in one request (a POST request with the list in its body). One request for five genes is quicker – and kinder to the service – than five.

symbols = ["CYP2C19", "CYP2C9", "CYP2D6", "CYP3A4", "CYP1A2"]
r = ask("POST", "https://rest.ensembl.org/lookup/symbol/homo_sapiens", json={"symbols": symbols},
        headers={"Content-Type": "application/json", "Accept": "application/json"})
genes = pd.DataFrame([{"symbol": s, "ensembl": g["id"], "chromosome": g["seq_region_name"], "start": g["start"], "end": g["end"]}
                      for s, g in r.json().items()])
genes

3. A polite loop

When there is no batch request, loop – with a short pause between requests. Here: the reviewed UniProt entry of each gene.

rows = []
for s in symbols:
    data = get_json("https://rest.uniprot.org/uniprotkb/search",
                    params={"query": f"gene_exact:{s} AND organism_id:9606 AND reviewed:true",
                            "fields": "accession,length", "format": "json"})
    for e in data["results"]:
        rows.append({"symbol": s, "uniprot": e["primaryAccession"], "length": e["sequence"]["length"]})
    time.sleep(0.2)          # the service is shared
genes.merge(pd.DataFrame(rows), on="symbol")

4. A table straight into pandas

Services that can answer as TSV or CSV are read by pandas directly:

from urllib.parse import urlencode

query = urlencode({"query": "ec:1.14.14.1 AND organism_id:9606 AND reviewed:true",
                   "fields": "accession,gene_primary,protein_name,length", "format": "tsv"})
p450 = pd.read_csv("https://rest.uniprot.org/uniprotkb/search?" + query, sep="\t")
print(len(p450), "human enzymes with EC 1.14.14.1")
p450.head()

5. NCBI with Biopython

Biopython's Bio.Entrez speaks NCBI's E-utilities: search, then fetch. NCBI asks who is asking – put your own address in Entrez.email.

from Bio import Entrez, SeqIO

Entrez.email = "your.name@example.ac.uk"
found = Entrez.read(Entrez.esearch(db="nuccore", term="CYP2C19[gene] AND human[orgn] AND refseq[filter] AND mRNA[filter]"))
record = SeqIO.read(Entrez.efetch(db="nuccore", id=found["IdList"][0], rettype="gb", retmode="text"), "genbank")
cds = next(f for f in record.features if f.type == "CDS")
print(record.id, record.description)
print(cds.qualifiers["protein_id"][0], cds.extract(record.seq).translate()[:60], "…")

6. The first reads of a big file

A sequencing run is usually far too big for the browser – the reads of SRR098401 are 7.5 GB per file. A Range request downloads only the start of a compressed file, and zlib unpacks as much as arrived:

import zlib

url = "https://ftp.sra.ebi.ac.uk/vol1/fastq/SRR098/SRR098401/SRR098401_1.fastq.gz"
r = ask("GET", url, headers={"Range": "bytes=0-4999999"})      # the first 5 MB
text = zlib.decompressobj(16 + zlib.MAX_WBITS).decompress(r.content).decode()
lines = text.split("\n")
keep = (len(lines) - 1) // 4 * 4                                # whole reads only
with open("SRR098401_start_1.fastq", "w") as fh:
    fh.write("\n".join(lines[:keep]) + "\n")
print(keep // 4, "reads")

!head -n 8 SRR098401_start_1.fastq

The first reads of a run come from the edge of the flow cell, and are often poor (all N): for a fair sample, use reads from the middle of the file, or trim and filter first (Reads).

7. In R

jsonlite::fromJSON() reads a URL's JSON into lists and data frames; readLines() reads text. R in the browser can make only GET requests: for POST requests (batches, BLAST, VEP with positions), use Python.

library(jsonlite)

symbols <- c("CYP2C19", "CYP2C9", "CYP2D6")
genes <- do.call(rbind, lapply(symbols, function(s) {
  g <- fromJSON(paste0("https://rest.ensembl.org/lookup/symbol/homo_sapiens/", s, "?content-type=application/json"))
  Sys.sleep(0.2)
  data.frame(symbol = s, id = g$id, chromosome = g$seq_region_name, start = g$start, end = g$end)
}))
genes

Which services work from the browser

A web page – and so a notebook here – can only ask services that allow requests from web pages (CORS): Ensembl, UniProt, NCBI's E-utilities, the EBI services (ENA, InterPro, QuickGO, BLAST …), RCSB PDB, AlphaFold DB, STRING, Reactome, Europe PMC, PubChem and the UCSC API do. Many others – KEGG, most web pages – do not: the request fails with no answer. Download such files with your browser and upload them (Your files).

References

The papers to cite when you use these programs in your work – with the version you used (Programs and versions):

Biopython: Cock PJA, Antao T, Chang JT, Chapman BA, Cox CJ, Dalke A, et al. (2009) Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics 25(11): 1422–1423. doi:10.1093/bioinformatics/btp163
Ensembl VEP: McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GRS, Thormann A, et al. (2016) The Ensembl Variant Effect Predictor. Genome Biology 17: 122. doi:10.1186/s13059-016-0974-4

Keep and share your work

Everything you run in the workbench can become a notebook – a record that shows what you did and runs again. Your files live in this browser: download what you want to keep.

Everything as a notebook

Where you worked	How it becomes a notebook
Notebook	it is one: Download → Jupyter notebook (.ipynb) (or a `.py`/`.R` script, or a web page with the results)
Terminal	Save session → Open as a notebook here or Download as a Jupyter notebook: each command in a `%%bash` cell, with its output; or `export-session NAME.ipynb`
Galaxy	the history's ⋯ → Open the history as a notebook (or download it): each upload and job as a cell that makes the same dataset again
Databases	Your queries → Save as a Python notebook (or R); Send to notebook for one query
AI agent	each run is `agent-record.ipynb` in its folder
Help recipes	Open this recipe as a notebook; Add to notebook on any code box
3D viewer	not a notebook: Save session keeps the structures, styles and views; the command line keeps most steps as commands – copy them into a script that makes the picture again

Download

One notebook: Download in the Notebook tab.
One file or folder: its ⋯ in the Files tab, or download FILE in the terminal (a folder becomes a .zip).
Everything: My files → Download all my files (.zip). Do this at the end of each session if you work on a shared computer – see where your files are kept.

Before you share or hand in a notebook

Restart the language and Run all: the notebook must run from top to bottom without the cells you deleted on the way.
Put a text cell at the top: what the notebook does, with which data, when.
Keep the data with it: copy the input files into a folder of your own, next to the notebook, and download the folder.
Record the versions (below) and the dates of database queries.
Never leave a key, a password or personal data in a notebook – also not in its outputs.

The versions

Results depend on the versions of the programs and packages. A last cell that prints them makes a notebook complete:

import sys, numpy, pandas, matplotlib
print("Python", sys.version.split()[0], "· numpy", numpy.__version__, "· pandas", pandas.__version__, "· matplotlib", matplotlib.__version__)

%%bash
samtools --version | head -n 1
bcftools --version | head -n 1
echo "minimap2 $(minimap2 --version)"

sessionInfo()

Running a notebook elsewhere

A downloaded notebook opens in Jupyter, JupyterLab, VS Code and Google Colab (R notebooks need the R kernel, IRkernel). What the other computer needs:

The data, at the paths the notebook uses. /data/examples exists only here: download it (Files tab: the ⋯ of Examples (read-only) → Download as a .zip) and change the paths.
The packages: pip install pandas matplotlib biopython requests scipy; in R, install.packages(c("ggplot2", "jsonlite")), and from Bioconductor BiocManager::install(c("DESeq2", "edgeR")).
The programs of %%bash cells – with conda, from Bioconda, in the versions used here (or the nearest ones Bioconda has):

conda create -n workbench -c conda-forge -c bioconda \
    samtools=1.17 minimap2=2.22 bowtie2=2.4.2 fastp=0.20.1 \
    seqtk=1.4 bedtools=2.31.0 kalign3=3.3.2 muscle=3.8.1551 fasttree=2.1.11
conda create -n workbench-bcftools -c conda-forge -c bioconda bcftools=1.10
conda create -n workbench-freebayes -c conda-forge -c bioconda freebayes=1.3.6
conda create -n workbench-bwa -c conda-forge -c bioconda bwa=0.7.19 bwa-mem2=2.2.1 hisat2=2.2.2 samtools
conda create -n workbench-rnaseq -c conda-forge -c bioconda subread=2.1.1 kallisto=0.52.0 salmon=2.8.0
conda create -n workbench-stringtie -c conda-forge -c bioconda stringtie=3.0.3
conda create -n workbench-mosdepth -c conda-forge -c bioconda mosdepth=0.3.14
conda create -n workbench-cutadapt -c conda-forge -c bioconda cutadapt=5.2
conda create -n workbench-somatic -c conda-forge -c bioconda lofreq=2.1.5 somatic-sniper=1.0.5.0
conda create -n workbench-strelka -c conda-forge -c bioconda strelka=2.9.10
conda create -n workbench-meta -c conda-forge -c bioconda kraken2=2.17.2 bracken=3.1
conda create -n workbench-annotation -c conda-forge -c bioconda prodigal=2.6.3 hmmer=3.4 bustools=0.45.1
conda activate workbench

bcftools 1.10 needs an older htslib than samtools 1.17, so it has an environment of its own (or take a current bcftools into the first one – its output can differ in details); Bioconda's BWA 0.7.19 needs a newer zlib than samtools 1.17 and FreeBayes 1.3.6, so it goes with a current samtools. featureCounts is in the package subread. StringTie 3.0.3 and mosdepth 0.3.14 each need their own newer HTSlib, and so an environment each. Cutadapt 5.2 needs Python 3.13, where Bowtie 2 2.4.2 needs Python 3.9 or older; LoFreq 2.1.5 needs samtools 1.24 or newer; Strelka 2.9.10 needs Python 2.7 – so each has an environment of its own. Kraken2 and Bracken share one (Kraken2 2.17.2 needs zlib 1.3.2 or newer); a Kraken2 database is not in the package – build one, or download a ready-made one. Prodigal, HMMER and bustools share one too; on a computer HMMER reads .gz files and uses several cores (--cpu). Bioconda has no Kalign 3.3.1; 3.3.2 is the nearest.

Lines that only the workbench understands – igv, open, download – do nothing useful elsewhere: leave them out.

Reference

Facts to look up: the programs and their versions, the keyboard shortcuts, what is different in the browser, what to do when something goes wrong, what happens to your data, and whose work the workbench is built from.

Programs and versionsEvery program, language and library, with its version
Keyboard shortcutsThe keys of every tab
What is different in the browserMemory, speed, storage, internet, what cannot run
When something goes wrongCommon problems and what to do
Your data and privacyWhat is kept where, and what is sent to whom
Credits and licencesThe software and data the workbench is built from

Programs and versions

Results depend on versions: cite the ones you used, with the papers of the programs (To cite; in full at the end of each program's help page). help in the terminal prints this list too; NAME --version one program's.

The terminal's bioinformatics programs

The real programs, compiled to WebAssembly (most by the Biowasm project; BWA, BWA-MEM2, HISAT2, FreeBayes, featureCounts, StringTie, kallisto, Salmon, mosdepth, Cutadapt, MACS3, LoFreq, SomaticSniper, Strelka2, Kraken2, Bracken, bustools, Prodigal, HMMER and PLINK 2 by the genomics-wasm project). Galaxy runs the same ones.

Program	Version	For	To cite
fastp	0.20.1	read quality control and trimming	Chen et al. 2018
Cutadapt (with dnaio 1.2.4 and xopen 2.1.0)	5.2	removing adapters and primers; quality trimming. It runs in a Python of its own: Pyodide 314.0.7, CPython 3.14.2	Martin 2011
seqtk	1.4	FASTA/FASTQ conversion, sampling, statistics	github.com/lh3/seqtk
wgsim	0.3.1-r13	simulating reads	github.com/lh3/wgsim
minimap2	2.22	mapping reads (short and long), aligning sequences	Li 2018; Li 2021
Bowtie 2 (`bowtie2`, `bowtie2-build`)	2.4.2	mapping short reads	Langmead and Salzberg 2012
BWA (`bwa index`, `bwa mem`)	0.7.19	mapping short reads	Li 2013; Li and Durbin 2009
BWA-MEM2 (`bwa-mem2 index`, `bwa-mem2 mem`)	2.2.1	mapping short reads (bwa mem's algorithm, rewritten for speed)	Vasimuddin et al. 2019
HISAT2 (`hisat2`, `hisat2-build`, `hisat2-inspect`)	2.2.2	mapping RNA-seq reads across introns	Kim et al. 2019
featureCounts (Subread)	2.1.1	counting RNA-seq reads per gene (several threads)	Liao et al. 2014
StringTie	3.0.3	assembling transcripts from RNA-seq alignments, and their abundance	Shinder et al. 2026; Pertea et al. 2015; Kovaka et al. 2019
kallisto	0.52.0	transcript abundances by pseudoalignment	Bray et al. 2016
Salmon	2.8.0	transcript abundances by selective alignment (several threads)	Patro et al. 2017
bustools	0.45.1	single-cell RNA-seq: sorting, correcting and counting the BUS records of `kallisto bus` into a matrix of cells and genes	Melsted et al. 2021; Melsted et al. 2019; Sullivan et al. 2025
MACS3 (`macs3 callpeak`, `hmmratac`, `bdgcmp` …; with hmmlearn 0.3.3)	3.0.5	peaks of ChIP-seq (with or without a control), ATAC-seq (HMMRATAC) and the like. It runs in the Python of Cutadapt, with NumPy 2.4.6, SciPy 1.18.0, scikit-learn 1.8.0, pandas 3.0.2, h5py 3.13.0 and AnnData 0.10.9; one process	Doherty et al. 2026; Zhang et al. 2008
SAMtools	1.17	SAM/BAM files: sort, index, view, statistics, depth, faidx	Danecek et al. 2021; Li et al. 2009
mosdepth	0.3.14	read depth per base, window or region	Pedersen and Quinlan 2018
BCFtools	1.10	variant calling and VCF/BCF files	Danecek et al. 2021
FreeBayes	1.3.6	variant calling (haplotype-based)	Garrison and Marth 2012
LoFreq (`lofreq call`, `indelqual`, `somatic` …)	2.1.5	variant calling at low frequencies; somatic variants of a tumour and its normal	Wilm et al. 2012; Albers et al. 2011
SomaticSniper (`bam-somaticsniper`)	1.0.5.0	somatic SNVs of a tumour and its normal	Larson et al. 2012
Strelka2 (`configureStrelkaSomaticWorkflow.py`, then `runWorkflow.py`)	2.9.10	somatic SNVs and indels of a tumour and its normal	Kim et al. 2018
Kraken2 (`kraken2`, `kraken2-build`, `kraken2-inspect`)	2.17.2	assigning reads to taxa (metagenomics); building and inspecting databases. `kraken2 --version` says 2.17.1, as the release does	Wood et al. 2019; Lu et al. 2022
Bracken (`bracken`, `bracken-build`)	3.1	the abundance of species (or genera …) from a Kraken2 report. Its Python programs run in the Python of Cutadapt; `bracken -v` says 3.0.1, as the release does	Lu et al. 2017
Prodigal	2.6.3	finding the genes that code for proteins in genomes of bacteria and archaea, and in metagenomes	Hyatt et al. 2010
HMMER (`hmmbuild`, `hmmsearch`, `hmmscan`, `hmmpress`, `phmmer`, `jackhmmer`, `nhmmer`, `nhmmscan`, `hmmalign`, `hmmfetch`, `hmmstat`, `hmmemit`, `hmmconvert`)	3.4	profile hidden Markov models: protein domains and families (as Pfam's), searches of proteins and DNA (one thread here)	hmmer.org; Eddy 2011
PLINK 2	2.0.0-b.1-dev (12 October 2026)	the genotypes of many people: quality control, allele frequencies, LD, kinship, association tests (GWAS) and scores (threads: a cross-origin isolated page; no `--pca` in this build)	Chang et al. 2015; cog-genomics.org/plink/2.0
HTSlib (`bgzip`, `tabix`)	1.17	compressing and indexing	Bonfield et al. 2021
bedtools	2.31.0	genomic intervals	Quinlan and Hall 2010
gffread	0.12.7	GFF/GTF files; transcript and protein sequences	Pertea and Pertea 2020
Kalign	3.3.1	multiple sequence alignment	Lassmann 2020
MUSCLE	3.8.1551	multiple sequence alignment	Edgar 2004
FastTree	2.1.11	phylogenetic trees	Price et al. 2010
seq-align (`needleman_wunsch`, `smith_waterman`, `lcs`)	2017.10.18	pairwise alignment	github.com/noporpoise/seq-align; Needleman and Wunsch 1970; Smith and Waterman 1981
LASTZ	1.04.52	aligning long DNA sequences	Harris 2007

Other tools of the terminal

Tool	Version
GNU coreutils (`sort`, `cut`, `head`, `wc` …)	8.32
GNU grep (`grep -P`: 3.11 with PCRE2 10.42)	3.7
GNU sed	4.8
GNU awk (`gawk`, `awk`)	5.1.0
GNU diffutils (`diff`, `cmp`)	3.10
jq	1.7
Graphviz (`dot`), through Viz.js 3.31.0	16.1.0
Snakemake, browser edition: the commonly used core of Snakemake, re-written to run here	9.27.0

The shell itself – bash's language, and ls, cp, curl, unzip … – is written for the workbench, to behave as on Linux.

Languages

Language	Version	With
Python (Pyodide 0.29.5)	3.13.2	NumPy 2.2.5, pandas 2.3.3, Matplotlib 3.8.4 when it starts; on import: Biopython 1.85, SciPy 1.14.1, scikit-learn 1.7.0, statsmodels 0.14.4, NetworkX 3.4.2, requests 2.32.4 and the other Pyodide packages; `%pip install`: pure-Python packages from PyPI
R (webR 0.6.0)	4.6.0	packages on first use, from the webR repository (CRAN) and r-universe (Bioconductor, the current release: DESeq2 1.52, edgeR 4.10, limma 3.68 in October 2026), and locfit 1.5-9.12, built for this site – see R

In a notebook, %pip list (Python) and sessionInfo() (R) show what is loaded.

In the other tabs

Tab	Made with
3D viewer	NGL Viewer 2.5.0
Genome browser	igv.js 3.8.9
Networks	Cytoscape.js 3.34.3, with the layouts fCoSE 2.2.0 and dagre (cytoscape-dagre 4.0.1); the analyses (statistics, hubs, MCODE, MCL, Louvain) written for the workbench and checked against NetworkX – see Network analysis
Files, notebook editor	CodeMirror 5.65.18
AI assistant	marked 18.0.13 and DOMPurify 3.4.16 (to show answers safely); the model is the one you connect
Galaxy	written for the workbench, in the style of Galaxy; its tools use the programs above, under the IDs and versions of the Galaxy Tool Shed

Their licences: Credits and licences.

Keyboard shortcuts

On a Mac, use ⌘ (Cmd) where it says Ctrl – except in the terminal: there Ctrl+C and Ctrl+L are the Ctrl key, as in the Mac's own Terminal (⌘+C copies).

Everywhere

`Ctrl`+`\`	hide or show the help: the workbench fills the window
`←` `→` on the tab bar	the next or previous tab
`←` `→` on the border between help and workbench	make the help narrower or wider (or drag the border)
`Esc`	close a dialog or menu

Help

`Enter` in the search box	open the first result
`↓` `↑`	move through the results
`Esc`	clear the search

Notebook – while editing a cell

`Shift`+`Enter`	run the cell, go to the next
`Ctrl`+`Enter`	run the cell, stay
`Alt`+`Enter`	run the cell, add a new one below
`Ctrl`+`S`	save now
`Ctrl`+`/`	comment or uncomment the lines
`Tab` / `Shift`+`Tab`	indent / unindent
`Esc`	stop editing: the cell stays selected

Notebook – a selected cell

`Enter`	edit it
`A` / `B`	a new cell above / below
`M` / `Y`	make it a text / code cell
`D` `D`	delete it
`↑` `↓` (or `K` `J`)	the cell above / below

Terminal

`↑` `↓`	earlier commands
`Tab`	complete a command or file name
`Ctrl`+`C`	while a command runs: the program that is running finishes, the rest is skipped; a second time: the program is stopped by force (what it was writing is lost). At the prompt: clear the line
`Ctrl`+`L`	clear the screen
`Home` / `End`	start / end of the line

To copy text from the terminal, select it with the mouse and press Ctrl+C – with text selected it copies instead of stopping.

Files

`Ctrl`+`S`	save the file (files are not saved by themselves: a dot on the tab of a file marks changes that are not saved)
`Ctrl`+`/`	comment or uncomment the lines
`↑` `↓` `Enter` in the list	move through the files, open one
`Shift`+`F10` in the list	the ⋯ menu of a file

3D viewer

`Esc`	clear the selection; stop measuring or box-selecting
`Delete`	remove the selected label or shape
`↑` `↓` in the command line	earlier commands

Networks

The keys work once you have clicked the network.

`Shift`+click, `Shift`+drag	add to the selection; select the nodes inside a box
Double-click a node	zoom to it and its neighbours
Right-click	the menu of a node, an edge or the background
`Ctrl`+`A`	select every node
`Esc`	select nothing
`Delete`	delete what is selected
`Ctrl`+`Z`, `Ctrl`+`Shift`+`Z` (or `Ctrl`+`Y`)	undo, redo
`F`	fit the whole network in view
`+` `−`	zoom in, out (or the mouse wheel)
`Enter` in Find nodes	select the nodes found and zoom to them

IDE: RStudio

The full list is in RStudio’s Help ▸ Keyboard Shortcuts Help, and on the page RStudio.

`Ctrl`+`Enter`	run the current line, or the selected lines, in the Console
`Ctrl`+`Shift`+`S`	source the whole script
`Ctrl`+`Shift`+`K`	knit or render the document
`Ctrl`+`Alt`+`I`	insert a chunk into an R Markdown document
`Esc` in the Console	stop the command that is running

IDE: Spyder

The full list is in Spyder’s Help ▸ Keyboard shortcuts, and on the page Spyder.

`F5`	run the whole file
`Ctrl`+`Enter` / `Shift`+`Enter`	run the current cell / run it and go to the next one
`F9`	run the selected lines, or the current line
`F12`	set or clear a breakpoint
`Ctrl`+`F5`	debug the file; then `Ctrl`+`F10` steps, `Ctrl`+`F12` continues and `Ctrl`+`Shift`+`F12` stops
`Ctrl`+`I`	show the documentation of the name at the cursor
`Ctrl`+`C` in the console	stop the code that is running

AI assistant

`Enter`	send
`Shift`+`Enter`	a new line in the message

What is different in the browser

Everything in the workbench runs on your own computer, inside this web page: the programs are the real ones, compiled to WebAssembly, and Python and R are the real languages. That brings some limits a server would not have.

Speed and memory

Each program runs on one processor core; options for threads (-t, -p, -@, -j, --threads, --cpu) are ignored, with a note. Three programs are different: featureCounts counts on up to four threads (-T 1 to -T 4), Salmon runs its own helper threads, and PLINK 2 reads its files on threads of their own and computes on one or two (--threads 2 at most) – see programs on several threads. Programs run more slowly than on a server, and the speed depends on your computer: keep the data small.
A program, Python and R each have a few gigabytes of memory at most (WebAssembly is 32-bit: 4 GB is the ceiling, and browsers often allow less). Mapping reads to a whole human genome does not fit; regions, a few genes, bacterial and viral genomes do.
Python has no threads, no multiprocessing and no subprocess; input() does not work in a notebook, though it does in Spyder (in the IDE tab). R cannot run system(): use a %%bash cell.

Programs on several threads

featureCounts, Salmon and PLINK 2 start threads that share the program's memory. A browser allows that only in a page that is cross-origin isolated: the page and its workers come with two headers, Cross-Origin-Opener-Policy and Cross-Origin-Embedder-Policy, which also keep pictures of other sites out unless those sites allow them. Static hosting (GitHub Pages) cannot send the headers, so the workbench adds them itself with a small service worker (isolation-sw.js; it stores nothing): the first time you open the workbench, the page loads once more, at once. serve.py sends the headers itself.

If featureCounts, Salmon or PLINK 2 say that the page is not isolated, reload the page. It cannot be isolated when it was not opened over https (or from localhost), or when the browser does not allow service workers there – some private windows do not.
The other programs, Python and R do not need isolation, and work either way. In the IDE tab, Spyder’s debugger and input(), and stopping code while it runs, need it too (see Spyder).
To turn isolation off – if a picture of another site does not show, say – open the workbench with ?isolation=off at the end of its address; open it without that to turn it on again.

Programs that are partly here

Cutadapt is a Python program: each command starts a Python of its own (Pyodide 314.0.7, CPython 3.14.2), which takes a few seconds, and the first one loads about 13 MB. It trims on one core.
LoFreq: lofreq call-parallel and lofreq vcfplot are not here – use lofreq call. lofreq somatic runs its steps one after another.
Kraken2 works with small databases: the database must fit in the browser's memory – a few hundred megabytes at most. The databases used for real samples (Kraken2's standard database is tens of gigabytes) do not; kraken2-build cannot download (--download-taxonomy, --download-library, --standard, --special) and has no --add-to-library – a database is built from sequences, a map and a taxonomy that you put in its folder (see Metagenomics). Reads compressed with bzip2 cannot be read. Kraken2 and Bracken run on one thread.
MACS3 is a Python program too: each command starts the same Python, with NumPy, SciPy, scikit-learn and pandas – a few seconds; the first command loads about 45 MB. It runs as one process (callvar -m N works through its N parts one after another, with the same result) and keeps the reads in memory: the reads of a chromosome or of a small genome are fine, a whole human experiment of tens of millions of reads is too much for a browser tab. Its messages say [RSS unavailable] where Linux gives the memory it uses.
Bracken: its two Python programs start a Python of their own for each command, as Cutadapt does (a few seconds). Where Bracken's script reports an error but ends with status 0 – a database or a report that is not there – it ends with status 2 here.
HMMER runs on one thread: --cpu is left out, with a note. It cannot read files compressed with gzip (.gz) – HMMER starts the program gzip to read them, and a web page cannot start programs: unpack them first (gunzip -k FILE.gz). Easel's own programs (esl-reformat, esl-sfetch …) are not here. HMMER needs WebAssembly SIMD, which all current browsers have. Its line “# CPU time” says 0.00u 0.00s: a web page does not tell a program how much processor time it used (the “Elapsed” time is right). Profile databases as big as the whole of Pfam do not fit in the browser's memory: use the InterPro card of the Databases tab for those.
bustools sorts on one thread (sort -t 4 becomes -t 1, with a note) and in 64 MB of memory unless told otherwise – 1 GB at most: a larger -m becomes -m 64M, with a note. Its matrices are those of bustools on Linux; the .ec.txt file that count writes without --genecounts has the same equivalence classes, but those the matrix does not use can be numbered differently.
PLINK 2 is built without the LAPACK library, so --pca and --pheno-svd are not available (PLINK says so and stops). It works in 640 MiB of memory, the least it accepts: the terminal adds --threads 1 --memory 640 require unless you give them, and a larger --memory becomes 1024 (more threads than 2 become 2), with a note. --script and --rerun, which read options from a file, are not available: plink2 $(cat FILE) does the same. Its genotypes must fit in that memory – thousands of people and variants, not a whole biobank. It needs a cross-origin isolated page.
Strelka2: the somatic workflow, not the germline one. Strelka's workflow scripts are Python 2, which is not here: the terminal runs runWorkflow.py itself, with Strelka's own programs – -m local only, the tasks one after another. python2 runs these scripts and nothing else.

Files

Your files are kept in this browser on this computer, not on a server – see where your files are kept.
Browsers give a site a share of the free disk space – usually several gigabytes; My files shows how much is used. Files of hundreds of megabytes work, but slowly. Python and R are not given files over 300 MB, and a file over 1 GB is not kept in the browser at all: it is there while the page is open – download it if you need it.
Large intermediate files fill the space fast: pipe programs together (minimap2 … | samtools sort -o mapped.bam -) instead of writing a SAM file, and delete what you no longer need.

Installing

No sudo, apt or background jobs (&). Programs cannot be added to the terminal, and there is no need to: every program is installed already, so there is no conda here either. (To set up the same programs on another computer, see Keep and share your work.)
Python packages install with %pip install (pure-Python packages, and those Pyodide has built); R packages with install.packages() (those built for webR). Both last until the page is reloaded – except in RStudio, which remembers the R packages you installed there and fetches them again when library() needs them.
Some packages have no build for the browser yet – mostly those that need system libraries; the error says so. (DESeq2 and edgeR work: locfit, which they need, comes with the workbench.)

The internet

Python is downloaded when a Python notebook is open, R when an R notebook is open – and the page opens the last notebook you had open when it starts (a Python one on the first visit) –, and again for Spyder or RStudio when you first open one of them in the IDE tab; packages when they are imported or installed: tens of megabytes, of which the browser keeps a copy for later visits.
A web page can only ask web services that allow requests from web pages (CORS). The big bioinformatics services do – Ensembl, UniProt, NCBI, EBI, the PDB, AlphaFold DB and the others of the Databases tab; KEGG and most web pages do not. This holds for notebooks, curl and wget alike.
R in the browser can only make GET requests; use Python for POST requests.
The Databases tab, the AI assistant, the hosted genomes of the genome browser, downloads from the PDB, and the Networks tab's STRING, IntAct, NDEx and enrichment need the internet all the time.

The tabs

Galaxy offers a selection of tools, runs one job at a time, and has no user accounts or shared histories.
Snakemake is a browser edition of its commonly used core: rules, wildcards, expand, -n, --dag, --summary; no cluster execution.
The IDE tab: RStudio and Spyder are imitations, made for teaching, of the programs of the same names. Each has an R or a Python of its own, separate from the notebooks’. Reports are knitted to HTML only (pandoc and LaTeX, which PDF and Word need, do not run in a browser); Shiny apps need the page’s service worker; there are no projects, no Git and no plug-ins (The IDE tab).
The 3D viewer needs WebGL (3D graphics in the browser): if it does not start, turn on hardware acceleration in the browser's settings.
The Networks tab draws networks of a few thousand nodes smoothly; larger ones (up to 60,000 nodes and 300,000 edges) open, but are slow to draw and to lay out. Statistics and hub rankings follow every shortest path: seconds for a few thousand nodes, about half a minute for 10,000. Cytoscape's apps (plugins) are not here; the most used of them – NetworkAnalyzer, cytoHubba, MCODE, MCL, enrichment – have their equivalents in the tab (Network analysis).
The workbench is made for a laptop or desktop screen; on a phone or tablet it works, cramped.

Use a current Chrome, Edge, Firefox or Safari.

When something goes wrong

The usual problems, and what to do about them. Error messages in the notebook and the terminal often carry a hint of their own – and Ask the AI assistant about this error explains them (AI assistant).

Starting

Python or R takes long to start: The first start of a visit downloads the language (tens of megabytes); later visits use the browser's copy and are quicker. The label at the top right of a notebook says what it is doing.
A tab says it could not start: Reload the page (your files are kept). For the 3D viewer: turn on hardware acceleration in the browser's settings. Use a current Chrome, Edge, Firefox or Safari.
The page is slow or stops answering: A long computation is running: Stop it (notebook), or press Ctrl+C in the terminal – once, and what follows the running program is skipped; twice, and the program itself is stopped. Close other heavy tabs of the browser.

Files

My files are gone: They are kept in one browser on one computer – is it the same browser, the same computer, the same web address, and not a private window? If the browser's data for the site was cleared, they cannot be recovered: download what you want to keep (Keep and share your work).
"No such file or directory": The file is not where the code looked. A notebook works in its own folder (~/notebooks): use full paths such as /data/examples/… or /home/user/… (Python's open() does not understand ~). %pwd or getwd() shows the folder; the Files tab shows the names.
"Read-only file system" or "Permission denied": /data/examples is read-only: its files cannot be changed, and nothing new can be written into it – not even what a program writes next to its input (samtools faidx writes its index next to the FASTA file). Copy what you need into your home folder and work there.
A PDF "cannot be read as a PDF": The Files tab draws PDF files itself, so no PDF plug-in of the browser is needed. A file it cannot read is damaged, is not a PDF, or is not finished: a plot that R or Python is still drawing into is written only when it is closed – dev.off() after pdf("plot.pdf") in R, plt.savefig("plot.pdf") or the end of with PdfPages(…) in Python. Download the file to open it with another program.
The space is running out: My files (top right) shows how much is used. Delete large intermediate files – SAM files, unsorted BAM files, trimmed reads – or download and delete them.

Programs and code

"command not found": The program is not in the workbench: help lists what is. Python packages: %pip install NAME; R packages: install.packages("NAME").
featureCounts, salmon or plink2: "… this page is not" cross-origin isolated: These three programs work on several threads, which a browser allows only in a cross-origin isolated page. Reload the page: the workbench isolates it as it starts. It cannot when the page was not opened over https (or from localhost), or in a browser window that allows no service workers (some private windows) – see programs on several threads.
A variant caller says an index is missing, or that the chromosomes do not match: Index the BAM files (samtools index) and the reference (samtools faidx) – again after a program has written a new BAM file, as lofreq indelqual does. The reference must be the one the reads were mapped to, with the sequence names of the BAM header (samtools view -H file.bam | grep @SQ). See Somatic variants.
bracken: "ERROR: …/database100mers.kmer_distrib does not exist": Bracken needs the distribution for the length of your reads, and it takes 100 when -r is not given. Give the length: -r 151 for the example data (ls DB/*.kmer_distrib shows which lengths a database has); bracken-build -l LENGTH makes one for another length. See Metagenomics.
Kraken2 leaves most reads unclassified, or names a species that cannot be there: Kraken2 finds only what is in its database: a read of an organism that is not there stays unclassified, or goes to a relative that is – E. coli reads to Shigella, say. Look at what the database holds (kraken2-inspect), raise --confidence, and let Bracken's threshold (-t) leave out species with few reads. See Metagenomics.
prodigal: "Sequence must be 20000 characters": In its usual mode (-p single) Prodigal learns the genome from the sequences it is given, and fewer than 20,000 bases are too few to learn from. Give it the whole genome – all its contigs in one file – or, for a short sequence or a metagenome, use -p meta.
hmmscan: "Failed to open binary auxfiles … use hmmpress first": hmmscan searches a collection of profiles that hmmpress has prepared: run hmmpress FILE.hmm first, in a folder of your own – it writes four files next to the profiles, which it cannot do in /data/examples. When hmmpress says "Delete old hmmpress indices first", the collection is pressed already; after changing its profiles, press it again with hmmpress -f FILE.hmm.
hmmsearch or hmmscan: "File format problem in trying to open HMM file": The file given as the profile is not one – an alignment, or sequences. Build the profile from the alignment first (hmmbuild NEW.hmm ALIGNMENT.sto), and mind the order: the profiles come first and the sequences second, for hmmsearch and hmmscan alike.
plink2: "--pca requires plink2 to be built with LAPACK": This build of PLINK 2 has no LAPACK, the library its principal components (--pca) and --pheno-svd need: they are not available here. The rest of PLINK works. See Genotypes and association.
plink2 --glm: "invoked without --covar … usually an analytical mistake": PLINK asks for covariates in an association test – ancestry (principal components), sex, age – because a test without them is usually confounded. Give them with --covar FILE, or, when you are sure, add the modifier: --glm allow-no-covars.
plink2: "Failed to open NAME.psam" (or .pgen, .pvar, .bed …): --pfile NAME reads three files, NAME.pgen, NAME.pvar and NAME.psam (--bfile NAME: .bed, .bim, .fam): give the name without its ending, and make the files first (plink2 --vcf FILE.vcf.gz --make-pgen --out NAME).
macs3 takes several seconds before it does anything: Each command of MACS3 starts a Python of its own with NumPy, SciPy, scikit-learn and pandas – the first time about 45 MB are loaded. After a command, the next one's Python gets ready in the background: a command typed a little later starts at once.
macs3 callpeak: "MACS3 needs at least 100 paired peaks at + and - strand to build the model": MACS3 works out the length of the fragments from strong peaks, and found too few: the reads are few, the enrichment is weak, or the genome size is far too large for the reads (-g hs for the reads of one chromosome). Give the fragment length yourself – --nomodel --extsize 200, with the length of the library's fragments – or use the reads as pairs (-f BAMPE, -f BEDPE): MACS3 then sees the fragments. See Peaks.
The peaks of MACS3 do not show in the genome browser: Choose the genome the reads were mapped to – hg19 for the example ChIP-seq (igv --genome hg19 …), sacCer3 for the yeast ATAC-seq: the chromosome names of the peaks must be those of the genome.
A program says it is out of memory, or stops without a message: The data is too big for the browser: take a part of it (samtools view of a region, seqtk sample) – see What is different in the browser.
Python or R forgot my variables: The language was restarted – by Restart, Stop, or a reload of the page. Run the notebook's cells again from the top (Run all).
The terminal does not let me type: The AI agent is working in it: wait, or press Stop in the AI assistant tab.

The internet

"Failed to fetch", "no answer", a network error: Either the internet connection is off, or the site does not allow requests from web pages (CORS – KEGG and most web pages do not). Download such files with your browser and upload them in the Files tab.
HTTP 429, "too many requests": Wait a minute; in code, pause between requests (Databases from code).
HTTP 400 or 404: The service did not understand the query, or there is no such identifier: check the spelling and the version (NM_000769.4).
HTTP 500 or 503: The service is busy or down – try again later.

The other tabs

The genome browser shows nothing: Is the right genome chosen? The chromosome names of the file must be those of the genome (human_CYP2C19 for the example genome, chr10 for hg19/hg38). A BAM file needs its .bai (samtools index), a FASTA genome its .fai (samtools faidx). Reads appear only when you zoom in.
A Galaxy dataset is red: Open its details (): the program's messages say what went wrong. Often an input of the wrong kind, or a setting; runs the tool again with the same settings, to change one.
Networks: STRING does not know some of my genes: STRING knows proteins: long non-coding RNAs and pseudogenes have none, and Ensembl IDs that Ensembl has retired are not found. Check the species. Versions at the end of Ensembl IDs (ENSG00000120129.5) are taken off before the IDs are sent.
Networks: few nodes found a row of my table: The two columns must hold the same kind of names – gene names with gene names, Ensembl IDs with Ensembl IDs. A network from STRING keeps the names you sent in its column query, and STRING's own names in name: choose the one that matches your table.
Networks: the network is a hairball: Raise STRING's confidence (0.7), keep fewer genes (a smaller p-value), hide the nodes without edges (Select → Nodes without edges), or look at the clusters one at a time (Analyse → Clusters, then new network on a cluster).
RStudio or Spyder takes a long time to open, or says it could not be opened: The first time, the IDE downloads R or Python, which can take a minute or so on a slow connection. If it fails, check your internet connection and press Try again, or reload the page; your files are kept.
My variable is missing in RStudio or Spyder (or in the notebook): Each IDE has an R or a Python of its own, separate from the notebooks’, so a variable only exists where you created it. Run the code that creates it in the IDE as well (The IDE tab).
Spyder: input() or the debugger doesn’t work: Both need the page to be cross-origin isolated, and the console says so when it starts if it isn’t. Reload the page, or open it in an ordinary (not private) window (Spyder).
RStudio: "there is no package called …": In RStudio, as on your own computer, a package has to be installed once before library() can load it: run install.packages("NAME"), or use Tools ▸ Install Packages… (RStudio).
The AI assistant says "Not connected", the key is not accepted, or it is busy: Check the service and key in AI settings and press Test. Free keys allow only a few requests per minute: wait a minute.

Reporting a problem

Say which browser and computer you use, what you did, and what happened – with the exact message (copy it, or take a screenshot). Leave out any personal data.

`Shift`+`Enter` or ▶ Run	run the cell and go to the next one
`Ctrl`+`Enter`	run the cell and stay on it
`Alt`+`Enter`	run the cell and add a new one below
Run all	run every cell from the top, in order

Option	Means
`-q 20`	a base is good with Phred quality ≥ 20
`-u 40`	drop a read with more than 40 % bad bases (the default)
`-l 50`	drop reads shorter than 50 bases after trimming
`--cut_tail`	trim low-quality bases from the 3′ end in a sliding window
`--detect_adapter_for_pe`	find the adapters of paired-end reads from the reads themselves

Option	Means
`-a SEQ` / `-A SEQ`	a 3′ adapter of read 1 / read 2: the adapter and what follows it are removed
`-g SEQ` / `-g ^SEQ`	a 5′ adapter or primer (with `^`: only at the very start of the read)
`-q 20`	trim bases below quality 20 from the 3′ end (before the adapter is searched)
`-m 30`	discard reads (pairs) shorter than 30 bases after trimming
`-e 0.1`, `-O 3`	errors allowed in a match (a fraction of its length); the shortest match that counts
`--poly-a`, `--trim-n`	trim poly-A tails; trim N bases from the ends
`--discard-untrimmed`	keep only the reads in which the adapter was found (amplicons, primers)
`--json report.json`	the report as JSON too (for MultiQC, or `jq`)

Tab	What it is for
Notebook	Python and R notebooks (`.ipynb`) – Notebooks
Terminal	a bash shell with bioinformatics programs – The terminal
Files	your files: open, edit, upload, download – Your files
Galaxy	tools with forms, a history, workflows – Galaxy
Databases	UniProt, PDB, Ensembl, NCBI, BLAST … – Databases
3D viewer	protein structures – The 3D viewer
Genome browser	reads, variants, genes – The genome browser
Networks	interaction networks: style, analyse, export – The Networks tab
IDE	RStudio for R scripts and reports, Spyder for Python scripts and debugging – The IDE tab
AI assistant	chat and agent – The AI assistant

`Shift`+`Enter`	run a notebook cell and go to the next
`Ctrl`+`S`	save – a file you edit in the Files tab is saved only then (or with Save); notebooks also save themselves
`↑` `↓` in the terminal	earlier commands
`Tab` in the terminal	complete a file or command name
`Ctrl`+`C` in the terminal	stop: the program that is running finishes, the rest of the command is skipped; press it again to stop the program itself
`Ctrl`+`\`	hide or show the help

`↑` `↓`	earlier commands (they are kept when you reload the page)
`Tab`	complete a file name or a command
`Ctrl`+`C`	stop: the program that is running finishes and the rest of the command is skipped; press it again to stop the program itself (what it was writing is lost; `python` and `Rscript` start their language again, which forgets its variables). With text selected, it copies
`Ctrl`+`L`	clear the screen
`Home` / `End`	to the start / end of the line

For	Programs	Guide
reads	fastp, cutadapt, seqtk, wgsim	Reads
mapping	minimap2, bowtie2, bwa, bwa-mem2, hisat2 (RNA-seq), samtools, mosdepth	Mapping reads
RNA-seq counts	featureCounts, stringtie, kallisto, salmon	RNA-seq: counts and abundances
single-cell RNA-seq	bustools	Single-cell RNA-seq
ChIP-seq and ATAC-seq: peaks	macs3	Peaks
variants	bcftools, freebayes, bgzip, tabix	Variant calling
somatic variants (tumour and normal)	bam-somaticsniper, lofreq, configureStrelkaSomaticWorkflow.py	Somatic variants
metagenomics: who is in a sample	kraken2, kraken2-build, kraken2-inspect, bracken, bracken-build	Metagenomics
genotypes of many people: quality control, LD, association (GWAS)	plink2	Genotypes and association
genes and protein domains	prodigal; hmmbuild, hmmsearch, hmmscan, hmmpress, phmmer, jackhmmer, nhmmer, nhmmscan, hmmalign, hmmfetch, hmmstat, hmmemit, hmmconvert	Genes and protein domains
intervals and annotation	bedtools, gffread	Intervals
alignments and trees	needleman_wunsch, smith_waterman, kalign, muscle, fasttree, lastz	Alignments and trees
text and data	grep, sed, awk, sort, cut, jq, curl …	Text tools
languages	python, Rscript, snakemake, dot	Scripts and pipelines

Program	Your sequence	Searches
`blastp`	protein	proteins
`blastn`	DNA	DNA
`blastx`	DNA, translated in all six reading frames	proteins
`tblastn`	protein	DNA, translated in all six frames

When you…	…this is sent	to
open the workbench	requests for the page and its programs; and for the language of the notebook that opens (the last one you had open – Python on the first visit), which starts in the background	the web server that hosts the workbench (as for any web page); jsDelivr (Python) or webR's servers (R)
open a notebook of the other language, open RStudio or Spyder in the IDE tab, or install a package	requests for R and the package	webR's servers, the R package repositories (r-wasm.org, r-universe.dev, ropensci.org), PyPI
search in the Databases tab, or ask a service from code	the query you typed (for BLAST, also your e-mail address and the sequence)	the service: EMBL-EBI, NCBI, RCSB PDB, UCSC, STRING, Reactome …
make a network from STRING, IntAct or NDEx, or ask for the enrichment of genes, in the Networks tab	the names of the genes or proteins, or the search words – from a table of genes, the names of those that pass the filters, never the table's values	STRING (string-db.org), EMBL-EBI (IntAct), NDEx (ndexbio.org), g:Profiler (biit.cs.ut.ee)
load a structure, a hosted genome or a genome region	the identifier or region	RCSB PDB, AlphaFold DB; for hosted genomes the IGV team (igv.org) and UCSC (hgdownload.soe.ucsc.edu, or its copy on Amazon's servers)
ask the AI assistant	your question and what you tick under the box: the open notebook's code, your file names, your last commands, a summary of the open network (its size, columns, most connected nodes, the names of the selected ones and the results of analyses), your code in the IDE (from the IDE tab: the script open in RStudio or Spyder, the end of its console, and the names and types of its variables, but not their values). From a notebook cell: its code, with its error (up to 2,500 characters) or the beginning of its output (up to 1,200 characters); from a failed terminal command: the command, its folder and the end of what it printed (up to 3,000 characters); from an error in RStudio or Spyder: the code that failed (or the script that was run), the error message and the IDE’s working directory	the AI service you connected, under its terms
give the AI agent a task	the task, its commands and what they print (up to 3,000 characters a step) – including what it reads from your files	the AI service you connected

Component	Licence
Biowasm and Aioli (Robert Aboukhalil) – the WebAssembly builds of most programs, and the runner	MIT
fastp, seqtk, wgsim, minimap2, SAMtools, BCFtools, HTSlib, bedtools, LASTZ, jq	MIT (HTSlib also modified BSD)
gffread (with GCLib)	MIT (GCLib: Artistic 2.0)
Bowtie 2, Kalign	GNU GPL 3
BWA, HISAT2 – compiled for the browser by the genomics-wasm project	GNU GPL 3
featureCounts (Subread) – compiled for the browser by the genomics-wasm project	GNU GPL 3 or later
FreeBayes – compiled for the browser by the genomics-wasm project	MIT; it contains code under the GNU GPL 2 or later, so the program is distributed under the GPL
BWA-MEM2, StringTie, mosdepth (with Nim's runtime) – compiled for the browser by the genomics-wasm project	MIT (their parts: MIT, BSD, zlib – see each program's folder)
kallisto (with Bifrost) – compiled for the browser by the genomics-wasm project	BSD 2-Clause
Salmon (with its Rust crates and Rust's standard library) – compiled for the browser by the genomics-wasm project	BSD 3-Clause; the crates MIT, Apache 2.0 and other permissive licences – listed with the program
LoFreq, SomaticSniper (with SAMtools 0.1.6) – compiled for the browser by the genomics-wasm project	MIT (their parts: MIT, BSD, zlib and others – see each program's folder)
Strelka2 (with parts of Boost, HTSlib, CodeMin and RapidJSON) – compiled for the browser by the genomics-wasm project	GNU GPL 3 or later (its parts: Boost Software Licence, MIT and others – listed with the program)
Cutadapt, dnaio, xopen – built for the browser by the genomics-wasm project, with the Python they run in: Pyodide 314.0.7 (CPython 3.14.2)	MIT; Pyodide MPL 2.0 (Python: PSF licence; the libraries in it their own – listed with the program)
Kraken2 (Derrick Wood and others; with kseq and code from KrakenUniq) – compiled for the browser by the genomics-wasm project	MIT
Bracken (Jennifer Lu) – built for the browser by the genomics-wasm project; its Python programs run in the Pyodide 314.0.7 of Cutadapt	GNU GPL 3 or later
MACS3 (Tao Liu and the MACS3 project team; with fermi-lite, ksw, swalign and SIMDe) and hmmlearn – compiled for the browser by the genomics-wasm project; they run in the Pyodide 314.0.7 of Cutadapt, with Pyodide's builds of NumPy, SciPy, scikit-learn, pandas, h5py (and HDF5), joblib, threadpoolctl, python-dateutil, pytz, tzdata, six, packaging and pkgconfig, and with AnnData, array-api-compat and natsort	BSD 3-Clause (fermi-lite, ksw and swalign MIT; SIMDe MIT and CC0; the Python packages BSD, MIT and Apache 2.0 – listed with the program, with the libraries compiled into them: OpenBLAS, LAPACK, Boost, HDF5, zlib); Pyodide's recipes for its packages MPL 2.0, their source is with the program
bustools (BUStools: Páll Melsted, A. Sina Booeshaghi, Lior Pachter and others; with CRoaring, Eigen, kseq, CppNumericalSolvers, LBFGSpp and zlib) – compiled for the browser by the genomics-wasm project	BSD 2-Clause (its parts: CRoaring Apache 2.0; Eigen MPL 2.0, its source is with the program; the others MIT and zlib)
Prodigal (Doug Hyatt, University of Tennessee / UT-Battelle) – compiled for the browser by the genomics-wasm project	GNU GPL 3 or later
HMMER (Sean R. Eddy and the HMMER team; with Easel and libdivsufsort) – compiled for the browser by the genomics-wasm project, with a small change of the workbench (see its folder)	BSD 3-Clause (Easel: BSD 2-Clause; libdivsufsort: MIT)
PLINK 2 (Shaun Purcell, Christopher Chang, Benjamin Demaille and others; with libdeflate, zstd, SIMDe, SFMT, QD and zlib) – compiled for the browser by the genomics-wasm project	GNU GPL 3 or later (its library parts also LGPL 3 or later; libdeflate and SIMDe MIT, zstd BSD, SFMT and QD BSD-style, zlib its own)
FastTree	GNU GPL 2 or later
MUSCLE 3.8, seq-align	public domain
GNU coreutils, grep, sed, gawk, diffutils	GNU GPL 3 or later
Pyodide (CPython for the browser) and its packages	MPL 2.0 (Python: PSF licence; each package its own)
webR (R for the browser)	GNU GPL 3 (R: GNU GPL 2 and 3); its JavaScript: MIT
locfit (R package, built for this site with genomics-wasm's recipe); DESeq2, edgeR and the other R packages are downloaded when used	locfit: GNU GPL 2 or later; DESeq2: LGPL 3 or later; edgeR: GNU GPL 2 or later; each package its own
Graphviz, through Viz.js	EPL 1.0 (Viz.js: MIT)

Library	Licence
NGL Viewer (3D viewer)	MIT
igv.js (genome browser)	MIT
Cytoscape.js (Networks; the Cytoscape Consortium), with the layouts cytoscape-fcose, cose-base and layout-base (iVis@Bilkent) and cytoscape-dagre (with dagre and graphlib)	MIT
CodeMirror 5 (editors)	MIT
marked (Markdown)	MIT
DOMPurify (safe display of AI answers)	Apache 2.0 or MPL 2.0
IPython and jedi (Spyder’s console and code completion), from Pyodide’s packages	BSD 3-Clause (IPython), MIT (jedi)
pyflakes (Spyder’s code analysis)	MIT
Shinylive (Posit): Shiny apps in RStudio reach their R in the way that Shinylive’s do – written again for the workbench	MIT
JSZip (zip files)	MIT or GNU GPL 3
PDF.js (PDF files in the Files tab)	Apache 2.0