Next Generation Sequencing for Microbiome Analysis - Getting Started

There are a number of online pipelines for the analysis of ribosomal samples for carrying out microbiome studies. Two of the most frequently used ones are:

• RDP
• MG-Rast

Another possibility is using a pipeline installed on a local machine. The weakness of this approach is that you need to keep the r-RNA databases up-to-date to make sure you get the best results. Possible locally installed pipelines are Pangaea and Qiime (pronounced chime).

An initial challenge can be that the files you get from Next Generation sequencing might not be compatible with the pipelines either because they have a different quality measure to that expected, or because the pipelines cannot work with fastq files.

Transforming Fastq Files

One tool for transforming fastq files is the FASTX-Toolkit. Another simpler and faster method but with much less functionality is to use BioPython (you need to install it first but that is relatively simple on Linux systems like Ubuntu). For the FASTX toolkit if you are using Ion Torrent data do not forget to include the -Q 33 flag to show that a quality string different to the default is being used.

The BioPython code for creating the fasta files from the fastq files is:

SeqIO.convert("filename.fastq", "fastq", "output.fasta", "fasta")

To create the accompanying quality file the command is:
SeqIO.convert("filename.fastq", "fastq", "output.qual", "qual")

R in Ubuntu

Keeping up to Date

One of the most frustrating things about using distributions and not installing from source is making sure that the programs are up to date. This is especially true of a rapidly changing program such as R and Bioconductor. So the first thing that is needed is to add the repository to those searched for updates.

You can do this by adding the following to the /etc/apt/sources.list file using your favourite editor (do not forget to use sudo)

deb http://cran.ma.imperial.ac.uk/bin/linux/ubuntu precise/

deb-src http://cran.ma.imperial.ac.uk/bin/linux/ubuntu precise/

For R to install the rgl library properly you also need to make sure that you have all of the headers installed (gl.h and glu.h). You can do this by installing the following libraries.

sudo apt-get install build-essential libsdl1.2debian libsdl1.2-dev libgl1-mesa-dev libglu1-mesa-dev libsdl-image1.2 libsdl-image1.2-dev

When you install the packages it is best if you install as root as some core packages need to be updated and so run R as sudo

sudo R

install.packages("multicore")
install.packages("sem")
install.packages("leaps")

install.packages("e1071")

install.packages("aplpack")

Bioconductor

You need curl and xml to connect to databases.

sudo apt-get install libcurl4-gnutls-dev

sudo apt-get install libxml++2.6-dev

source("http://bioconductor.org/biocLite.R")

biocLite()

biocLite("ShortRead")

biocLite("Rsamtools")

biocLite("nucleR")

biocLite("SRAdb")

biocLite("edgeR")

Annotating Next Generation Sequencing

To annotate a vcf file using vcftools you first need to create a file containing the annotations.

This should have the following format:


#CHR FROM TO Annotation
5 53719508 53936990 ENSGALG00000011590; gene_name=TMEM179; gene_type=KNOWN_protein_coding
5 54024011 54038102 ENSGALG00000011608; gene_name=INF2; gene_type=KNOWN_BY_PROJECTION_protein_coding
5 54073641 54096083 ENSGALG00000011618; gene_name=ADSSL1; gene_type=KNOWN_protein_coding


I got this file from editing a GFF file for the region of interest. This file then needs to be zipped with bgzip and indexed with tabix. You will need a different annotation file for each feature if you edit a gff as otherwise they will not appear in order in the file.


bgzip annotate.gff
tabix -p gff annotate.gff.gz


The the annotations can be added to the file.


cat input.vcf.gz | vcf-annotate -a annotate.gff.gz -d key=INFO,ID=ANN,Number=1,Type=Integer,Description='My custom annotation' -c CHROM,FROM,TO,INFO/ANN > out.vcf


Once the file is annotated you can use grep to pull out the lines with the required annotation.

Next Generation Sequencing

The problem is combining variant calls from different species to a reference genome in order to find the variants between the two species.

First you need to edit the filtered vcf files to remove the sample name - as the program is intended to find differences between samples of the same species.


vi varX.flt.vcf
vi varY.flt.vcf


Next you need to use bgzip and tabix from Heng Li to get a compressed and indexed datafile.


bgzip varX.flt.vcf
bgzip varY.flt.vcf
tabix -p vcf varX.flt.vcf.gz
tabix -p vcf varY.flt.vcf.gz


Next you can use the vcftools function vcf-isec to find the complements of the two datasets. These will be the variants that are unique to the different species.


vcf-isec -c varX.flt.vcf.gz varY.flt.vcf.gz | bgzip -c > unique_varX.vcf.gz
vcf-isec -c varY.flt.vcf.gz varX.flt.vcf.gz | bgzip -c > unique_varY.vcf.gz


You can also create a Venn diagram of the overlap of variants between the different species.


vcf-compare var0.flt.vcf.gz var15.flt.vcf.gz > venn.out


And also look at the overlap in variants


vcf-isec -o -n +2 var15.flt.vcf.gz var0.flt.vcf.gz | bgzip -c > overlap_var15.vcf.gz


Permutating miRBase

You can also user seqinr to permutate sequences.

For the miRBase example the code is:

permutation(mature[[1]],modele='base',frame=0)


GC Content Analysis of miRBase

I want to create a dataset for supervised learning based on miRBase. So I need to know the statistical properties of miRBase. Having used seqinr to import miRBase into R, I need to carry out some analysis. GC content is a first place to start as many of the sequences I am familiar with have a high AT content and so an AT bias might be a factor affecting the learning process.

There are 19724 sequences in the database. I could have used length(mature) instead of hard coding the number. You also need to initialise out as a variable first.

> for (x in c(1:19724)) out <- c(out, GC(mature[[x]])) > hist (out)


So I am happy with this as it is normally distributed and so I do not need to take any special care in the GC/AT content for my training sets.

Plan to automate miRbase

Need to check for a new release of miRbase and then download it

It is located at:
ftp://mirbase.org/pub/mirbase/

Need to check the data-stamp on the directory that CURRENT points to and then download it if it is newer than the date-stamp of the last update on the local machine. On average 6 monthly.

Need to get the mature.fa.gz file from this directory

Then this needs to be moved to the local /usr/share/ncbi/data directory
Then it needs to be gunziped
Then it needs to be converted to a blast database with the formatdb command.

Grep gga from mature.fa > gga_mature.fa selects out those from Gallus gallus.

To use
megablast -d /usr/share/ncbi/data/mature.fa -i sequences.fa -o sequences.out -W 7 -D 2 -p 100 [-e 1000]