Skip to content

Latest commit

 

History

History
270 lines (200 loc) · 9.75 KB

lab-5.md

File metadata and controls

270 lines (200 loc) · 9.75 KB
Raw
tags
ggg, ggg2022, ggg201b

Lab 5 outline - variant calling odds and ends; transition to assembly (and annotation). - GGG 201(b) lab, Feb 4, 2022

hackmd-github-sync-badge

(Permanent link on github)

Contents:

[toc]

preparing for class

Before or at the very beginning of class, please log into farm and run:

mamba create -n assembly -y megahit prokka quast
  • thank you!

revisiting the homework

Re homework #1, let's just revisit the steps:

  • logged into farm
  • git clone: made a copy of a set of files stored in your GitHub account (your "repository")
  • srun - allocated a specific computer from the larger farm cluster
  • ran and edited and tested your workflow
  • git commit - save changes locally
  • git push - send saved changes back to your github repo. <-- this is what hands in your assignment

I have access your github repository and will use that to look at your handed-in assignment.

variant calling vs de novo assembly

revisiting variant calling: assumptions and scope:

  • which do we want? SNPs, short indels, long indels/structural variation
  • accurate short reads work well for SNPs
  • inaccurate long reads vs accurate long reads
  • assumes that we have a reference that includes the main things
  • value: VEP, population structure, GWAS

but variant calling is almost always reference based. It's even in the name - "variant"!

Where do our reference genomes come from?

De novo assembly.

(Diagram stuff ensues HERE.)

running our first assembly

We're going to assemble a paired-end E. coli sample from the Long Term Evolution Experiment - SRA run SRR2584857. It's too big for us to download quickly, so I've provided it locally.

get the data

mkdir ~/ggg201-week5
cd ~/ggg201-week5

cp ~ctbrown/data/ggg201b/SRR2584857_*.fastq.gz .
ls -lh

and you should see something like

>total 360M
>-rw-rw-r-- 1 datalab-01 datalab-01 180M Feb  4 07:00 SRR2584857_1.fastq.gz
>-rw-rw-r-- 1 datalab-01 datalab-01 180M Feb  4 07:00 SRR2584857_2.fastq.gz

allocate a node

srun --nodes=1 -p high2 -t 2:00:00 -c 4 --mem 6GB --pty /bin/bash

here we're asking for 4 CPUs, 6 GB of RAM, for 2 hours.

using a conda environment w/software

we're going to use the following software:

  • megahit, a de novo genome assembler for microbial samples; alternatives here would be SPAdes.
  • quast for genome assembly evaluation.
  • prokka for bacterial and archaeal genome annotation.

In brief, megahit takes a pile of reads and gives you a genome; quast gives you statistics about that genome; and prokka annotates the genome with genes.

At the beginning of class I asked you to install all of these; you don't need to run this again if you already have:

mamba create -n assembly -y megahit prokka quast

and now let's use that conda environment:

conda activate assembly

running our first assembly

megahit -1 SRR2584857_1.fastq.gz -2 SRR2584857_2.fastq.gz -f -m 5e9 -t 4 -o SRR2584857_assembly

this will take about 3 minutes, and produce a directory, SRR2584857_assembly/.

let's look at the contents of that directory:

% ls -l SRR2584857_assembly
total 3340
-rw-rw-r-- 1 datalab-01 datalab-01     230 Feb  4 07:28 checkpoints.txt
-rw-rw-r-- 1 datalab-01 datalab-01       0 Feb  4 07:28 done
-rw-rw-r-- 1 datalab-01 datalab-01 4562183 Feb  4 07:28 final.contigs.fa
drwxrwxr-x 2 datalab-01 datalab-01      70 Feb  4 07:28 intermediate_contigs
-rw-rw-r-- 1 datalab-01 datalab-01  120250 Feb  4 07:28 log
-rw-rw-r-- 1 datalab-01 datalab-01     948 Feb  4 07:25 options.json

The key file here is final.contigs.fa; the rest are log files or temporary files.

Let's look at final.contigs.fa:

head SRR2584857_assembly/final.contigs.fa

...well, that's a lot of DNA sequences...

How do we evaluate assemblies?

Let's start by getting some high level statistics. That's what quast is for.

Run quast

Let's copy the assembled file out of the subdirectory to make it more convenient to work with...

cp SRR2584857_assembly/final.contigs.fa SRR2584857-assembly.fa

and then run quast:

quast SRR2584857-assembly.fa

This will produce a directory quast_results/ with a lot of files in it. Let's look at the basic text report:

cat quast_results/latest/report.txt

assembly metrics and stats explained

You should see something like:

All statistics are based on contigs of size >= 500 bp, unless otherwise noted (e.g., "# contigs (>= 0 bp)" and "Total length (>= 0 bp)" include all contigs).

Assembly                    SRR2584857_assembly
# contigs (>= 0 bp)         126                
# contigs (>= 1000 bp)      92                 
# contigs (>= 5000 bp)      68                 
# contigs (>= 10000 bp)     65                 
# contigs (>= 25000 bp)     52                 
# contigs (>= 50000 bp)     33                 
Total length (>= 0 bp)      4557069            
Total length (>= 1000 bp)   4542425            
Total length (>= 5000 bp)   4474835            
Total length (>= 10000 bp)  4451946            
Total length (>= 25000 bp)  4249832            
Total length (>= 50000 bp)  3540628            
# contigs                   103                
Largest contig              326752             
Total length                4549884            
GC (%)                      50.72              
N50                         98799              
N75                         52265              
L50                         16                 
L75                         31                 
# N's per 100 kbp           0.00

A few top-level observations:

  • assemblies are often fragmented, especially if they use short-read sequencing. This is because of repeats and/or low coverage.
  • you may have an estimate of the total length of your genome from other sources (experimental biology, for example :). Typically your assembly will be either a bit SHORTER than that, for haploid samples (because of repeat collapse/loss) or MUCH MUCH LONGER (because of diploidy or polyploidy).

Given that, most of these numbers should be self-explanatory, except for the N50, N75, L50, and L75.

To calculate these, you need to rank-order the contigs by size (largest to smallest, say) and then pick two ranks:

  • the rank at which the sum of the length of all contigs below that rank equals 50% of the length of the whole assembly
  • the same, but for 75%

Then:

L50: the number of contigs below the 50% rank N50: the contig length at the 50% rank L75: the number of contigs below the 75% rank N75: the contig length at the 75% rank

These exist because assembly is always interested in producing long contigs from shorter ones. Again, that's in the name :). So you want measures that reflect how much of the content of your assembly is in long bits.

Let's consider some extremes:

  • if your entire genome is assembled into one contig, then your L50 is 1, and your N50 is the genome length.
  • if your reads completely failed to assemble, then your L50 is the number reads, and your N50 is the size of an individual read.

annotating the assembly

Right now the assembly is just a pile of FASTA sequences. That's not so useful. Let's run the prokka bacterial / archaeal genome annotator:

prokka --prefix SRR2584857_annot SRR2584857-assembly.fa

This will produce a directory SRR2584857_annot; let's take a look:

ls -lh SRR2584857_annot

you will see:

-rw-rw-r-- 1 datalab-01 datalab-01 1.3M Feb  4 07:44 SRR2584857_annot.err
-rw-rw-r-- 1 datalab-01 datalab-01 1.5M Feb  4 07:44 SRR2584857_annot.faa
-rw-rw-r-- 1 datalab-01 datalab-01 4.1M Feb  4 07:44 SRR2584857_annot.ffn
-rw-rw-r-- 1 datalab-01 datalab-01 4.5M Feb  4 07:41 SRR2584857_annot.fna
-rw-rw-r-- 1 datalab-01 datalab-01 4.5M Feb  4 07:44 SRR2584857_annot.fsa
-rw-rw-r-- 1 datalab-01 datalab-01 9.3M Feb  4 07:44 SRR2584857_annot.gbk
-rw-rw-r-- 1 datalab-01 datalab-01 5.5M Feb  4 07:44 SRR2584857_annot.gff
-rw-rw-r-- 1 datalab-01 datalab-01  62K Feb  4 07:44 SRR2584857_annot.log
-rw-rw-r-- 1 datalab-01 datalab-01  15M Feb  4 07:44 SRR2584857_annot.sqn
-rw-rw-r-- 1 datalab-01 datalab-01 900K Feb  4 07:44 SRR2584857_annot.tbl
-rw-rw-r-- 1 datalab-01 datalab-01 289K Feb  4 07:44 SRR2584857_annot.tsv
-rw-rw-r-- 1 datalab-01 datalab-01  113 Feb  4 07:44 SRR2584857_annot.txt

and check out the summary information:

cat SRR2584857_annot/*.txt

You should see something like:

>organism: Genus species strain
>contigs: 126
>bases: 4557069
>CDS: 4205
>rRNA: 8
>repeat_region: 2
>tRNA: 76
>tmRNA: 1

so that's a summary of what prokka found. Note that repeat regions in particul are going to be underestimated in number.

What are all the other files?

Well, SRR2584857_annot/SRR2584857_annot.log contains a copy of the text output from prokka; you can page through it with more SRR2584857_annot/SRR2584857_annot.log (use 'q' to quit).

As for the rest, there are a variety of file formats in which annotations are stored. The ones I have found most useful are:

  • SRR2584857_annot/SRR2584857_annot.faa - contains protein FASTA sequences for the predicted genes
  • SRR2584857_annot/SRR2584857_annot.ffn - contains DNA sequences for the predicted genes

although there's a bunch of other potentially useful output in there if you plan on submitting a bacterial genome to Genbank :).

You can read more about prokka here.

next class

we'll build a bit of a snakefile.

we'll talk about parameter exploration.

and we'll discuss how assembly works underneath.