Final Stages of Using Qiime for Ion Torrent Data

Once you have the rarefaction curves for within the sample you can now calculate them for between the bar-coded samples. This will give a measure of distance between the samples that can then be used for clustering.

beta_diversity_through_plots.py -i otus/otu_table.biom -m mapping.txt -o wf_bdiv_even4/ -t otus/rep_set.tre -e 4

Where the last term is a cutoff for the minimum number of reads in a sample (must be even) - this is the sampling depth. You need to try a series of cut-offs to find out which works best, but keep it less than the median. If you choose the median then values below the median are ignored (you will lose half of you data but if you are expecting a large number of blank results this might be OK) so only do this if you have lots of skew and the mean is much larger than the median.

This will give a series of 2D and 3D Principle Components plots that can be used to look for any clustering of the samples. Clusters should be visually identifiable and well defined.

Qiime can also carry out a jack-knife analysis to test the clustering dependence on omitting some of the samples. This produces a phylogenetic tree where like samples are on closely related nodes. Clusters are represented by clades (a group of leaves with a common root).

jackknifed_beta_diversity.py -i otus/otu_table.biom -t otus/rep_set.tre -m mapping.txt -o wf_jack50 -e 50

Finally you can make Biplots to show the most important axes between the clusters.

make_3d_plots.py -i wf_bdiv_even4/unweighted_unifrac_pc.txt -m mapping.txt -t wf_taxa_summary/otu_table_L3.txt --n_taxa_keep 5 -o 3d_biplot

The 5 tells the program how many top level taxa to display for comparison to the clusters.

Using Qiime for Microbiome Analysis of Ion Torrent Data

First read my last post about getting Fasta and quality files from the Fastq files. If there are multiple samples in the run then there will be a single file containing all of the data with the different bar-codes attached.

The process to follow is the same as for 454 data.

You sample should be made up of read containing:

1. An adaptor (A) sequence
2. A barcode sequence
3. A linker primer sequence 
4. The actual target sequence
5. A reverse primer sequence
6. An adaptor (B) sequence

First you need to create a mapping file telling Qiime about the barcodes, linkers adaptors and any experimental conditions you might want to add.

• #SampleID BarcodeSequence LinkerPrimerSequence Treatment ReversePrimer Description
• 1 CTAAGGTAAC CCTACGGGAGGCAGCAG control ATTACCGCGGCTGCTGG oral swab
• 
  

This file is easiest to create in a spreadsheet and then save it as a tab delimited text file. You should then check this file to make sure that it is in the right format.

check_id_map.py -m Mapping.txt -o mapping_output

If you are doing a properly blinded experiment then you should not know what the treatments are and so this should be just given a number or a simple code.

The first stage of the Qiime pipeline is splitting the sample by barcodes removing any low quality or ambiguous reads. 

split_libraries.py -M3 -l12 -L250 -b variable_length -m Mapping.txt -f Example.fna -q Example.qual -o split_library_output

The command line arguments allow a mismatch in the primer of 3 a minimum sequence length of 12 and maximum of 250 and allow for variable barcode lengths. This minimum length is too short and for a real run with quality control it should be > 120.

It is important to check the output to make sure you are getting the number of sequences you expect from the data. The output is in the directory split_library_output.

Look at the histograms.txt file. If you get something like this then you have a problem.

# bins raw sequence lengths, length of sequences that pass quality filters before processing, and lengths of sequences that pass quality filters post processing.
Length Raw Before After
0 4602 0 28
10 130658 0 58
20 15425 10 55
30 46920 32 46
40 118695 58 28
50 95443 57 43
60 67658 39 102
70 65570 27 124
80 78656 54 99
90 101635 117 82
100 94073 114 108
110 99076 100 80
120 93229 76 116
130 88529 114 145
140 79631 94 506
150 79257 133 1991
160 114626 124 539
170 408446 1709 3027
180 768141 810 27
190 293266 1429 1
200 1051049 2103 1
210 31571 5 2
220 2318 0 0
230 1175 1 0
240 623 2 0
250 534 0 0
260 309 0 0
270 315 0 0
280 213 0 0
290 136 0 0
300 75 0 0
310 47 0 0
320 32 0 0
330 45 0 0
340 29 0 0
350 33 0 0
360 51 0 0
370 5 0 0

There might be a problem because the quality control is removing almost all of your data before it puts it into the different splits and removes the barcode sequences. In this case there is a lot of junk sequence because of short reads witjout primers. Looking at the top of the split_library_log.txt file makes this clear

Number raw input seqs 3932096
Length outside bounds of 12 and 250 26356
Num ambiguous bases exceeds limit of 6 0
Missing Qual Score 0
Mean qual score below minimum of 25 1222172
Max homopolymer run exceeds limit of 6 7327
Num mismatches in primer exceeds limit of 3: 2491654

Next Qiime will cluster the sequences together based on a threshold of 97% identity. It will then take a representative member of each cluster and compare it to the taxonomic database to create a list of Operational Taxonomic Units (OTUs). Qiime carries out all of the necessary steps including alignments and creating a clustering tree of the OTUs using the following command:

pick_de_novo_otus.py -i split_library_output/seqs.fna -o otus

You can see a summary of the output with the following command:

print_biom_table_summary.py -i otus/otu_table.biom

This is an example output for a sample with a small number of sequences

Num samples: 45
Num observations: 1313
Total count: 7180.0
Table density (fraction of non-zero values): 0.0461
Table md5 (unzipped): c3d9bcf2d50c7477f7eadd115f5a0b99

Counts/sample summary:
 Min: 3.0
 Max: 1831.0
 Median: 53.0
 Mean: 159.555555556
 Std. dev.: 334.228834154

The OTUs can be visually displayed as a heat-map showing their relative abundance.

make_otu_heatmap_html.py -i otus/otu_table.biom -o otus/OTU_Heatmap/

This produces a webpage that can be hovered over or clicked on to provide more details about the samples.

Alternatively you can create a network for visualisation that can be viewed using Cytoscape (Networks are useful in some contexts but the relatedness of nodes here is unclear and so I prefer not to use this form of visualisation).

make_otu_network.py -m mapping.txt -i otus/otu_table.biom -o otus/OTU_Network

The OTUs are then used to cluster the samples. Samples which have the same OTUs present should be clustered together.

summarize_taxa_through_plots.py -i otus/otu_table.biom -o wf_taxa_summary -m mapping.txt

These can then be viewed as area charts under the wf_taxa_summary/taxa_summary_plots folder. The area charts are a bit of an awkward representation as they are shown as continuous plots and so the bar charts which are discrete should be used.

Ecologists calculate the within sample (alpha) and between samples (beta) diversities, to show the range of species found in the samples. If you want to use the Shannon Index you have to tell Qiime that you want to calculate it using the following command.

echo "alpha_diversity:metrics shannon,PD_whole_tree,chao1,observed_species" > alpha_params.txt

The diversity and rarefaction are then calculated using the following script.

alpha_rarefaction.py -i otus/otu_table.biom -m mapping.txt -o wf_arare/ -p alpha_params.txt -t otus/rep_set.tre

The curve is assymptotic to having a complete sampling of all of the species in the sample

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]

Useful tools for linux

ps -eo pid,pmem,cmd

Shows the memory usage of the processes that are currently running.

updatedb

Updates the locate database.

locate String

Very fast find routine.

MySQL

If the tar contains the binaries


groupadd mysql
useradd -r -g mysql mysql
cd /usr/local
ln -s /usr/local/src/msyql-VERSION mysql
cd mysql
chown -R mysql .
chgro -R mysql .
scripts/mysql_install.db --user=mysql
chown -R root .
chown -R mysql data
bin/mysqld_safe --user=mysql &
.bin/mysql_secure_installation
cp /usr/local/mysql/support-files/mysql.server /etc/init.d/mysql.server


Postgresql


./configure --prefix=/usr/local/pgsql
gmake
gmake install
adduser postgres
mkdir /usr/local/pgql/data
su - postgres
/usr/local/pgsql/bin/initdb -D /usr/local/pgsql/data
/usr/local/pgsql/bin/pg_ctl -D /usr/local/pgsql/data -l logfile start


Testing

/usr/local/pgsql/bin/createdb test
/usr/local/pgsql/bin/psql test


Configuring Apache/PHP

./configure --prefix=/usr/local/apache --enable-so
make
make install
vi /usr/local/apache/conf/httpd.conf
Add the places from the webpages - usually this will not actually be in the same place as you are running apache and it is most commonly /var/www/html

/usr/local/apache/bin/apachectl -k start

Upgrading

copy the config.nice file from the working version of the server to the new version directory. Use a new prefix and a different port to test the configuration - change the Listen directive to specify the new port.
./config.nice
make
make install
usr/local/apache/bin/apachectl -k graceful-stop
usr/local/apache/bin/apachectl -k start
add a call to apachectl to the rc.local file so that it autostarts.

PHP

./configure --prefix=/usr/local/php \
--with-apxs2=/usr/local/apache/bin/apxs \
--with-config-file-path=/usr/local/php \
--with-gd \
--with-mysql=/usr/local/mysql/ \
--with-mysqli=mysqlnd \
--with-pdo-mysql=mysqlnd \
--with-pgsql=/usr/local/pgsql

make
/usr/local/apache/bin/apachectl stop
make install

Edit httpd.conf to add
LoadModule php5_module modules/libphp5.so
AddType application/x-httpd-php .php
<Files *.php>
SetOutputFilter PHP
SetInputFilter PHP
</Files>


Edit php.ini for pgsql and mysql lines removes the ; but not on the dll lines.

Setting up MySQL

When you add MySQL to your Linux (in my case Fedora) system it has a root account but no password! This is not a very good idea and so the first thing that you need to do is create a root account and password. You do this using the mysqladmin tool as root.

# mysqladmin -u root password new_password


Where new_password is the new password (not just the word new_password).


You can now login to the MySQL using the root username and password and create a new account. 


# mysql -u root -p
...
mysql> CREATE USER 'drupal'@'localhost' IDENTIFIED BY 'password';


You need to the exit (use quit) and use mysqladmin again to create the database.
# mysqladmin -u root -p create drupalbase
# mysql -u root -p
...
mysql> GRANT SELECT, INSERT, UPDATE, DELETE, CREATE, DROP, INDEX, ALTER, LOCK TABLES, CREATE TEMPORARY TABLES ON drupalbase.* TO 'drupal'@'localhost' IDENTIFIED BY 'password';

Dante-R

Dante-R is a program for normalising mass spectroscopy data. It is downloaded as an executable installer from the pacific northwest national laboratory.

Note that this is not an R zip file that you have to manage from within R - this is an executable. The one tricky point with installation is making sure that R is in your path. In my case on Windows 7 the path was not amended to contain R and so I had to add it manually.

To do this I added the following:
;C:\Program Files (x86)\R\R-2.10.1\bin

You get to the path from the Control Panel.

Once the path is added the Danter-R installer will then install all of the required R packages - in particular GTK support.

Dante-R has been used for unlabelled mass-spec data where there are a large number of samples, in cases where there are fewer samples as in the case of processing labelled data the methods might not be robust [Private discussion with the author] and so this has to be checked.

Annotation

One of the best features of Bioconductor is the built in annotation of microarray probes. This makes it veryeasy to query a large number of databases to test out any potential biological explanations for a particular expression profile.

The following code snippet will take the top values above a p-value threshold for the corrected multiple testing (if Benjamini-Hochberg is used this will be the false discrovery rate) and produce an HTML file that summarises the significant results.

library("KEGG.db")
library("GO.db")
library("annaffy")
library("XML")
library("annotate")
library("hgu133a.db")

afarms1<-topTable(fit2farms, coef=1, adjust="BH", n=100, p.value=0.01)
gnfarms <-as.character(afarms1$ID)
symfa1 <-getSYMBOL(gnfarms, "hgu133a")
symfa1
atest <- aafTableAnn(gnfarms, "hgu133a.db", aaf.handler())
saveHTML(atest, file="reportfarms1.html")

Linear Models and Testing

Perhaps the best reason for using R for building complex statistical models is the strength of the linear model packages, which make building complex models reasonably straight forward, even if you have multiple groups. This is important following cluster analysis where you are unlikely to only end up with two clusters.

The package for constructing linear models from microarray data is Limma.

The following code constructs a model for three different clusters (the 0 cluster is used to label all the arrays that did not fall into one of the three identified clusters). First a design matrix is constructed which will define that comparisons will be carried out. The basic model is then fitted, then a second model with more complex interactions is designed using the contrast matrix. This model is then fitted using empirical Bayes (this is most effective at small sample size).

library(limma)
f <- factor(f2LCgcrma1$Cluster, levels =c("0", "1", "2", "3"))
design <-model.matrix(~0 + f)
colnames(design) <- c("cl0", "cl1", "cl2", "cl3")
fitgcrma <- lmFit(f2LCgcrma1, design)
contrast.matrix <- makeContrasts(cl3-cl2, cl3-cl1, cl2-cl1, levels=design)
fit2gcrma <- contrasts.fit(fitgcrma, contrast.matrix)
fit2gcrma <- eBayes(fit2gcrma)
tab1gcrma<-topTable(fit2gcrma, coef=1, adjust="BH")
tab2gcrma<-topTable(fit2gcrma, coef=2, adjust="BH")
tab3gcrma<-topTable(fit2gcrma, coef=3, adjust="BH")

Finding the Clusters

When is a cluster real or an artefact of the analysis or noise in the data? You want to find consistency of the dendrograms between methods and also for different levels of clustering. Confidence in the clustering falls as there are more contradictions between methods and between filtering cut-offs. Here are some ideas about what you should be looking at:

1. What is the effect of the different normalisation methods? - compare rma with gcrma for example
2. What is the effect of filter cut-off? - try samples of 2000 and 500 genes
3. Is the clustering affected by using different measures? Manhattan vs Euclidean
4. Is the clustering affected by the clustering algorithm? - agglomerative or divisive

Adding the Clustering to the Expression Data

The easiest way to add the clustering to the Expression Data is to edit a separate file which can be used to update the phenoData element of filtered exprtession set. This will be needed for subsequent statistical testing and it is easier to do it this way rather than trying to edit individual phenoData elements. The file has to be a text csv file that can be created in any spreadsheet program. This is then read into R and attached to each of the normalised and filtered data files.


pd<- read.AnnotatedDataFrame("pheno_cluster.csv",sep="\t",quote = "\"'")
phenoData(fLCgcrma1)<-pd

Clustering the Arrays

If you are looking at biological variation you may want to cluster the arrays. When analysing a study of 54 tissue culture arrays from lung cancer cell lines the arrays were clustered after filtering using euclidean distance matrices for the log2 expression levels.


dgcrma1 <- dist(log2(t(exprs(fLCgcrma1))),method="euclidian")
hcgcrma1 <- hclust(dgcrma1, method = "average")
plot(hcgcrma1)

The cluster dendrograms from the different normalisation methods can then be compared, and a concensus devised. In this case RMA performed slightly worse than FARMS and GCRMA as the conserved clusters from the other two methods showed some variation in RMA. The factor that makes the biggest difference to the ordering is the number of genes in the clustering and so dendrograms are very sensitive to the filtering step.