Preprocessing of NGS reads: Trimming and filtering

The last few MGRC seminars I attended has been emphasizing the importance of pre-processing of NGS reads. SNPs detection relies heavily on read quality. Softwares like Maq and Samtools use quality score. But it's a different story when it comes to assembly because assembly doesn't use quality score. Some think that trimming is only necessary in very high coverage data.

Only the past few months, more tools on read quality assessment are available for the wider community. A recent IlluminaGAII 1.3+ Pipeline documentation reported that a run of bases with a quality score 2 or symbol 'B' is unreliable and should be removed (View discussion here). With increasing understanding and the tools available, I think the community can make a better decision to trim or not to trim.

Here is a list of trimming tools/scripts I found online:
1. HT Sequence Analysis with R and bioconductor
2. HTSeq by Simon Anders.
3. Softtrim.R by Jeremy Leipzig
4. TrimBWAstyle.pl by Joe Fass
5. Solexa_Sig2
6. Biopython

Bioconductor ShortRead package has been available since 2008. It requires a bit of R knowledge. It's only recently I realized that it has a function to trim reads to desired length. It only can perform right trimming.

HTSeq is able to trim adaptors but not removing low quality reads. The ht-seq-qa script is particularly useful. It generate a nice plot to show distribution of quality score over read position.

Here's how I perform trimming:
First, I installed Bioconductor ShortRead package. Make sure you use the latest version of R. Please note that trimming using R is RAM intensive. In my case, it took an hour to process a 1Gb file using 12Gb RAM computer. The program sometimes get killed or return an error message saying "Error: cannot allocate vector of size 1.5 Gb".

Then, convert the quality score of your fastq file to Sanger Phred score using IllQ2SanQ.pl script from UC Davis. Next, I used Leipzig's Softtrim R script. I can choose the min quality score, min read length of trimmed reads and enable left trimming. This script works better for me. Click here to see how to use his script. Lastly, I separated the trimmed reads into paired and single reads using a Python script.

Here's how my trimmed reads look like:

Before trimming

After trimming

*Plots are generated using htseq-qa script by Simon Anders

Read more...

May 2010: NGS Seminar & lecture

There are two important Next Generation Sequencing events this month:

MGRC Lecture
Title: Application of Next-Generation Sequencing Technologies as Shared Research Resources
Date : 25 May 2010
Time : 2 - 5 pm
Venue : Plenary Theatre, Level 3, KL Convention Centre

Illumina South Asia Pacific Seminar
Date: 31 May 2010
Time: 8.45am
Venue: Cempaka Room, Hotel Equatorial, Bangi
Organizer: Sciencevision

Read more...

Bioinformatics for Biologist: Installing NGS software

I remember I have difficulty install Velvet, the first Next Generation Sequencing software I use. As a biologist, I just don't understand how to install despite reading the manual many times. Once I did, NGS software installation becomes simple. Occasionally I run into problems. Maq and MapView are among the hardest softwares to install. Here, I'm gonna share how to install a NGS software.

Installation
The NGS softwares are compressed in a tar.gz file. Once downloaded, right click and select "Extract here". A folder with the same name will appear. It's your choice to install in which directory. Preferably in a hard disk (not Desktop or Documents) as the file generated will take a lot of space.

To install, ALWAYS read the manual or README.txt for instructions. Instructions differ for each software. Some softwares like Samtools are ready to run once unzipped. The two basic commands to install are:
> make
> make install

Here's an example, I just downloaded BWA 0.5.5 from http://sourceforge.net/projects/bio-bwa/files/. After unzipping it, I found no instruction to install. I thought it's ready to use. But command "./bwa" doesn't work. So I tried "make install" and "make". Only the latter works. During installation, there will be a long list in the Terminal to check if the required files/libraries/systems are present.

Once the installation completed, I type "./bwa". A help file with a list of options appeared. In bwa-0.5.5 folder, you will notice many files are added during installation. Mission accomplished.


Troubleshooting
When you fail to install it, the Terminal will show "ERROR" or "Stop". If you wanna know what went wrong, you need to check the lines that Terminal produces during Installation. Maybe you didn't install a library or some GNU environment required. Once missing components are identified, go to System > Administration > Synaptic package manager to install them. Try installation again.

Tutorial
Let's try installing maq-0.7.1. I have problem installing this version eventhough I was certain I have all the components installed. I even tried to install an older version 0.6.8 but it didn't work. Took me awhile to figure out what's missing.

For installation instruction, read http://maq.sourceforge.net/maq-man.shtml#install. To download, http://sourceforge.net/projects/maq/files/maq/

Why don't you try it on your own first. Just download the zip file, unpack it and simply type on Terminal:

./configure; make; make install

Alternatively, you can also type 3 separate commands:
>./configure
>make
>make install

If you manage to install it, congratulations! If you don't, welcome to the club!

OK. Let's cut the chase and tell you what I did to fix the problem. I realized a library known as "zlib1g-dev" is missing. After installing that library, I type:
>./configure
>make
>sudo make install
Then, type in administrator password and Wala!

Since I already install the library needed, I can't show you an example of failed installation. Oh this is the worst tutorial ever!

Read more...

Bioinformatics for biologist: Using terminal

One of the first things in Bioinformatics is learning how to use Operating Systems (OS) other than Windows such as Linux and Mac. I'm currently dual-boot 64-bit Ubuntu with Windows Vista. What's the different between 64-bit and normal 32-bit? 64-bit OS enable more than 4 GB RAM usage. But don't bother buying 64-bit Windows because most NGS tools are built and tested on 64 bit Linux OS. My advice is to get a free 64-bit Linux distro which is suitable for you.

And, most tools run on command line instead of Graphic User Interface(GUI). In fact, many bioinformatians think it's a stupid idea to create a program in GUI. Who needs GUI? Biologist does.

The command line is run using a shell. The shell in Ubuntu is called Terminal. To open Terminal, click on Applications menu -> Accessories -> Terminal.

Let's look at some of the basic operations I use daily:

cd
To change directory. Directory can be refered to as Folder in Windows term.
cd ..
Back to previous directory
cd /
Back to root directory

ls
list the file/directory in the current directory.
ls *.pl
List all the file/directory ending with ".pl". For example, fastq2fasta.pl will be listed.
ls fastq*
List all the file/directory starting with "fastq".
It's easy to find a file if you cannot remember part of the file name.


rm
Remove a file or directory

wc -l filename
Count the number of lines in a file. Useful to find out the number of sequences in a FASTA file. # sequences = # lines divided by 2

split -l 20000 filename
Create smaller files consisting of 20000 lines each. Type "split --help" for more split options.

more filename
Read the content of the file from beginning. This is very useful when the file size is very big (>100MB). Opening it using a text editor will crash the system. Press ENTER to read more lines.

top
Display information about system, running processes and RAM usage. To exit, press 'q'.

sudo
Root user or administrator. Any command starting with "sudo" requires administrator password.
sudo rm -rf ~/.local/share/Trash/files/*
Administrative privilege used to empty trash.

When typing a long command on terminal, you would prefer to paste that command there. Note that Terminal is different. To copy something from Terminal, hightlight and press shift+ctrl+c. To paste onto Terminal, press shift+ctrl+v. Alternatively, you can right-click your mouse and choose the copy/paste option.

When a command or application finish running, you will see the directory address with "~$" sign again. To stop any running application, press shift+ctrl+z.

Find out more on https://help.ubuntu.com/community/UsingTheTerminal

Read more...

ABI SOLiD softwares: SOLiD Bioinformatics Workshop Part III

During the workshop, Dr. Sun presented some unpublished data on de novo genome assembly of rice blast fungus with a size ~40Mb. Partek Genomic Suite which is a third party software that can handle SOLiD data has been demonstrated. Generally, it performs gene expression analysis. It can run on all OS and requires minimum RAM.

We were then introduced to SOLiD Software Community. SOLiD suers can download a variety of data analysis tools from http://solidsoftwaretools.com for free. Tools available includes secondary and tertiary analysis softwares. Datasets such as human small RNA, whole transcriptome, E.coli single reads and PE reads can be found in that website. I got a chance to play around with SOLiD de novo assembly tool using E.coli dataset. The software is user friendly. I imagined working with color space is a very complicated task but it's easier than I expected. This software automatically display base read after assembly. There is a Find option to locate the reads which I find pretty useful. Another advantage is that the softwares can be installed on all types of OS (which is great for Windows user like me). Click here to read application note on SOLiD software suite.

Here's a list of other SOLiD softwares
• Academic and open source - Velvet, MAQ, SOCS, Bowtie, BFAST, SHRiMP, PASS
• commercial softwares - NextGENe, CLC Bio, ZOOM, Partek
(Most resources come from seqanswers.com and lh3)


Related posts:
SOLiD Bioinformatics Part I : Introduction
SOLiD Bioinformatics Part II: de novo assembly

Read more...

De novo NGS Assembly: SOLiD Bioinformatics workshop part II

Sorry it took so long to post Part II of SOLiD bioinformatics workshop. Things that were puzzling to me two months ago all make sense now. I'm in debt with Dr. Sun JianDong for his kind explanation.

Dr. Sun spoke about de novo assembly. As we know, de novo assembly is fragment assembly without reference. Genome sequences and existing ESTs databases can be used as reference although the former has more advantages. Most assemblers used for Sanger sequencing cannot handle short reads generated by Next Generation Sequencing. Because of the length, short reads must be produced in large quantities and at greater coverage depths (Velvet 2008). Since the introduction of NGS, the science community has been very concerned with the development of algorithms that are suitable for these reads. Generally, there are TWO types of de novo NGS assemblers.

Hamilton Path

• also known as overlap-layout-consensus approach (Batzoglou 2005).
  
• each node is represented by one read and each detected read overlap as an arc between the appropriate nodes (Zerbino & Birney 2008)
• not suitable for short reads
• Examples: SSAKE, SHARCGS, SHORTY, Edena
  

Eulerian approach

• De Bruijn graph
• less complex and more accurate
• very sensitive to errors and low quality reads
• reads are mapped to the path based on k-mers.
• Examples: Velvet, Euler-SR
  

Programs that use De Bruijn graph are more accurate and faster compared to Hamilton path. However, these programs cannot tolerate the slighest error rates like 0.3%. (Chaisson et al 2009) Therefore, erroraneous ends must be trimmed and errors correction must be done. Latest version of Euler SR: Euler USR can assemble error-prone reads. Among all NGS assemblers, Velvet is still the widely used programs due to its fast and efficient assembly.

Next, he explained several metrics used in assembly evaluation:

• N50 contig length - Longer size is wanted
  
• Number of contigs - less number of contigs is desirable
  
• Length of contigs - longer contigs is better
  
• Coverage - higher coverage, better assembly.
  

(Note: N50 contig length = The size of contig such that 50% of the assembly is contained in contigs size of N50 or greater)

Here's another definition of N50. If we sort contigs from the largest to the smallest, and start covering the genome in that order. N50 is the length of contigs that just covers 50th percentile. Although longer contigs are better, N50 size deteriorates rapidly as low coverage region in the contigs increases in the attempt to obtain longer contigs. For microbial genome sequencing, the best assembly should only give one contig because the chromosome is circular. But in reality, that never happen.

Recent papers have revealed that paired-end reads are extremely useful in improving assembly and resolve some repeat problems. Read length also helps assembly but until it reaches a barrier according to Chaisson et al (2009). The barrier for E.coli is 35nt while assembly quality for yeast genome doesn't improve much after exceeding 60nt. This comes as a total suprise. The same goes with coverage issues. Assembly doesn't improve after coverage reaches certain threshold. Some repeats and sequence complexity cannot be resolved by high coverage.

And how exactly does what I mentioned has anything to do with SOLiD bioinformatics? Nothing ...because most assemblers mentioned aboved (except Velvet) cannot read colour space data. hehe... I think I have to write another post to come to the points of this SOLiD Bioinformatics Workshop.

Related posts:
SOLiD Bioinformatics Part I : Introduction
SOLiD Bioinformatics Part III: SOLiD softwares

Reference:
Zerbino & Birney 2008. Velvet: Algorithms for de novo short read assembly using De Bruijn graphs
Chaisson et al 2009. De novo fragment assembly with short mate-paired reads: Does the read length matter?

Read more...

MGRC seminar Part II: Plant Genomics

This talk is presented by Dr. Jane Rogers. She was the Head of Sequencing in Wellcome Trust Sanger Institute before joining The Genome Analysis Centre, JIC. She started her presentation with a brief introduction on genome sequencing methods and the types of genome sequences based on sequence quality. The sequence quality measure based on human genome project is showed below:


The main problem in genome sequencing is closing the gaps between large contigs which is thought to be caused by large repeats. According to Dr. Rogers, there are about 340 gaps in the human genome project completed in 2003. Recent studies showed that these gaps are caused by copy number variations (CNV) and different clones that are used for sequencing. This is a fascinating discovery!

Out of nearly 50 plant genome sequencing projects, only complete Arabidopsis and rice genome are obtained. The rest produce draft genome sequences. So why do we need to sequence so many plant genome? The ultimate goal is crop breeding. Genomics assisted selection program is expected to take at least 9 years time before phenotypic selection can be carried out.

Polyploidy plants are the most difficult species to sequence due to huge genome size, repetitive regions and complexity. One example is hexaploid wheat. 80% of this genome is repetitive. Physical mapping of wheat genome using BAC clones has been initiated in 2005. Now Dr. Rogers and her team are interested in sequencing chromosome 3DL of wheat.

Most sequencing projects involve plants with small genome size. Besides Arabidopsis, poplar and Medicago trunculata is used as model plant for tree and legume. Medicago trunculata is a nitrogen fixing plant. The original size of this plant is thought to be 200 Mb but may reach up to 300 Mb. Almost 2600 BAC clones have been sequenced up until Oct 2008. Additional shotgun sequencing using 454 and Illumina platforms are used. The scientists are interested to study nodulation pathway of this plant using genomic approach.

Another important application of genomics is to understand plant resistance. She elaborated on the recent work of Prof. Jonathan Jones which is highly cited. His team sequenced two strains of Arabidopsis white rust using Illumina platform. The 36 bp short reads are assembled using Velvet. Comparison of these two genomes has successfully identified candidate effectors in host resistance mechanism.

Read more...

SOLiD bioinformatics workshop Part I

I almost missed this workshop because I was told by my colleague that it’s just another seminar on SOLiD. Luckily I emailed the sales rep and received the agenda the night before. I think Analisa Resources Sdn. Bhd did a great job organizing this workshop. Not to mention the wonderful buffet lunch at Hotel Equatorial. :p

Illumina Solexa and SOLiD are two NGS platforms that came up almost the same time. Despite similar read length, SOLiD can generate twice the amount of sequence at a slightly lower cost. Judging from the number of publications, it is assumed that most people favour Solexa over SOLiD. But, the scenario is changing.

SOLiD differs from other platforms in which the reads are mapped in colour space. Its 2-base encoding system allows better SNPs detection by differentiating a SNP from sequencing error. Two colour change will indicate a SNPs while single colour change indicates an error.

Despite what paper reviews said about SOLiD, here’s the things they can’t tell you
• Since Invitrogen has acquired Applied Biosystems (not the other way around) last year, Applied Biosystems now provide end-to-end solutions from sample preparation to sequencing.
• SOLiD can run two slides per run. Each slide can be divided up to 8 chambers. With its barcoding system, up to 320 samples can be run each time (with 20 barcodes).
• Both 454 and SOLiD uses emulsion PCR concept but the difference is SOLiD places the bead randomly on the flow cell (That’s how they claim the throughout is increased). But now, they are looking into ways to increase bead density by arranging beads on the slide.
• The new SOLiD version 3 allows users to view the sequencing process via the internet and control the process remotely.
• SOLiD will give forward and reverse reads (This is new to me. I wonder if other platforms can do the same)
• Besides higher accuracy 99.96%, SOLiD offers different size of insert size up to 4 kb to enable scaffolding. Scaffolding is important to form larger contigs.
• Unlike other sequencing platform, SOLiD user can download data analysis softwares from the website and sample data. http://www.solidsoftwaretools.com

Related posts:
SOLiD Bioinformatics Part II: de novo assembly
SOLiD Bioinformatics Part III: SOLiD softwares

Read more...

The Next Generation Sequencing Status in Malaysia

Since the introduction of Next Generation Sequencing (NGS) in 2005, the science community in Malaysia has been increasingly aware of it due to the rapid development of NGS technology in various applications. Here are the major local players in Malaysia’s NGS field to look out for:

Synamatix Sdn. Bhd. It is the most experienced company in NGS bioinformatics. It is founded in 2002. Synamatix has successfully established itself overseas with the introduction of SXoligoSearch. This software will align Illumina reads against a range of Refseq RNA or NCBI genome builds for a number of organisms. Synamatix also provides outsourcing NGS service. MGRC, a wholly-owned subsidiary of Synamatix Sdn Bhd. has been actively engaging in talks with potential NGS users and conducting NGS-related seminar and workshop.

Novocraft Technologies Sdn. Bhd. It is another Malaysia-based software company that has gained BioNexus status recently. Its product Novoalign is a fast and efficient tool for Illumina Solexa reads. Unlike SXoligoSearch, Novoalign is free for academic and non-profit use. The company is the official vendor in Malaysia for Cofactor Genomics. Cofactor Genomics is a US-based company that offers Illumina Solexa Seqeuncing. Besides cheaper and better quality of sequencing, it offers transparency in sequencing price and sample queue to its customers.

Two other companies are Science Vision Sdn. Bhd. and Analisa Resources Sdn. Bhd., the official representative for Illumina and Applied Biosystems respectively. Both have something in common. Besides organizing NGS-related seminar and workshop and offering sequencing service via outsourcing, they work closely with their Singaporean counterparts.

After so much talk and expectation, the first NGS platform has finally landed in Malaysia! UiTM has bought one unit of Illumina GAII. I'm sure the local researchers will be looking forward to cheaper Solexa sequencing from UiTM.

Read more...