Instructions
Use the human or yeast dataset to investigate how the choice of aligner, variant caller or settings affects the results. Choose a comparison, run the pipelines and examine sites where they disagree.
Plan your comparison
Decide what you will change. Keep the reads, reference and other settings the same so that you can explain any differences. Consider where the tools might behave differently, such as repeats, indels or regions with low coverage.
Run the pipelines
Use the terminal, the Galaxy forms, or both. They use the same tools and files. Align the reads, sort and index the alignment, then call variants. Give each result a different file name so that you can compare them later.
Compare the results
Select the VCF files in the Compare tab. Check which alleles and genotypes differ. Two callers can make the same error, so agreement alone does not show that a call is correct.
Examine a difference
Inspect the alignment with samtools view, depth or mpileup. Use bcftools query to examine allele depths and quality scores. Are the reads supporting a variant on both strands? Could a repeat or an alignment difference explain the call?
Example comparisons and commands
The Example pipeline panel provides commands you can edit. You could compare minimap2 and Bowtie2 with the same bcftools calling model, or compare the multiallelic and consensus models using the same BAM file.
exactSNP uses a different calling algorithm. It reports SNPs and some simple indels but does not provide sample genotypes. Its SNP QUAL is min(40, −log10(p)), and simple indels have QUAL=1. These scores are not Phred-scaled, so choose its filters separately from bcftools filters.
Replace the file names and region below with those from your analysis.
samtools flagstat results/alignment.bam
samtools view results/alignment.bam contig:start-end
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\n' results/calls.vcf
awk '!/^#/ {print $1, $2, $4, $5}' results/calls.vcf
Available tools
The programs run in your browser on the files you supply. The Galaxy interface provides forms for running these tools; it is a teaching version of Galaxy.
Bowtie2 can take several minutes on a full dataset in this browser version. For an initial run, use minimap2 or a smaller paired subset of reads. Reducing the number of reads changes coverage and may affect the calls.
Run times depend on the browser build and your computer. They should not be used to rank the speed of the standard installed programs. The small reference sequences also leave out possible alignments elsewhere in the genome, which can affect mapping quality.
Create and edit text files
Type nano notes.txt to create a file, or use an existing file name to edit it. The editor opens inside the terminal. Selecting Text editor in Galaxy switches to the same terminal editor.
nano notes.txt cat notes.txt
Use Ctrl+O, then Enter, to save. Use Ctrl+X to exit. The shortcut labels at the bottom of the terminal are also clickable. If you exit with unsaved changes, press Y to save, N to discard, or Ctrl+C to keep editing.
The editor is useful for notes, small FASTA files, BED regions and lists of file names. The Tab key inserts a tab character for tab-separated files. It supports UTF-8 text up to 2 MiB; compressed and binary files cannot be edited.
When you exit, the command prompt and previous terminal output return. This browser editor uses nano-style controls. Type nano --help for its supported commands. Saved files are available to the terminal and Galaxy. Use Download work before reloading or closing the page to keep them.
Terminal commands
Type help for the available commands. The terminal supports pipes, input and output redirection, quoted arguments, simple filename patterns, variables and command chaining. You can use awk, grep, sed and other text utilities to process your files.
Shell loops, process substitution, background processes, package installation and arbitrary executables are not supported.
About the data
The datasets contain public sequencing reads from the NA12878 human exome and BY4741 yeast. The yeast reads were selected for mitochondrial alignment. The human FASTQs were recovered from a BAM file whose quality scores had already been recalibrated. No variants have been added. See Data sources for the accessions and selection methods.
The human examples use a diploid calling model. The yeast mitochondrial examples start with a haploid model. This assumes one allele at each site, but a sample can contain different mitochondrial sequences. Mitochondrial copy number is not the same as ploidy. Check possible low-frequency alleles for sequencing errors, ambiguous alignments and heteroplasmy.
The mitochondrial reference is linear, although the genome is circular. Reads crossing the origin may be clipped or fail to align.
VCF positions refer to the supplied FASTA sequences. The human reference consists of two genome regions with local coordinates; use the offsets in the data documentation to convert them to genome positions.
The comparison measures agreement between callers. It does not calculate accuracy, precision or sensitivity because no validated benchmark is included. The reads were selected using existing alignments, so reads that failed to map in those pipelines may be missing.
Your report
Describe your question, the commands and settings used, and the differences between the results. Examine the reads supporting at least one disagreement and explain any uncertainty. Use Download work to save your input and output files, notes and command history.