Exercises series 1

Before diving into the exercises, please read the General Guidelines

In this first series of exercises, we will implement the basic version of our workflow:

• align reads to a reference genome
• sort and index the resulting alignment files
• call genomic variants from the sorted alignment files
• process the variant calling results to generate a table of substitution counts

---

Exercise 1

Context

We will start by aligning reads to the reference assembly using the BWA mem algorithm, which is arguably the standard aligner for short reads.

Task

Implement a rule to align reads to the already indexed assembly using BWA.

Details

• The reference assembly file is data/genome.fa
• For now, your rule should only align the existing reads file: data/samples/A.fastq
• The command to align reads to an indexed assembly with BWA is: bwa mem <assembly_file> <reads_file>
• BWA mem outputs alignments in the SAM format; usually, it is advised to use the BAM format (a compressed binary version of SAM). SAM files can be converted to BAM with the software Samtools; in this case, we can pipe the output of BWA to samtools to obtain a BAM file with the following command: bwa mem <assembly_file> <reads_file> | samtools view -b > <output_file>. By default, the results of samtools view are output to stdout, we need to redirect it to the actual output file with >
• The output file will have the extension .bam (don't forget to specify this in the output directive!).

Relevant Wiki section

Defining rules

---

Exercise 2

Context

We have implemented a rule to align reads to a reference, but this rule only works for the reads file data/samples/A.fastq. We would like to align another reads file data/samples/B.fastq, and maybe even more reads files in the future.

Task

Generalize the alignment rule to handle any reads file present in data. From now on, rules should always be written to handle any relevant input from previous rules.

Relevant Wiki section

Wildcards

Hint

Click to expand

Use a wildcard to specify the reads file name.
Don't forget the output file extension: '.bam'

---

Exercise 3

Context

We aligned the reads to the reference and generated BAM files. Many analyses require that BAM files are sorted by genomic coordinates, because sorted BAM files can be processed much more easily than unsorted ones.

Task

Implement a rule to sort the BAM files with Samtools.

Details

• The command to sort BAM files with Samtools is samtools sort -O bam <alignment_file> > <output_file>
• The -O parameter specify the output format. By default, the results are output to stdout, which we redirect to the output file with >.
• The output file will have the extension .sorted.bam (to differentiate them from the alignment bam files).

---

Exercise 4

Context

Just like we did with the reference assembly, it is common to index alignment files to speed up access to specific regions in the alignment.

Task

Implement a rule to index the sorted BAM files with Samtools.

Details

• The command to index BAM files with Samtools is samtools index <alignment_file>
• For this rule, you don't have to specify the output file in the command, as it will be automatically generated. You'll still need to specify the output in the output directive! The index file will have the extension .bam.bai.

---

Exercise 5

Context

We aligned reads from two samples (A and B) to the reference assembly, and obtained two sorted and indexed alignment files. We will now call variants from these two files together; the output will be a VCF file.

Task

Implement a rule to call genomic variants from the alignment files using BCFtools.

Details

• The command to generate the VCF file is: bcftools mpileup -f <assembly_file> <alignment_files> | bcftools call -mv - > <output_file>
• Variant calling processes all alignment files together; use concepts from this section to do that. Your rule should also make sure that the alignment files are indexed
• The -v flag in bcftools call is used to output only variant sites. The -m flags specify to use a specific variant calling model (multiallelic and rare variants caller)
• The output file will have the extension .vcf

Relevant Wiki section

The expand syntax

Hints

Process all alignment files

- Use the expand syntax to collect all alignment files.
- Add a list of samples to your Snakefile for convenience.

Ensure alignment files are indexed

Add the output of the indexing rule as input to the variant calling rule.

---

Exercise 6

Context

BCF tools calls variants using a probabilistic model. An important parameter of this model is the expected substitution rate. We would like to add a parameter to the variant calling rule to easily specify this substitution rate.

Task

Add a non-file parameter to the variant calling rule to control prior substitution rate.

Details

• The prior substitution rate is given to bcftools call with the parameter -P.
• The default value is 0.001

Relevant Wiki section

Non-file rule parameters

---

Exercise 7

Context

We now have a VCF file containing all variants for our samples. We would now like to quantify each type of substitution (e.g. A->T, C->A, ...) to detect potential biases.

We provide the external script create_substitution_table.py at the root of the workflow folder.

If you are ahead, you can have fun writing the Python script yourself.

This script should parse the output of the variant calling rule and output a table of substitution counts for each combination of nucleotides. If you decide to implement the script yourself:

• Only consider SNPs, i.e. single nucleotide substitutions, ignore insertions and deletions
• A suggestion of final output would be a tabulated file (tsv) containing a matrix with the following structure:

	A	T	G	C
A	-	x	x	x
T	x	-	x	x
G	x	x	-	x
C	x	x	x	-

Task

Implement a rule to quantify SNP substitutions from the VCF file generated by the variant calling rule.

Details

• Call the external Python script create_substitution_table.py from the rule using the script directive.

Relevant Wiki section

Executing Python code

Hints

Parsing of VCF file

VCF files are complex and not entirely standardized. Fortunately, for this exercise, we only need very little information that is easily accessible: the reference nucleotide and the variant nucleotide. The first lines of the file are comments, starting with "##". You can skip these lines when parsing the file. Then, there is a header line starting with a single "#", that you can skip as well (you could skip all lines starting with "#").
Then, for each line in the file, we want to split the line with the separator TAB ("\t"); the reference allele will be in the 3rd field and the variant allele will be in the 4th field (0-based indexing).
To save the results while processing the file, you could store the counts for each substitution in a dictionary, or even a Defaultdict if you know how to use them.

---

Exercise 8

Context

We have a basic implementation of our entire workflow. We have seen how to run Snakemake to generate a given target, and in this case we would run snakemake <output_of_substitution_count_rule>. We would like to be able to automatically generate the final output when calling snakemake without arguments.

Task

Implement a special rule so that the final output is generated by default when running snakemake without specifying a target.

Hints

Special rules

Remember how the default target is selected if not specified in command line (see the Executing Workflows section). Also, keep in mind that a rule is not required to have an output!

Structure of the special rule

You can define a rule at the top of the Snakefile that will take as input the final output that you want to generate (from the last rule you wrote). This rule should not have any output or any shell command, only an input directive.

---

Exercise 9

Context

We have now implemented all the steps of our pipeline. We would like to visualize the process with Snakemake's built-in DAG visualization feature.

Task

Visualize the entire workflow’s Directed Acyclic Graph.

Details

• Save the results as a PNG or SVG picture.
• Specify the right targets to visualize the entire workflow!

Relevant Wiki section

Executing workflows

Hint

Command

Command: snakemake --dag | dot -Tpng > dag.png

Targets

To visualize the entire workflow, you need to specify the results of the substitution count rule.