1 Overview of the tutorial

Using the data from the Nanopore sequencing you conducted,

Now the following sessions, we aim to:

Compare the read length and quality between the experimental conditions
Investigate how the read cleaning process affect the read yields
Identify the genomic variants from the sequence data
Interpret how the genomic variants affects the animal biologically

And we will learn:

How to use Orion and conduct genome analysis
Quality check, read filtering, mapping to the reference genome and variant calling
How to interpret summary statistics of Nanopore sequence data
How to identify and interpret genetic variants

Using Orion

Prepare your computational tools with Conda

2 Reads quality check

Overview:

Quality check -> Trimming of low quality reads -> Quality check
Compare the overall reads quality between conditions

2.1 Connect to Orion and the prepare the tools

Go to https://orion.nmbu.no/ at NMBU or with VPN.

In the Terminal/Command prompt, go to your directory.

Review: the concept of current directry

cd your_directory

Let’s make a directory for analysis and enter in it.

mkdir analysis # make directory "analysis"
cd analysis # set the current directory "analysis"

Now, you will inspect the fastq file from your experiment, which contains Nanopore read information.

2.2 Check the read quality by Nanoplot

2.2.1 Browse the inside of the read (fastq) file

Review: look into a file content in a command line

Let’s learn how a fastq file (sequencing reads) looks using a sample file.

zcat /mnt/courses/BIO326/EUK/result_2024/pig_demo_fasta_BC01.fastq.gz | more

Now you are seeing the content of a fastq file. (.gz = compressed)

Each entry in a FASTQ files consists of 4 lines:

A sequence identifier with information about the sequencing run. (run time, run ID, cflow cell id … )
The sequence (the base calls; A, C, T, G and N).
A separator, which is simply a plus (+) sign.
The base call quality scores. These are Phred +33 encoded, using ASCII characters to represent the numerical quality scores.” quality score sheet

by Illumina Knowledge

2.2.2 Get basic stats of the fastq file

“zcat”-> look inside

“wc” -> word count

“-l” -> line

zcat /mnt/courses/BIO326/EUK/result_2024/pig_demo_fasta_BC01.fastq.gz | wc -l

2.2.3 Discussion Point

Now you got the number of lines in the fastq file.

How many sequence reads are in the fastq file?

What is the quality of bp between 1 and 5? What is the quality of bp between 21 and 25? Why do you think they are different?

Need Help?

We see that there are 48000 lines in the fastq file.

As we learned that “each entry in a FASTQ files consists of 4 lines”, one read is corresponding to four lines. So in this file we have 48000/4 = 12000 reads.

Hint for the 5 first bp:

zcat /mnt/courses/BIO326/EUK/result_2024/pig_demo_fasta_BC01.fastq.gz | sed -n 4p | cut -c 1-5

2.3 Quality check by Nanoplot

The original fastq files may contain low quality reads. In this step, we will use “Nanoplot” to see the quality and lentgh of each read.

Make a slurm script to conduct the quality check on the sample file like below and run it.

Review: make a slurm script and run it by sbatch


#!/bin/bash
#SBATCH --job-name=Nanoplot  # sensible name for the job
#SBATCH --mail-user=yourname@nmbu.no # Email me when job is done.
#SBATCH --mem=12G 
#SBATCH --ntasks=1   
#SBATCH --cpus-per-task=8
#SBATCH --mail-type=END

##Activate conda environment
module load Miniconda3 && eval "$(conda shell.bash hook)"


conda activate /mnt/courses/BIO326/EUK/condaenvironments/condaEUK

echo "Working with this $CONDA_PREFIX environment ..."


## run nanoplot

NanoPlot -t 8  --fastq /mnt/courses/BIO326/EUK/result_2024/pig_demo_fasta_BC01.fastq.gz --plots dot  --no_supplementary --no_static --N50 -p before_BC01

Nanoplot will generate the result files, named “before”xxx. Lets look into them…

Review: File transfer between Orion and your computer


# taking too long?
qlogin 

cp /mnt/courses/BIO326/EUK/result_2024/before_BC01NanoPlot-report.html  before_BC01NanoPlot-report.html

Open “before_BC01NanoPlot-report.html” on your local computer

##Filtering by Nanofilt


#!/bin/bash
#SBATCH --job-name=Nanofilt  # sensible name for the job
#SBATCH --mail-user=yourname@nmbu.no # Email me when job is done.
#SBATCH --mem=12G 
#SBATCH --ntasks=1   
#SBATCH --mail-type=END

##Activate conda environment
module load Miniconda3 && eval "$(conda shell.bash hook)"

conda activate /mnt/courses/BIO326/EUK/condaenvironments/condaEUK

echo "Working with this $CONDA_PREFIX environment ..."

## run nanofilt
gunzip -c /mnt/courses/BIO326/EUK/result_2024/pig_demo_fasta_BC01.fastq.gz | NanoFilt -q 10 -l 500 | gzip > cleaned_pig_demodata_BCO1.fastq.gz

-l, Filter on a minimum read length

-q, Filter on a minimum average read quality score

In this case, we are removing reads lower than quality score 10 and shorter than 500 bases.

If you are ambitious, please adjust the filtering parameters and see how they change the result.

(In that case, do not forget to name the result files differently.)

2.4 Compare the sequences before and after cleaning

Run Nanoplot again on the cleaned sequences.

Need help?


#!/bin/bash
#SBATCH --job-name=Nanoplot  # sensible name for the job
#SBATCH --mail-user=yourname@nmbu.no # Email me when job is done.
#SBATCH --mem=12G 
#SBATCH --ntasks=1   
#SBATCH --cpus-per-task=8
#SBATCH --mail-type=END

##Activate conda environment
module load Miniconda3 && eval "$(conda shell.bash hook)"


conda activate /mnt/courses/BIO326/EUK/condaenvironments/condaEUK

echo "Working with this $CONDA_PREFIX environment ..."


## run nanoplot
NanoPlot -t 8  --fastq /mnt/courses/BIO326/EUK/result_2024/cleaned_pig_demodata_BCO1.fastq.gz --plots dot  --no_supplementary --no_static --N50 -p after_BC01

Open “after_BC01NanoPlot-report.html” on your local computer.


# taking too long?
qlogin 

cp /mnt/courses/BIO326/EUK/result_2024/after_BC01NanoPlot-report.html after_BC01NanoPlot-report.html

2.4.1 Discussion Point

Did you see the difference of read and quality distribution between before and after the filtering?

2.5 Your mission

Do the quality check and filtering, and compare the read length and quality between the four experimental conditions.

3 Mapping to the reference genome

3.1 run Minimap and map the reads to the reference genome

cow genome

#!/bin/bash
#SBATCH --job-name=Sniffles  # sensible name for the job
#SBATCH --mem=18G 
#SBATCH --ntasks=1   
#SBATCH --cpus-per-task=8

##Activate conda environment
module load Miniconda3 && eval "$(conda shell.bash hook)"


conda activate /mnt/courses/BIO326/EUK/condaenvironments/condaEUK

echo "Working with this $CONDA_PREFIX environment ..."


minimap2 -t 8 -a -a /mnt/courses/BIO326/EUK/cow_genome/Bos_taurus.ARS-UCD1.3.dna.toplevel.fa.gz  /mnt/courses/BIO326/EUK/result_2024/cleaned_cow_demo_data.fastq.gz > cow.sam


# convert the sam file to bam format
samtools view -S -b cow.sam > cow0.bam

## sort the bam file
samtools sort cow0.bam -o cow.bam

# index the bam file
samtools index -M  cow.bam

# Variant Calling using Sniffles
sniffles --input  cow.bam --vcf cow.vcf

bam and sam format

Now you got the variant file!

4 Investigate the variants

Celian will explain how to read a vcf file
Go to mentimeter: https://www.menti.com/alvrtwtcrkj5 (or use code 6182 8854)

Copy the vcf in your directory

qlogin 

cp /mnt/courses/BIO326/EUK/result_2024/cow_Control.vcf  /cow_Control.vcf

Look inside the vcf

# INFO field

grep '^##' cow_Control.vcf | tail -n 20

# variants
grep -v '^##' cow_Control.vcf | more

Important parameters

1 16849578 : location of the variant

SVTYPE=DEL;SVLEN=-60 : size and type of the variant

0/1 : genotype

(you can open a vcf file in notepad, excel etc.)

Now you have variants! Lets see which genes are affected by the variants.

For that you will select a random variant to investigate.

#Check the number of variant in the file

NBVAR=$(grep -v '^##' cow_Control.vcf | wc -l)

## sample a random number

RANDOMVAR=$(echo $((RANDOM % $NBVAR + 1)))

## let's check the variant sampled

grep -v '^##' cow_Control.vcf |  sed -n ${RANDOMVAR}p

4.1 Estimate the effect of variants

Go to VEP (Variant Effect Predictor)

Variant Effect predictor tells us where in the genome the discovered variants are located (genic, regulatory, etc…)

Select “cow” as the reference species.

Upload: cow_Control.vcf - downloaded from Orion or the section above as the file to investigate.

There are 1481 variants; 344 genes are affected by these varaints.

Let’s closely investigate your variant !

Find your variant by downloading the .txt file

hint: Open the file and use Ctrl+F and paste the ID of your variant, which can be find in the vcf file

If you don’t find your variant, what does that mean ? Try another one !

What does your variant look like in the sequence ?

We can look at it on IGV

Select the cow reference genome - ARS-UCD1.2/bosTau9

Then import the vcf from your computer

You can select chr1 and put the location of your variant to look at it !

hint: remember the location of your variant is in the vcf file

That’s it ! You are now expert in genome analysis ! If you’re considering undertaking a similar bioinformatics/genomics project for your master’s thesis, you can contact our lab

Bio326_2024

MS and colleagues

2024-01-29