Lecture 4

Copy files for part 1 from `/home/hgen_share/lec4_pt1`

Set up directory and copy over files

module load StdEnv/2020 gcc/9.3.0 openmpi/4.0.3 hdf5/1.10.6 minimap2/2.17 fastqc/0.11.9 samtools/1.11 mugqic/qualimap/2.2.1 mafft/7.471 emboss/6.6.0 fasttree/2.1.11 mugqic/pycoQC/2.5.2
export PATH="/home/hgen_share/Anaconda/bin:$PATH"
export PATH="/home/hgen_share/Anaconda/envs/longread_env/bin:$PATH"
export PATH=/home/hgen_share/ont-guppy-cpu/bin:$PATH

Take a look at `ex.fast5`

View the file with h5ls

h5ls ex.fast5 | less

Show structure of one read using h5ls

h5ls -r ex.fast5/read_00017b75-18fc-40ae-8a21-a5b3e49dc753

Basecall using Guppy

This will take a while with just 1 CPU thread, and since the output is provided already, you can terminate it prematurely with Ctrl-C

conf="dna_r9.4.1_450bps_hac"
guppy_basecaller \
  --compress_fastq \
  -i ex.fast5 \
  -s basecall/ \
  --cpu_threads_per_caller 1 \
  --num_callers 1 \
  -c dna_r9.4.1_450bps_hac.cfg

(Extra) Assess quality with pycoQC & FastQC

Run pycoQC with the provided summary

pycoQC -f sequencing_summary.txt -o pyco.html

Run FastQC with the provided sequences

fastqc ex.fq

Download pycoQC html report

Compare with another report here

(Extra) Rectify quality issues by trimming

Trim adapters with Porechop

porechop -i ex.fq -t 1 -v 2 -o chopped.fq

Is it better now? Run fastQC on trimmed fastq

(Extra) Assess error profile after alignment

Align trimmed reads using minimap2

minimap2 -t 1 -x map-ont -a wuhCor1.fa chopped.fq | samtools sort -o ex.bam
samtools index ex.bam

Summarize alignment results with Qualimap

qualimap bamqc -bam ex.bam --java-mem-size=1G -nw 5000 -nt 1 -c -outdir qualimap/

Download results to your local computer

Compare with the report from a typical Illumina run

Hybrid assembly with Unicycler

Set up variables and copy over files for part 2

export PATH=/home/hgen_share/utils/SPAdes-3.15.5-Linux/bin:$PATH
export PATH=/home/hgen_share/utils/ncbi-blast-2.13.0+/bin:$PATH
export PATH=/home/hgen_share/utils/racon/build/bin:$PATH
export PATH="/home/hgen_share/Anaconda/bin:$PATH"
export PATH="/home/hgen_share/Anaconda/envs/longread_env/bin:$PATH"
cp /home/hgen_share/lec4_pt2/* ./

Run Unicycler using boths short- and long-reads. For hybrid assemblies, the direct long-read bridging step of the pipeline can take a while to complete. Thus, the output is provided for you. You can terminate it prematurely with Ctrl-C

unicycler -1 Illumina_F.fq -2 Illumina_R.fq -l ONT.fq -t 20 --min_fasta_length 10000 --linear_seqs 1 -o unicycler

The output can be found in /home/hgen_share/lec4_pt2/unicycler

Evaluate the assembled genome's quality

Use QUAST to create a summary report

module load mugqic/Quast/5.0.2
quast.py -t 1 -o quast -R wuhCor1.fa -g wuhCor1.gff3 unicycler/assembly.fasta

Compare with a report on eukaryotic assemblies

BLAST the top 2 contigs (contig1.fa and contig2.fa) found in the unicycler output directory

Place newly assembled viral genomes in a global context

Set up variables for nextstrain environment

export PATH="/home/hgen_share/Anaconda/bin:$PATH"
export PATH="/home/hgen_share/Anaconda/envs/nextstrain_V2/bin:$PATH"

Copy the tutorial dataset for SARS-Cov2 provided by nextstrain

cp -r /home/hgen_share/ncov ./

Take a look at the /data sub-directory. Inspect the data files with less to get an impression of what data needs to be provided for nextstrain. Unzip the fasta file using gunzip

Go back to the main /ncov directory and run the basic nextstrain workflow

snakemake --cores 14 --profile ./my_profiles/getting_started

On your own computer, not the cluster, establish a tunnel to server. Here you should use your own username.Also, replace the hashtag with any two numbers to land on a port that is not busy

ssh -L 40##:localhost:40## studXX@workshop2021a.vhost37.genap.ca

Start the auspice server to view the results using your credentials. If the PORT is busy, change to another port by exporting the port number as indicated below

export PORT=40##
auspice view --datasetDir /home/studXX/projects/lec4/ncov/auspice

You can visit the results page at http://localhost:40## in your browser. If you're done inspecting the results, terminate the auspice server with Ctrl+C

Include our own recently assembled data. Use consensus_sequences.fasta and metadata.tsv from /home/hgen_share/lec4_pt2. Append the fasta file to the example sequences (you can also inspect the fasta file with less)

cat consensus_sequences.fasta >> data/example_sequences.fasta

Do the same for the metadata file. Append it to example_metadata.tsv

Run the same steps as above. Where does our data reside in the phylogeny?