Bioinformatics Bootcamp 2024

Schedule

• a computer with internet connection
• an account on super-computing Wales
• basic knowledge of DNA and genome structure

Lecture slides

https://drabarlow.github.io/bioinformatics_bootcamp/

Bootcamp worksheet

https://drabarlow.github.io/bioinformatics_bootcamp/bootcamp_worksheet.html

Github repo

https://github.com/drabarlow/bioinformatics_bootcamp

Supercomputing Wales shared project

• please request access to: scw2141 Bioinformatics training

Axel Barlow, Lecturer in Zoology at Bangor University

Interests

• Population and evolutionary genomics of animals. Palaeogenomics of extinct animals. Conservation genomics of native species

Bioinformatics experience

• Working in the field for ten years
• Entirely self taught post-PhD
• No proper bioinformatics/computer science background
• Knowledge of bash and R.

• £16m investment, part-funded by the European Regional Development Fund (ERDF) through Welsh Government
• Provide university research teams access to HPC
• Consortium of Cardiff, Swansea, Bangor and Aberystwyth
• 2 Supercomputers:
  • Cardiff HPC System - Hawk
  • Swansea HPC System - Sunbird
• Hawk: 280 nodes, 12,736 cores, 68.224 TB memory

• Scientific linux OS
• compute nodes
  • 134x Intel nodes with 20 cpus + 192 Gb RAM each
  • 64x AMD nodes with 32 cpus + 256 Gb RAM each
• high memory nodes
  • 26 nodes with 384 Gb RAM each
• GPU nodes
  • 26 Nvidia P100 GPUs on 13 nodes
  • 30 Nvidia V100 GPUs on 15 nodes
• Storage
  • 1192TB (usable) scratch space
  • 420TB of home directory space

• / [root] is uppermost level of filesystem
• Everything is contained in /
• Directories exist within the filesystem, they can contain files and other directories
• We specify a path through this hierarchy using forward-slashes (note back-slashes on DOS)
• Our current directory is called the working directory

/home/b.xlb21brx/
/scratch/b.xlb21brx/

• We can navigate through the filesystem (change working directory)
• Or we can specify the patch to directories or files remotely

1. Sample preparation
2. Bind DNA to flowcell, generate clusters
3. Sequencing by synthesis
4. Data analysis (in the machine)

https://www.youtube.com/watch?v=fCd6B5HRaZ8

*Indexes allow multiple samples to be sequenced at the same time

[Not an exhaustive list]

• Whole genome sequencing (pure DNA sample from a single individual)
• Reduced representation genome data (RADseq, targeted SNPs, single individual)
• Poolseq (multiple individuals)
• Transcriptome (RNA sample from single tissue/individual)
• Metabarcoding (PCR amplicon, multiple individuals/species)
• Metagenomics (whole genomes, multiple individuals/species)

Short reads from a single individual can be mapped to a reference genome assembly

• The current market leader
• Massive output
• Many applications (genome resequencing, RADseq, transcriptomes, metabarcoding)
• Cheap (£10 per Gb)
• Major limitation is the read length

• Nikolas Basler, Achim Klittich, Axel Barlow
• An environment for organising, processing, and archiving Illumina data
• BEARCAVE philosophy
  • All users can access all data
  • Avoid data redundancy
  • All samples processed using identical software programs and parameters
  • Incorporates sample metadata
  • Documents results of data processing
  • Easy to use wrapper scripts for programs
  • Publicly available
  • Safeguards in place to ensure consistency
• Consequently BEARCAVE is not for everyone, and has idiosyncrasies in use

• fastq is the standard output format for data from Illumina (and other) platforms

@A00551:758:HKTVJDSX7:4:1101:3595:6872 1:N:0:CCTGAGATGT+GGTCTAGTTG
CTGAATATGGATTTTAATTGAATCCTAAGATATTATAGCATCTTTCACTCCCTGTCCTGTGCATGTCAGA
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

• Line 1: info on sequencer, flowcell, cluster position, indexes (sometimes)
• Line 2: called bases
• Line 3: a +
• Line 4: quality scores on Phred scale
• 10 = 90% accuracy; 20 = 99% accuracy; 30 = 99.9% accuracy
• Recoded as single character: F = 37; ? = 30; 5 = 20; + = 10

• DNA can be fragmented
• The fragment lengths have a distribution

• decompress fastqs
• trim adapter seqs using Cutadapt
  • 30 bp min length
  • min overlap 1 bp
• merge overlapping read pair using FLASH
• recompress files and clean up
• save appropriate log files

Expected output in /BEARCAVE2/trimdata/*processing/

• merged reads *_mappable.fastq.gz [big file]
• unmerged R1 *_mappable_R1.fastq.gz [big file]
• unmerged R2 *_mappable_R2.fastq.gz [big file]
• trim report *_trim_report.log and merge report *_merge_report.log

• decompress fastqs
• merged (SE) and PE data processed separately
  • mapping using bwa mem algorithm
  • PCR duplicates identified and removed using samtools
  • Reads with poor mapping quality (Q30) removed using samtools
• SE and PE data merged
• mapping log file generated
• file cleanup and renaming

Expected output in /BEARCAVE2/mapped*/*processing/

• mapped filtered data *.bam [big file]
• bam index *.bam.bai
• mapping log *_mapping.log

• Widely used program
• MANY population genetics analyses possible
• Tends to work directly from bam files (unlike plink, admixtools, etc)
• MANY filters available
• Genotype likelihood approach is a particular speciality

• Several spin off programs that use GLs
  • NGSadmix
  • PCangsd
  • NGSrelate
  • realSFS

Covariance matrix

• All indviduals
• Allele frequency covariance matrix
• Used for PCA

Distance matrix

• All individuals
• absolute genetic distance between populations
• used as input for NJ algorithm

Heterozygosity

• Your single adder
• Calculate GLs, ML estimation of SFS along a sliding window using realSFS

Suppose you're a survey company. To carry out your survey you need all the people seated in a classroom, which you have to build. You're not sure how many, so you build an ordinary classroom, with 5 rows of 6 desks for 30 people, after 30 people file in you notice there's a 31st. You build a second 30-person classroom right next to the first, and now you can accept 60 people, but then you notice a 61st. So you ask them to wait, and you build two more classrooms, so now you've got a nice 2x2 grid of 30-person classrooms, but the people keep coming and soon enough the 121st person shows up and there's not enough room. So you build a big 5-story building next door with 50-person classrooms, 5 on each floor, for a total of 50 x 5 x 5 = 1,250 desks, and you have the first 120 people file out of the old rooms into the new building, and you hire some wreckers to demolish the old classrooms and recycle some of the materials, and the people keep coming. And when you're all done with all this, the only "survey question" you're going to ask is "How many rows are there?"

Meanwhile, Bob's discount survey company, who can only tell you how many people he surveyed, is down there on the streetcorner, and the people are filing by, and Bob is jotting down tally marks on his clipboards, and the people, once surveyed, are walking away and going about their business, and Bob isn't wasting time and money building any classrooms at all.

An abridged version of https://stackoverflow.com/questions/30948366/why-is-unix-terminal-faster-than-r

Most people disagree (in some cases strongly)

• Rstudio is terrible (except for R markdown)
• Base R is really good
• ggplot2 code is hellishly complex
• tidyverse is not the way to teach R to beginners
• It's boring if everyone's graphs look the same
• It's OK to type commands directly into the terminal
• It's OK to modify a raw data file
• It's OK to combine other tools and languages
• Data doesn't have to be tidy

Objects

• Contain data and results
• Created with <-
• Stored in active memory
• Names include letters numbers _ .
• Names must begin with a letter

Functions

• Carry out operations on objects
• Often generates new objects
• function()
• ?function

Vector

• List of values of the same type
• Numbers, strings, or logical values
• Can be generated using c()
• Indexing vector objects my_vector[]

Matrix

• 2D data of same type in rows and columns
• Indexing matrix objects my_matrix[row, column]

Dataframe

• Rectangular "table"
• Mixture of data types
• Set of vectors of equal length
• Extract columns with $, which can then be indexed like vectors

List

• Set of components with different structures
• Extract named components with $

• Input is sliding window estimates of heterozygous sites