HWperl

Materials: the lecture, Connecting from Windows, Editors

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 and last lines when you run ./fastq-quality.pl < /tasks/perl/reads-tiny.fastq should be

For example position 1 occurs in all reads and the qualities are B,A,B,A,D, which have ASCII values 66, 65, 66, 65, 68. After subtracting 33 from each and computing average we get 33. On the last position (21), we have only two reads with qualities A, which translates to value 32.

Running

Run the following commands, which runs your script on the large FASTQ file and selects every 10th position. The output is stored in reads-qualities.tsv and printed using cat

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

This input is a sample of real sequencing data from Illumina technology. What trends (if any) do you see in quality values with increasing position? Answer this question in your protocol.

Submitting

Submit fastq-quality.pl, reads-qualities.tsv

Task D (FASTQ trim)

Goal

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 the standard input and write trimmed FASTQ file to the 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

Testing 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

Use the following snippet, which stores command-line arguments 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]+$/;

HWperl

Contents

Files and setup

Submitting

Task A (series)

Task B (FASTQ to FASTA)

Task C (FASTQ quality)

Task D (FASTQ trim)

Navigation menu

Personal tools

Namespaces

Variants

Views

More

Search

Navigation

Tools