HWperl

Materials: the lecture, Connecting from Windows, Editors

Files and setup

We recommend creating a directory (folder) for this set of tasks:

mkdir perl  # make directory
cd perl     # change to the new directory

The folder /tasks/perl/ contains input files, an example script working with FASTQ files and correct outputs for some of your scripts; see further details below. In this first homework we will provide almost all commands you need to run.

Write your protocol starting from an outline provided in /tasks/perl/protocol.txt. Make your own copy of the protocol and open it in an editor, e.g. kate:

cp -ip /tasks/perl/protocol.txt .  # copy protocol
kate protocol.txt &                # open editor, run in the background

Submitting

• Directory /submit/perl/your_username will be created for you
• Copy the required files to this directory, including the protocol named protocol.txt
• You can modify these files freely until deadline, but after the deadline of the homework, you will lose access rights to this directory

Task A (series)

Running the script from the lecture

• Consider the program for counting series in the lecture 1, save it to file series-stats.pl
• Open editor running in the background: kate series-stats.pl &
• Copy and paste text to the editor, save it
• Make the script executable: chmod a+x series-stats.pl
• Try running it on the small input: ./series-stats.pl < /tasks/perl/series-small.tsv
• You should get the following output:

Black Mirror 3
Game of Thrones 3

Extending the script

• Extend the script to compute the average rating of each series (averaging over all episodes in the series)
• Each row of the input table contains rating in column 5.
• Output a table with three columns: name of series, the number of episodes, the average rating.
• Use printf to print these three items right-justified in columns of a sufficient width; print the average rating to 1 decimal place.

Running the script

• If you run your script on the small file, the output should look something like this (the exact column widths may differ):

./series-stats.pl < /tasks/perl/series-small.tsv
        Black Mirror        3        8.2
     Game of Thrones        3        8.8

• Run your script also on the large file: ./series-stats.pl < /tasks/perl/series.tsv > series-stats.txt
• Check the output, e.g. by running cat series-stats.txt

Submitting

• Submit your script, series-stats.pl and the output series-stats.txt

Task B (FASTQ to FASTA)

Introduction

• In the rest of the assignment, we will write several scripts for working with FASTQ files introduced in the lecture. Similar scripts are often used by bioinformaticians working with DNA sequencing data.
• We will work with three input files:
  • /tasks/perl/reads-tiny.fastq a tiny version of the read file, including some corner cases
  • /tasks/perl/reads-small.fastq a small version of the read file
  • /tasks/perl/reads.fastq a bigger version of the read file

Goal

• Write a script which reformats FASTQ file to FASTA format, call it fastq2fasta.pl
  • FASTQ file should be on standard input, FASTA file written to standard output
• FASTA format is a typical format for storing DNA and protein sequences.
  • Each sequence consists of several lines of the file. The first line starts with ">" followed by identifier of the sequence and optionally some further description separated by whitespace
  • The sequence itself is on the second line, long sequences can be split into multiple lines
• In our case, the name of the sequence will be the ID of the read with the initial @ replaced by > and each / replaced by (_)
  • you can try to use tr or s operators (see also the lecture)
• For example, the first two reads of the file /tasks/perl/reads.fastq are as follows (only the first 50 columns shown)

@SRR022868.1845/1
AAATTTAGGAAAAGATGATTTAGCAACATTTAGCCTTAATGAAAGACCAG...
+
IICIIIIIIIIIID%IIII8>I8III1II,II)I+III*II<II,E;-HI...
@SRR022868.1846/1
TAGCGTTGTAAAATAAATTTCTAGAATGGAAGTGATGATATTGAAATACA...
+
4CIIIIIIII52I)IIIII0I16IIIII2IIII;IIAII&I6AI+*+&G5...

• These should be reformatted as follows (again only first 50 columns shown, but you include entire reads):

>SRR022868.1845_1
AAATTTAGGAAAAGATGATTTAGCAACATTTAGCCTTAATGAAAGACCAGA...
>SRR022868.1846_1
TAGCGTTGTAAAATAAATTTCTAGAATGGAAGTGATGATATTGAAATACAC...

Start programming

• You can start by modifying our script /tasks/perl/fastq-lengths.pl, which prints the length of each read in the FASTQ input file. You can use the following commands to start:

# copy our script to your folder under the new name
cp -i /tasks/perl/fastq-lengths.pl fastq2fasta.pl
# open in editor in background
kate fastq2fasta.pl & 

Running the script

• Run your script on the tiny file, compare with our precomputed correct answer:

./fastq2fasta.pl < /tasks/perl/reads-tiny.fastq > reads-tiny.fasta
diff reads-tiny.fasta /tasks/perl/reads-tiny.fasta

• Command diff prints differences between files. Here the output of diff should be empty. Otherwise try to look at the input and output files and fix your program to obtain the same output as we have.

• Also run your script on the small read file

./fastq2fasta.pl < /tasks/perl/reads-small.fastq > reads-small.fasta

Submitting

• Submit files fastq2fasta.pl, reads-small.fasta, reads-tiny.fasta

Task C (FASTQ quality)

Goal

• Write a script fastq-quality.pl which for each position in a read computes the average quality
• FASTQ file should be on standard input. It contains multiple reads, possibly of different lengths
• As quality we will use ASCII values of characters in the quality string with value 33 subtracted
  • For example character 'A' (ASCII 65) means quality 32.
  • ASCII value can be computed by function ord
• Positions in reads will be numbered from 0
• Since reads can differ in length, some positions are used in more reads, some in fewer
• For each position from 0 up to the highest position used in some read, print three numbers separated by tabs "\t": the position index, the number of times this position was used in reads, the average quality at that position with 1 decimal place (you can again use printf)

Example The first 4 lines when you run ./fastq-quality.pl < /tasks/perl/reads-tiny.fastq should be

0	5	28.2
1	5	33.0
2	5	33.0
3	5	33.2
4	4	33.5

• The last two lines when you run ./fastq-quality.pl < reads-small.fastq should be

99      86      5.5
100     86      8.6

Run the following command, which runs your script on the larger file and selects every 10th position.

./fastq-quality.pl < reads.fastq | perl -lane 'print if $F[0]%10==0'

• What trends (if any) do you see in quality values with increasing position?
• Submit only fastq-quality.pl
• In your protocol, include the output of the command and the answer to the question above.

Task D (FASTQ trim)

Write script fastq-trim.pl that trims low quality bases from the end of each read and filters out short reads

• This script should read a fastq file from standard input and write trimmed fastq file to standard output
• It should also accept two command-line arguments: character Q and integer L
  • We have not covered processing command line arguments, but you can use the code snippet below
• Q is the minimum acceptable quality (characters from quality string with ASCII value >= ASCII value of Q are ok)
• L is the minimum acceptable length of a read
• First find the last base in a read which has quality at least Q (if any). All bases after this base will be removed from both the sequence and quality string
• If the resulting read has fewer than L bases, it is omitted from the output

You can check your program by the following tests:

• If you run the following two commands, you should get file tmp identical with input and thus output of the diff command should be empty

./fastq-trim.pl '!' 101 < reads-small.fastq > tmp  # trim at quality ASCII >=33 and length >=101
diff reads-small.fastq tmp                         # output should be empty (no differences)

• If you run the following two commands, you should see differences in 4 reads, 2 bases trimmed from each

./fastq-trim.pl '"' 1 < reads-small.fastq > tmp   # trim at quality ASCII >=34 and length >=1
diff reads-small.fastq tmp                        # output should be differences in 4 reads

• If you run the following commands, you should get empty output (no reads meet the criteria):

./fastq-trim.pl d 1 < reads-small.fastq           # quality ASCII >=100, length >= 1
./fastq-trim.pl '!' 102 < reads-small.fastq       # quality ASCII >=33 and length >=102

Further runs and submitting

• ./fastq-trim.pl '(' 95 < reads-small.fastq > reads-small-filtered.fastq # quality ASCII >= 40
• Submit files fastq-trim.pl and reads-small-filtered.fastq
• If you have done task C, run quality statistics on the trimmed version of the bigger file using command below. Comment on the differences between statistics on the whole file in part C and D. Are they as you expected?

# "2" means quality ASCII >= 50
./fastq-trim.pl 2 50 < reads.fastq | ./fastq-quality.pl | perl -lane 'print if $F[0]%10==0'

• In your protocol, include the result of the command and your discussion of its results.

Note: in this task set, you have created tools which can be combined, e.g. you can first trim FASTQ and then convert it to FASTA (no need to submit these files).

Parsing command-line arguments in this task (they will be stored in variables $Q and $L):

#!/usr/bin/perl -w
use strict;

my $USAGE = "
Usage:
$0 Q L < input.fastq > output.fastq

Trim from the end of each read bases with ASCII quality value less
than the given threshold Q. If the length of the read after trimming
is less than L, the read will be omitted from output.

L is a non-negative integer, Q is a character
";

# check that we have exactly 2 command-line arguments
die $USAGE unless @ARGV==2;
# copy command-line arguments to variables Q and L
my ($Q, $L) = @ARGV;
# check that $Q is one character and $L looks like a non-negative integer
die $USAGE unless length($Q)==1 && $L=~/^[0-9]+$/;