Wikiomics:WinterSchool day1

From OpenWetWare
Jump to navigationJump to search

The second day materials

Winterschool program day 1

Before we start

Who are the teachers

  • Darek Kedra (darek dot kedra at gmail dot com)
  • Emilio Palumbo

We both work at CRG (Centre for Genomic Regulation) in Barcelona, Spain. Roderic Guigo (DK & EP) and Cedric Notredame (DK) groups. Our work place looks like this: CRG image

Our main research topics start with RNASeq (DK, EP), scientific software development (EP), genome assembly and annotation (DK), SNP calling (DK). As for the farm animals we have been working with sequences mostly from chicken and cow, plus few species (duck, pig, sheep).

Organization etc.

  • lets turn out /silence the cell phones
  • our aim today is to teach you the bare necessities about command line and mapping Illumina reads to genomic DNA
  • we go through the page for about 50-60mins, then we switch to exercises for 20-30mins
  • no exposure to command line interface / Unix is assumed (we start from novice level)
  • this page will stay online after the course "forever", and we will improve it over time (old version will be accessible through History link on the top)
  • you are welcome to register as a wiki user and comment in the Talk page/fix any errors/ ask for enhancements
  • in the boxed parts of the page lines starting with # are comments, and what below are examples of commands. After typing equivalents of these in your terminal, press Enter

Introduction to Linux and the command line

Why Linux?

  • runs on everything from cell phones to supercomputers
  • long history of stable GNU/Unix tools
  • most of the bioinformatics software was written and intended to run on Linux

start vagrant or Linux from inside your virtual box machine, log in (Emilio)

vagrant instructions

go to relevant [age on the github:


#after downloading
tar xfvz Ex_1_2.tar.gz

Second set of files to download to vagrant:


tar xfvz Ex34.tar.gz

cd ex34/

logging in, connecting to other servers with ssh / sftp

As with other computers, one requires username password combination to connect to a specific computer. This combination can be specific to each of the computers or shared between i.e. all workstations at a given location.

SSH is a name for secure, encrypted connections between computers. It consist of two components, ssh server running on a remote machine and a ssh client on your laptop / workstation. The client is included in the default installations of recent Linux and OS X (Mac), but on Windows one has to install it ( ). It is not required you have it for the course, but it is a good idea to have it on your computer if you want to access remote servers on the command line.

With properly configured ssh connection (on Linux and Mac) you can run not just command line programs but also graphical user interface. This is also possible on Windows, but it is way more complicated to set up.

#to connect to a remote server, command line (we can not do it here, since we have closed ssh ports at this network

#to connect withe the same username as on the workstation you are already logged in:

#to connect with ability to run remotely programs requiting GUI
ssh -X

#to download single file from a host to a local directory
scp .

#to transfer a file from a local machine to a remote host:
scp local_filename

#to establish sftp connection 

#after establishing the connection we can list the content of remote directory:

#download  files
get remote_filename

file and directories naming

Linux is case sensitive, do MyFile.txt is different from MYFILE.TXT or myfile.txt. Try to use some consistent naming schemes for your project directories, input data or result files. You can use very long, descriptive names, like these result files from ENCODE project:


Things to keep in mind:

  • never ever use space/tabs in your file/directory names
  • never use Unix special characters in file/directory names (!?"'%&^~*$|/\{}[]()<>:)

Linux directory structure/navigating

You can see the Linux directory structure here

As with Windows and Mac systems one can imagine the systems of directories (synonym of folders) as a giant tree, where each of the branches has an upper directory and may contain lower level directories (sub-directories). At the top of this branched hierarchy it is just one place from where everything starts (no C:, D: X: drives), and in Unix speak it is called root directory. Some examples of directory naming:





Individual users have separate directories in which they store their data, results etc.

Some shortcuts to remember (what is after # sign is a description):

#root directory:

#your home directory:

#one directory above (two dots):

#current directory (one dot)

There are few commands to move between directories, create new ones, and list what is located there:

#print working directory = show where you are

#list what is located in the current directory

#list what is located in the current directory with details
ls -l 

#go to this directory
cd /home/linus/projects/chicken 

# special shortcut to go back to the directory you have been before
cd - 

#go to the root directory
cd / 

#go to your home directory 
cd ~ 

#go todirectory one level above the current one
cd .. 

#create new directory
mkdir mynew_directory 

#remove some directory (it must be empty!)
rmdir myold_directory 

absolute vs relative directory naming

We can do all the operations on directories or files using two conventions:

# absolute directory path starting with root "/"

ls /home/linus/projects/chicken/genome.fasta

#relative directory path (lets assume we are in  /home/linus/projects/banana/)

ls ../chicken/genome.fasta

This is useful shortcut for saving typing and avoiding errors when i.e. working on different systems.


In order to restrict actions associated with a given file or directory, each entity has 3 flags (r = read, w = write, e = execute) for 3 groups of users (file owner, group to which file owner belongs, i.e. students, and "all" for all the remaining users on that computer). So we have 9 fields describing what each of these groups can do with the file. On the top of it (or rather in the front of the string) we have another flag to tell us about what kind of thing it is ("-" = just a regular file, d = directory, l = link, etc.)

By default, owner has read+write permissions, his group members can read but not write (modify or delete) his files, and the rest of the word should have no right to see content of his files. These permissions are visible when listing content of a directory with "ls -l":

-rwxr-xr-x 1 vagrant users    49 Nov  7 14:43
-rw-r----- 1 vagrant users  1812 Nov  7 14:39 myfile01.txt
---------- 1 vagrant users  9045 Nov  7 14:41 myfile02.txt
-r--r--r-- 1 vagrant users  2016 Nov  7 14:39 myfile03.txt
-rw------- 1 vagrant users 67863 Nov  7 14:40 myfile04.txt
-rw-rw-rw- 1 vagrant users  8125 Nov  7 14:40 myfile05.txt
-rw-rw---- 1 vagrant users  7233 Nov  7 14:40 myfile06.txt
lrwxrwxrwx 1 vagrant users    12 Nov  7 14:50 myfile07.txt -> myfile06.txt
drwxr-xr-x 2 vagrant users     0 Nov  7 14:50 test_dir1

To change permissions, we have to specify first the group we want to modify, action (add permissions or remove them) then the permission themselves, finaly the name of the file(s):

#add execute permission to for all
chmod a+x  

# remove write permission for owner, group and the rest
chmod a-w 

This is often useful for sharing data / results with other users on the same machine or when writing scripts/executing some programs downloaded from the net.

copy, rename/move files, create symbolic links

Often we need to organize files in directories from the existing ones, which requires moving thing around, making copies, renaming directories and files. To save space or to avoid using long path names, we can create links to other existing files or folders.

CAVEAT: both cp and mv commands will perform the task without warning if the destination file exist. It is one of the ways of destroying your data. To avoid this error use "cp -i" or "mv -i" which will ask you for confirmation when you are overwriting existing file.

#copy move rename pattern: 
cp/mv source destination

#make a copy of  file1.txt and name it file2.txt
cp file1.txt file2.txt

#make a copy of  file1.txt with the same name but in a different directory
cp file1.txt /some/directory/of/your/choice/ 

#notice the space and the dot "." at the end of the command. The command is copy the file to the current directory
cp /home/linus/projecs/linus/file1.txt . 

Instead of doing copy (original file stays intact) one can do mv (move/rename) where only the new copy remains and the old one is destroyed.

#rename  file1.txt to file2.txt
mv  file1.txt file2.txt

#move  file1.txt to another directory
mv file1.txt /some/directory/of/your/choice/ 

#notice the space and the dot "." at the end of the command. The command is copy the file to the current directory
mv /home/linus/projecs/linus/file1.txt . 

Linking is a way of organizing your data. I.e. you have one big genome fasta file, and want to map your next generation sequencing reads using different mappers, which will require using this one file as an input to several programs. You can store your fasta file in one directory (i.e. /home/vagrant/projects/chicken/genome1.fa) and create links to this file in several subdirectories /home/vagrant/projects/chicken/bwa_mapping/, /home/vagrant/projects/chicken/bowtie2_mapping/, etc.). The command:

cd /home/vagrant/projects/chicken/bwa_mapping/
ln -s  /home/vagrant/projects/chicken/genome1.fa .

cd /home/vagrant/projects/chicken/bowtie2_mapping/
ln -s  /home/vagrant/projects/chicken/genome1.fa .

view files (more/less, head, tail), count (wc)

There are two basic types of files, text files and binary files. A lot of file formats in bioinformatics belong to simple text files, so the basic viewing/processing them is essential. Never ever do not try to open big text files (say few GB FASTQ file) in an editor, be it an editor on Linux or Microsoft Word...

Basic checking of the file content:

#displays content of the file screen by screen
more myfile04.txt 

#same as above but you can scroll back
less myfile04.txt 

#displays just the first 10 top lines of the file
head myfile04.txt 

#displays the first 100 lines 
head -100 myfile04.txt 

#displays just the last 10 lines of the file
tail myfile04.txt  

#displays  the last 100 lines of the file
tail -100 myfile04.txt  

wc myfile04.txt #newline, word, and byte counts for the file

# just count the lines of the file
wc -l myfile04.txt 

Using these commands is very useful for basic sanity check, so we can be sure that the FASTA file looks OK, that we have X lines in the file etc. First practical NGS application:

Q: how many reads do we have in our FASTQ file A: count the lines in FASTQ file, divide by 4 (each sequence occupies 4 lines)

compressing / uncompressing files (gzip, tar)

The sequence files are often really big, and we save space and time during downloads by compressing them. There are several compression programs, but the most frequently used is gzip:

#we will get some_file.gz, and the original uncompressed one will be removed
gzip some_file

# reversing the process, we get some_file in an uncompressed form
gunzip some_file.gz 

Program sources (almost always) or data are distributed not as single file but bundled in one file with "tar" program, then compressed for easier transfers/storage.

#we downloaded program_123.tar
#this will unpack the content of the tar archive
tar xfv program_123.tar 

#we downloaded program_123.tar.gz
#this will unpack the content of the compressed by gzip tar archive
tar xfvz program_123.tar.gz 

#we want to create tar archive containing all data_001.txt to data_100.txt (assuming there are only 100such files in this directory)
#tar cfv result.tar files_to_be_archived

tar cfv my_data.tar data_*.txt
#to create gzipped archive
tar cfvz my_data.tar.gz data_*.txt

checking the integrity of the files

While with fast internet access we get what we wanted most of the time (meaning: downloading data files), from time to time we can get a corrupt file for various reasons. Also when performing data analyses the common source of endless checks and rechecks is a data file with the same name but slightly different content. A great way to assure that the transferred data were not corrupted, or that we and our collaborators are talking about the same input data or results, is to perform file "fingerprinting", or calculate md5sum hashes.

md5sum file_2_fingerprint

md5 bioinfo.txt
MD5 (bioinfo.txt) = 06d5cd9ffb9154d7912350250e99951b

In short, this string is strongly influenced even by minuscule changes in the file (both text and binary, including software programs), and for practical reasons if 2 hashes are identical then the content of the files is also identical.

For more info about md5sum:

pattern matching, redirection and piping

Instead of using full name of a file each time, we can substitute part of the file name with special characters. Say we want to list the names all the files with .txt suffix ( myfile01.txt, myfile02.txt, etc. ) in a directory. And lets say we have about 500 of these:

ls *.txt

On Linux we can capture the output of one simple command (assuming it produces a text output) and do something useful with it, instead of just reading it from the screen. If we want to count how many such files we got, we can pipe "|" (vertical line, not the letter "I") the output of the ls command to wc command:

#display the number of .txt files in this directory
ls *.txt | wc -l 

We can also redirect the output to a file, creating i.e. list of txt files with names starting with some letter (redirection with single ">":

#we create a new file listing all the names of m*.txt files, one name per line
ls m*.txt > list_of_txt_files_starting_with_m.fof   

Or we can keep on adding the text to some file, which will be created if it does not exist (redirection with double ">>"):

#file is created and we get just a*.txt file names in it
ls a*.txt >> list_of_txt_files_starting_with_a_or_m.fof  

#file exists, so we are just appending   m*.txt file names to it
ls m*.txt >> list_of_txt_files_starting_with_a_or_m.fof  

#checking what is inside
more list_of_txt_files_starting_with_a_or_m.fof

search for strings / replace strings (grep & sed)

We can search the text files using command called "grep". Few examples:

#general pattern:
grep some_word text_file_name

#print all lines containing the word ITALY
grep ITALY my_file.txt 

#"reverse grep" prints all lines not containing ITALY
grep -v ITALY my_file.txt 

#just counts the number of lines containing ITALY
grep -c ITALY my_file.txt 

#prints all FASTA headers (lines starting with ">" sign)
grep "^>" multiple_sequences.fasta

# prints number of sequences in the multiple fasta file
grep -c "^>" multiple_sequences.fasta 

Instead of opening an editor and doing search and replace, we can do it on a command line, using utility called sed, and operate both on files or use pipes.

# replaces each occurrence of ITALY with SPAIN
sed 's/ITALY/SPAIN/g' my_file.txt 

# we get all FASTA headers and then remove the ">" sign, getting just the proper FASTA IDs
grep "^>" multiple_sequences.fasta  | sed 's/>//g' 

awk in few minutes

Awk is a simple programing language used for basic text processing. It is still being used because often it is faster to write and execute a command in awk than write small small script in more advanced languages, such as perl, or python. The basic concept is that awk splits the lines of the text into individual words/numbers (== great for text files in a column form), which we can then access using $column_number notation

# prints just first column of the file
awk '{print $1}' my_data.txt 

# prints first column but also "_my_new_label" string just after it
awk '{print $1"_my_new_label"}' 

We can perform simple calculations, (assuming we have columns of numbers).

#this will simply sum all the numbers from the first column then print that number
awk '{ sum+=$1} END {print sum}' my_data.txt 

working with Linux processes

Every program running on Linux has a special number given to it by the operating system. We can view what is running on a given machine using following commands:

#list what the user executing the command is running

#list all processes running on the system and display them in a full form:
ps -aF
#ps -aF output
frodo   798     1  0  6069   696   2 Oct29 pts/54   00:00:00 dbus-launch --autolaunch 2b4168a54c58a6631d2304110000001c --binary-syntax --close-stderr
133003    3584  1532  1 161322 412996 0 08:50 pts/65   00:03:09 /users/GR/mb/dsoronellas/bin/lib64/R/bin/exec/R
guigo   5040     1  0  6069   812   9 Nov13 pts/56   00:00:00 dbus-launch --autolaunch 2b4168a54c58a6631d2304110000001c --binary-syntax --close-stderr
guigo  5206  4294  0 103493 22668  0 09:18 pts/64   00:00:06 gedit
guigo  5210     1  0  6069   812   1 09:18 pts/64   00:00:00 dbus-launch --autolaunch 2b4168a54c58a6631d2304110000001c --binary-syntax --close-stderr
henry  6341  6282  0 26524  1300  17 09:33 pts/67   00:00:00 /bin/bash scripts/ scripts/ result.txt results/Drosophylas
henry  6345  6341 99 32127 17092  17 09:33 pts/67   03:54:44 python scripts/ results/Drosophylas_Adu-RNAseq_Analysis/Ann_FlyBase/res_Drosophi

#show the machine load and the jobs taking most CPU 
#top output

top - 13:34:13 up 108 days,  2:51, 53 users,  load average: 2.06, 2.10, 2.13
Tasks: 1063 total,   3 running, 1043 sleeping,  16 stopped,   1 zombie
Cpu(s):  6.1%us,  0.1%sy,  0.0%ni, 93.8%id,  0.0%wa,  0.0%hi,  0.0%si,  0.0%st
Mem:  132119396k total, 130992952k used,  1126444k free,   380468k buffers
Swap:  4194296k total,   100488k used,  4093808k free, 126879876k cached

  PID USER      PR  NI  VIRT  RES  SHR S %CPU %MEM    TIME+  COMMAND                                                                                                              
20737 frodo     20   0 54764  35m  224 R 100.0  0.0  24137:30 ld-linux-x86-64                                                                                                      
 6345 toad  20   0  125m  16m 1604 R 93.0  0.0 238:58.55 python                                                                                                                
45922 linus  20   0 1131m 149m  58m S  1.0  0.1  41:01.38 soffice.bin                                                                                                           
53551 linus    20   0 20032 2144 1056 R  0.7  0.0   0:00.06 top                                                                                                                   
53586 sambu     20   0     0    0    0 Z  0.7  0.0   0:00.02 >                                                                                                      
  131 root      20   0     0    0    0 S  0.3  0.0  14:47.98 events/0  

Sometimes we can have a process running and taking all the resources, writing huge output files or just producing the wrong results and one which we can not just quit. We can kill these (terminate) using following:

#kill process based on its PID (process ID). You have to own the process or be a root
kill -9 20737

#kill all programs with given name
killall -9 firefox

Other very useful command is nohup. It prevents the program from being killed when the terminal from which we started is closing for some reason, like a falling remote ssh connection.

#the output goes to nohup.out
nohup ./ 

Where to go from there

more information

All the commands used above have multiple options impossible to cover during the short introduction. Once you know the name of the program, you can learn about it functions and usage reading installed manual:

man ls
man grep

The recommended reference book is:

Unix Power Tools, 3rd Edition

By Jerry Peek, Shelley Powers, Tim O'Reilly, Mike Loukides

Publisher: O'Reilly Media

Released: October 2002

Pages: 1156

shell scripts

Sometimes we need to do the same thing on a set of files, like submitting a hundred jobs to a computer cluster using the same program (lets say for mapping FASTQ files) but with different input data (FASTQ files) and differently named outputs (BAM files). For doing it we need to first create 100 individual "program files" (== shell scripts), here called "map_script" i.e. using python (next section). But instead of submitting the jobs one by one we can do instead on a command line:

for file in map_script*. #press enter
do #press enter
qsub $file #press enter
done #press enter


Despite that there is a lot of bioinformatics software written in Perl (ENSEMBL, GBrowse, BioPerl), Python is a good choice as a first (and maybe even only) programming language you need for day to day bioinformatics tasks. It is easier to learn and read than Perl, but still fast enough to deal with i.e. FASTQ files.

There are two versions of Python, 2.7 and new Python 3 series which differ in subtle ways (Python 3 is an improvement). Because of the legacy software, it is still safer for time being to use the older version. One thing which distinguishes Python from other languages is that spaces/tabs/line ends are syntactic elements crucial for the language. In short, instead if using () to mark the beginning and the end of some function, Python uses indentations (tab or spaces, by default 4) to mark code blocks. With almost any programming editor, keeping track of these is painless.

  • multiple file mover/renamer
#!/usr/bin/env python


This is just and example how to rename a bunch of files using external mv
Renames all mapping_result_001.bam to mapping_result_001.bwa_cow_ENS73.bam
import glob, os

#in the pattern below the ??? stand for any 3 characters which must occur at that position in file name  
input_file_pattern = "mapping_result_???.bam"

#glob.glob below is a way to get back a list of files matching the pattern. 

for input_file_name in glob.glob(input_file_pattern):
    #here we remove the last 3 characters from the input_file_name, then we add the "bwa_cow_ENS73.bam" to it
    output_file_name = input_file_name[:-3] + "bwa_cow_ENS73.bam"
    #here we construct the whole command to be executed. Observe extra space needed between names if the input and output files
    command = "mv " + input_file_name + " " + output_file_name
    #this is an idiom using os library to execute command

  • Fasta header fixer
#!/usr/bin/env python
Comment: this script simplifies names of sequences in FASTA file limiting them to just first name
Usage: some_fasta_file_I want_to_transform > some_fasta_file_I want_to_transform.fixed_file


import sys
input_fasta_fn = sys.argv[1]

#input_fasta_fn ="some_fasta_file_I want_to_transform"

for line in open(input_fasta_fn).readlines():
   if line[0] == ">":
       new_fasta_header = line.split()[0]
       print new_fasta_header
       print line[:-1]   

computer clusters

Individual workstations or even multi CPU servers are not good enough for running hundred of CPU-demanding jobs. For that we have computer clusters, where we use tens or hundreds of computer nodes. To distribute loads on such clusters we use special software maintained by system admins, i.e. Sun Grid Engine (SGE). It takes care that not all jobs are run at once the moment they are submitted by the user, and it distributes the jobs among the nodes balancing the loads. Because individual cluster installations at your institutions are likely different from each other (they may use different software than SGE), this section is mostly for convincing you that whenever you have some cluster you can use, you should investigate how you can use it and run your next generation sequencing jobs on them. There are few concepts you need to know about typical cluster configuration and how it needs to be taken into account:

  • clusters are often partitioned into smaller entities, so for example one person from one group can not block the use of the cluster for the whole institution. These partitions are being called queues.
  • cluster jobs may have differer in their requirements. I.e genome assembly jobs often require large memory servers, and may take a week or more to complete. On the other end we can have tasks which are trivial to split into say 50 separate jobs by the user, and which will finish in minutes or few hours. So cluster admins often create separate queues for such different tasks, and may reserve small group on nodes for each research group, so the most urgent/or the longest running tasks can be run that way.
  • It is a good policy to require that all individual jobs executing threads/processes on several CPUs in parallel should be labeled as such and wait their turn until node with required number of CPUs becomes free.
  • same thing is valid for memory usage
  • often jobs without information about number of CPUs/RAM requirements will have set some default values (i.e. 1 CPU) by the admins

Each job run on the cluster creates two files, capturing the output of the job which normally is printed on the screen. One is for standard messages ("mapping reads from file XYZ") other for error messages ("file XYZ not found, terminating"). These are crucial for debugging. The naming convention is as follows:

#two jobs we submitted to run

#These will receive job numbers once they start running, ie:

#the standard output files will be names:

#the standard error files will be names:

#if we re-run these scripts, the job numbers will be different, so we can get:

Caveat: these files may be created in your home directory, not where the scripts and data are, when you forgot to specify -cwd during submitting your jobs

#submit a single job to a cluster to the mygroup_cluster_queue queue (with default values). 
#Change to the current working directory (one from which the script was submitted
qsub -q mygroup_cluster_queue -cwd 

#monitor the status of your jobs on the cluster

#the typical output will look like this:
job-ID  prior   name       user         state submit/start at     queue                          slots ja-task-ID 
2178314 0.05226 LiverE14Al epalumbo     r     10/31/2013 20:00:31      8        
2178315 0.05226 LiverE14Al epalumbo     r     10/31/2013 20:00:31      8        
2178316 0.05226 LiverE14ha epalumbo     r     10/31/2013 20:00:31      8        
2178317 0.05226 LiverE14ha epalumbo     r     10/31/2013 20:00:31 rg-el6@node-hp0104rg.linux.crg     8   

#delete job with a number 23456 submitted to a cluster 
qdel 23456

Tip: it makes sense that when submitting different groups of jobs to a cluster you use different script names. If something goes wrong it will be easier to find jobs without a chance of running to completion, nd remove them.

Lets say we have 50 jobs bwa_map_001 to bwa_map_050 and another 30 jobs bowtie_map_001 to bowtie_map_030. We want to remove all our still running or submitted and waiting bowtie_map jobs, but do not touch bwa_map jobs. Here is the command:

qstat | grep myuser_name | grep bowtie_map | awk '{print $1}' | xargs qdel

Explanation of the steps:

  • qstat: we list our jobs running or waiting
  • grep myuser_name: just keep the lines listing the jobs, not qmon extra lines
  • grep bowtie_map : we keep lines containing that string (our jobs to delete). Keep in mind that qmon does not have space to display print the whole name, so we have to check first that such string is present in qmon listing
  • awk '{print $1}': we get just the first qmon column (job number)
  • xargs : special command telling the next command that the parameter will be passed from the pipeline
  • qdel : delete cluster job with a given job number


Format and quality checks

Already in the 90ties when all sequencing was being done using Sanger method, the big breakthrough in genome assembly was when individual bases in the reads (ACTG) were assigned some quality values. In short, some parts of sequences had multiple bases with a lower probability of being called right. So it makes sense that matches between high quality bases are given a higher score, be it during assembly or mapping that i.e. end of the reads with multiple doubtful / unreliable calls. This concept was borrowed by Next Generation Sequencing. While we can hardly read by eye the individual bases in some flowgrams, it is still possible for the Illumina/454/etc. software to calculate base qualities. The FASTQ format, (usually files have suffixes .fq or .fastq) contains 4 lines per sequence:

  1. sequence name (should be unique in the file)
  2. sequence string itself with ACTG and N
  3. extra line starting with "+" sign, which contained repeated sequence name in the past
  4. string of quality values (one letter/character per base) where each letter is translated in a number by the downstream programs

Here it is how it looks:

@SRR867768.249999 HWUSI-EAS1696_0025_FC:3:1:2892:17869/1

Unfortunately Solexa/Illumina did not follow the same quality encoding as people doing Sanger sequencing, so there are few iterations of the standard, with quality encodings containing different characters. For the inquisitive:

What we need to remember from it, that we must know which quality encoding we have in our data, because this is an information required by mappers, and getting it wrong will make our mappings either impossible (some mappers may quit when encountering wrong quality value) or at best unreliable.

There are two main quality encodings: Sanger and Two other terms, offset 33 and offset 64 are also being used for describing quality encodings:

  • offset 33 == Sanger / Illumina 1.9
  • offset 64 == Illumina 1.3+ to Illumina 1.7

For that, if we do not have direct information from the sequencing facility which version of the Illumina software was used, we can still find it out if we investigate the FASTQ files themselves. Instead of going by eye, we use a program FastQC. For the best results/full report we need to use the whole FASTQ file as an input, but for quick and dirty quality encoding recognition using 100K of reads is enough:

head -400000 my_reads.fastq > 100K_head_my_reads.fastq 
fastqc 100K_head_my_reads.fastq
#we got here 100K_head_my_reads.fastq_fastqc/ directory

grep Encoding 100K_head_my_reads.fastq_fastqc/fastqc_data.txt

Encoding	Sanger / Illumina 1.9

CAVEAT: all this works only on unfiltered FASTQ files.. Once you remove the lower quality bases/reads containing them, guessing which encoding format is present in your files is problematic.

Here is a bash script containing awk oneliner to detect quality encoding in both gzip-ed and not-compressed FASTQ files.


if [[ $file ]]; then
    if [[ $file =~ .*\.gz ]];then
    command="$command $file | "

command="${command}awk 'BEGIN{for(i=1;i<=256;i++){ord[sprintf(\"%c\",i)]=i;}}NR%4==0{split(\$0,a,\"\");for(i=1;i<=length(a);i++){if(ord[a[i]]<59){print \"Offset 33\";exit 0}if(ord[a[i]]>74){print \"Offs
et 64\";exit 0}}}'"

eval $command

Types of data

  • read length

from 35bp in some old Illumina reads to 250+ in MiSeq. The current sweet spot is between 70-120bp.

  • single vs paired

Just one side of the insert sequenced or sequencing is done from both ends. Single ones are cheaper and faster to produce, but paired reads allow for more accurate mapping, detection of large insertions/deletions in the genome.

Most of the time forward and reverse reads facing each other end-to-end are

  • insert length

With the standard protocol, the inserts are anywhere between 200-500bp. Sometimes especially for de novo sequencing, insert sizes can be smaller (160-180bp) with 100bp long reads allowing for overlap between ends of the reads. This can improve the genome assembly (i.e. when using Allpaths-LG assembler requiring such reads). Also with some mappers (LAST) using longer reads used to give better mappings (covering regions not unique enough for shorter reads) than 2x single end mapping. With paired end mappings the effects are modest.

Program for combining overlapping reads: FLASH:

For improving the assembly or improving the detection of larger genome rearrangements there are other libraries with various insert sizes, such as 2.5-3kb or 5kb and more. Often sequencing yields from such libs are lower than from the conventional ones.

  • stranded vs unstranded (RNASeq only)

We can obtain reads just from a given strand using special Illumina wet lab kits. This is of a great value for subsequent gene calling, since we can distinguish between overlapping genes on opposite strands.

quality checking (FastQC)

It is always a good idea to check the quality of the sequencing data prior to mapping. We can analyze average quality, over-represented sequences, number of Ns along the read and many other parameters. The program to use is FastQC, and it can be run in command line or GUI mode.

  • good quality report:

  • bad quality FastQC report

trimming & filtering

Depending on the application, we can try to improve the quality of our data set by removing bad quality reads, clipping the last few problematic bases, or search for sequencing artifacts, as Illumina adapters. All this makes much sense for de novo sequencing, were genome assemblies can be improved by data clean up. It has a low priority for mapping, especially when we have high coverage. Bad quality reads etc. will simply be discarded by the mapper.

You can read more about quality trimming for genome assembly in the two blog posts by Nick Loman:

Trimmomatic From the manual:

Paired End:

java -jar trimmomatic-0.30.jar PE --phred33 input_forward.fq.gz input_reverse.fq.gz output_forward_paired.fq.gz output_forward_unpaired.fq.gz output_reverse_paired.fq.gz output_reverse_unpaired.fq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

This will perform the following:

Remove adapters Remove leading low quality or N bases (below quality 3) Remove trailing low quality or N bases (below quality 3) Scan the read with a 4-base wide sliding window, cutting when the average quality per base drops below 15 Drop reads below the 36 bases long Single End:

java -jar trimmomatic-0.30.jar SE --phred33 input.fq.gz output.fq.gz ILLUMINACLIP:TruSeq3-SE:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

This will perform the same steps, using the single-ended adapter file

Tagdust (for simple unpaired reads)

Tagdust is a program for removing Illumina adapter sequences from the reads containing them. Such reads containing 6-8 bases not from genome will be impossible to map using typical mappers having often just 2 mismatch base limit. Tagdust works in an unpaired mode, so when using paired reads we have to "mix and match" two outputs to allow for paired mappings.

tagdust -o my_reads.clean.out.fq -a my_reads.artifact.out.fq adapters.fasta my_reads.input.fq

Error correction

For some applications, like de novo genome assembly, one can correct the sequencing errors in the reads by comparing them with other reads with almost identical sequence. One of the programs which do perform this and are relatively easy to install and make it running is Coral.


web site:

version: 1.4

It requires large RAM machine for correcting individual Illumina files (run it on 96GB RAM)

#Illumina reads
./coral -fq input.fq -o output.fq  -illumina  

#454 reads
./coral -fq input.454.fq -o output.454.fq  -454  

source of published FASTQ data: Short Read Archive vs ENA

While we will often have our data sequenced in house/provided by collaborators, we can also reuse sequences made public by others. Nobody does everything imaginable with their data, so it is quite likely we can do something new and useful with already published data, even if treating it as a control to our pipeline. Also doing exactly the same thing, say assembling genes from RNASeq data but with a newer versions of the software and or more data will likely improve on the results of previous studies. There are two main places to get such data sets:

go there 
put mouse RNASeq
417 public access sets
Click on it,
it looks like we got just: RNA (348)
Go there
put Bos taurus
see on the right table SRA Experiments 636

on the left:
DNA (171)
RNA (454)
metagenomic (13)

  • European Nucleotide Archive

Go there:
put RNAseq
see Experiment (5)
put RNA-seq
see Experiment (109)

Which one to use? ENA is easier as you get gzipped fastq files not SRA archives requiring extra processing, sometimes painful (at one stage the funding for SRA programs was cut). But NCBI tools may have better interface at times, so you can search for interesting data set at NCBI, then store the names of experiments and download fastq.gz from ENA.

Genomic fasta and gtf/gff gene annotation

ENSEMBL together with UCSC are two "clearing houses" when it comes to genome assembly and annotation. For mot of the genomes at least the genome versions (== sequences) are identical, but from time to time we can encounter a switch, like with bovine genome (UMD 3.1 @ENSEMBL and Baylor Btau 4.6.1 @UCSC). While it can be that you will have several releases of ENSEMBL with unchanged genome and annotation, it is still good idea to always include in the protocol / directories & file names the version used for the mappings/analyzes.

Side note

Since we are EU project, we tend to use ENSEMBL for the vast majority of the analyzes. There are some tools / data which are easier to use at UCSC, i.e. liftOver utility, for "mapping" genomic intervals from one mammalian genome (usually human or mouse) to another one (like cow). This is useful when we got large set of human genomic features (say ChIPSeq peaks) and want to find corresponding regions in cow. Check:


For getting the genome and gene annotations:

1. go to
2. put in white on blue bar search window (default text "Filter") i.e. "chicken". 
# you will get just one row with multiple entries, two of our interest at the moment:

* Gene sets (GTF)

3. Click on DNA (FASTA)
# in the directory:

we are interested in: 

Do not click on it as doing it 50x on slow WiFi will kill the connection.

4. go back in the browser to (with just chicken row)
5. click on GTF
# you get to
file to download (again NOT NOW!) is: Gallus_gallus.Galgal4.73.gtf.gz

Keep in mind that there are several species which were not included for various reasons in the main ENSEMBL (currently 73) release, but are accessible in preENSEMBL:

This time not every species has downloadable genome, there are just few available here:

As for GTFs it is even more unordered set, with all annotations laying around:

Still, in times of need this is probably better than no genome & no annotation.

grepping fasta headers

This was covered in the grep session, so just a quick reminder:

#count number of sequences in multiple fasta
grep -c "^>" multi_fasta.fa

#list sequence names in fasta, display screen by screen with "q" to escape:
grep "^>" | less 

fasta reformat from exonerate

Most likely well written program for short read mapping should cope with long fasta headers:

>GI12345|QW76543 bovine contig from Angus containing genes X, Y and Z

or sequences with one chromosome per line. Still, it is in general a good idea to simplify and reformat multiple fasta genomic draft sequences to avoid potential problems problems down the road. For fixing FASTA headers (assuming that the first part of the header is unique), use Python script from the above Python section.

#for this you have to have both python script and the raw fasta file in one directory
./ ref_gen.raw.fa > ref_gen.fixed.fa

For reformatting sequences so that all sequences will have 70 columns:

fastareformat input.fa > output.fa

GTF / GFF annotation formats

For storing information about multiple, often genomic features the most comonly used format are few versions of GFF, in particular GTF. The GTF (General Transfer Format) is identical to GFF version 2.

Section adopted from:

GTF Fields:

seqname - name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix.
source - name of the program that generated this feature, or the data source (database or project name)
feature - feature type name, e.g. Gene, Variation, Similarity
start - Start position of the feature, with sequence numbering starting at 1.
end - End position of the feature, with sequence numbering starting at 1.
score - A floating point value.
strand - defined as + (forward) or - (reverse).
frame - One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..
attribute - A semicolon-separated list of tag-value pairs, providing additional information about each feature.
Sample GFF output from Ensembl export:

===Example ===

X	Ensembl	Repeat	2419108	2419128	42	.	.	hid=trf; hstart=1; hend=21
X	Ensembl	Repeat	2419108	2419410	2502	-	.	hid=AluSx; hstart=1; hend=303
X	Ensembl	Repeat	2419108	2419128	0	.	.	hid=dust; hstart=2419108; hend=2419128
X	Ensembl	Pred.trans.	2416676	2418760	450.19	-	2	genscan=GENSCAN00000019335
X	Ensembl	Variation	2413425	2413425	.	+	.	
X	Ensembl	Variation	2413805	2413805	.	+	.

End of day 1