A new peak calling software: peakzilla

From Github https://github.com/steinmann/peakzilla

PEAKZILLA
---------
Peakzilla identifies sites of enrichment and transcription factor binding sites from transcription factor ChIP-seq and ChIP-exo experiments at hight accuracy and resolution. It is designed to perform equally well for data from any species. All necessary parameters are estimated from the data. Peakzilla is suitable for both single and paired end data from any sequencing platform.

Note that peakzilla is not suited for the identification of broad regions of enrichment (e.g. ChIP-seq for histone marks), we recommand using MACS instead: Zhang et al. Model-based Analysis of ChIP-Seq (MACS). Genome Biol (2008) 9(9):R137

web/GUI tools for biological data analysis

If you are not a programmer, you will be happy to look at these websites. Even you can program, I find it is sometimes convenient to use them :)

1. Galaxy https://main.g2.bx.psu.edu/
founded in PSU, it's been there for several years, and under active development.

2. GenePattern from Broad Institute http://www.broadinstitute.org/cancer/software/genepattern
I used it sometime ago for some microarray data analysis  when I did not anything about programming. Now, I use R/Bioconductor http://www.bioconductor.org/

3. Taverna http://www.taverna.org.uk/
I got to know it from a blog, never used it, but it looks promising

4. GenomeSpace http://www.genomespace.org/  it integrates many tools together including the above Galaxy and GenePattern

There are just way too many softwares out there for different kind of data analysis. First, get to know what are out there. Then, choose the right one meet your need for research.

Two nature protocols for RNA-seq analysis

One is from Simon Anders, his group wrote the DESeq

http://www.nature.com/nprot/journal/v8/n9/full/nprot.2013.099.html

Count-based differential expression analysis of RNA sequencing data using R and Bioconductor

The other one 

http://www.nature.com/nprot/journal/v7/n3/full/nprot.2012.016.html

Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks

I will go through a complete RNA-seq analysis following the guidelines.

There are other packages make RNA-seq analysis easier:

Blacktie https://pypi.python.org/pypi/blacktie/0.2.1

https://github.com/BradyLab/RNASeq

Nesoni http://www.vicbioinformatics.com/nesoni-cookbook/

Lots of things are waiting for me to try out!

How to make a heatmap based on ChIP-seq data by R

update on 04/20/2016.
I noticed that this post is the most frequently visited one. it has been almost 3 years
since I wrote this post. Now, there are various tools for this purpose. See a
post on biostars.

well, I recently just went through the whole process for making a heatmap based on a ChIP-seq data set. If you do not know the technique, google it :)  http://en.wikipedia.org/wiki/ChIP-sequencing

Often, you have a ChIP-seq data that are mapped to the reference genome ( a bam file). You want to plot the sequence tag intensity around certain features ( transcription start sites, gene body, enhancers, or any other genomic region you defined).

you can make an average plot http://crazyhottommy.blogspot.com/2013/04/how-to-make-tss-plot-using-rna-seq-and.html. I will need to re-write this one though, the code format is just too bad (I have to learn how to embed R and python code into the blog) ....and the Y axis is not normalized to counts per million.

you may also want to generate a heatmap with the same data.  see ngsplot for examples https://code.google.com/p/ngsplot/.  If you do not want to code R by yourself, try it. It has been improved a lot since last time I checked it. I once asked a question in the google group : https://groups.google.com/forum/#!topic/ngsplot-discuss/efHQ-P-14XM.

Seqmonk http://www.bioinformatics.babraham.ac.uk/projects/seqmonk/  from Simon Andrews can also plot this kind of figure very easily, I am just not satisfied with the picture quality, and I want more customized control of the picture.


I will just paste my code below, and it is heavily commented, you should be able to follow it fairly easily.
update on 05/05/2015, I put the code in a gist instead:

	# This R script is to generate the TF or histone modification heatmap
	# at certain genomic features (TSS, enhancers) from the ChIP-seq data
	# the input matrix is got from Homer software. alternative to R, use cluster3 to cluster, and visualize by # java Treeviewer
	# generate the matrix by Homer: annotatePeaks.pl peak_file.txt hg19 -size 6000 -hist 10 -ghist -d TF1/ # > outputfile_matrix.txt
	# see several posts for heatmap:
	# http://davetang.org/muse/2010/12/06/making-a-heatmap-with-r/
	# http://www.r-bloggers.com/r-using-rcolorbrewer-to-colour-your-figures-in-r/
	# 08/20/13 by Tommy Tang
	# it is such a simple script but took me several days to get it work...I mean the desired
	# color setting of the heatmap. Making a pretty figure takes time :)
	# use either heatmap.2 from gplots or pheatmap
	library(gplots)
	library(pheatmap)
	# you can compare the heatmaps generated by heatmap.2 and pheatmap
	setwd("/home/tommy/homer")
	d1 <- read.table("outputfile_matrix.txt", header=T)
	# have a look at the matrix
	head(d1)
	d1$Gene # the dataframe generated by homer has a column named Gene
	m1<- as.matrix( d1[,2:ncol(d1)]) # the first column is the TSS id,
	rownames(m1)<- d1$Gene # heatmap.2 works only on matrix, turn the dataframe to matrix, and # add the TSS id as the row name
	m1<- log2(m1+1) # log2 transform the raw counts
	#exam the data range
	range(m1)
	hist(m1) # histogram to look at distributions
	# it is from the ChIP-seq count data, many of them are 0s. I transfromed the raw data
	# by log2(m+1), so if the log2 value is 0, the raw number is also 0 (count)
	# to compare different heatmaps, I have to map the value to the color using the
	# break argument in pheatmap or heatmap.2 That's why I examed the data range.
	# http://stackoverflow.com/questions/17820143/how-to-change-heatmap-2-color-range-in-r
	# https://stat.ethz.ch/pipermail/bioconductor/2011-November/041866.html
	# http://seqanswers.com/forums/showthread.php?p=114275&posted=1#post114275
	bk = unique(c(seq(-0.1,3, length=100),seq(3,9.7,length=100)))
	hmcols<- colorRampPalette(c("white","red"))(length(bk)-1)
	# you can play around with the break points, 9.7 is the max of the matrix m1
	# if just use bk = c(seq(-0.1,1, length=100),seq(1,12.5,length=100)) without the unique function
	# the pheatmap gave me error message:
	# Error in cut.default(x, breaks = breaks, include.lowest = T) :
	# 'breaks' are not unique
	# It is caused by concatenating several seq() results together, which share the same boundaries.
	# try ?dist ?hclust
	png("test.png", width=300, height = 800) # width and height are in pixel
	heatmap.2(m1, col=hmcols, breaks = bk, Rowv= TRUE , Colv=FALSE, dendrogram="row", useRaster = TRUE, symkey=FALSE, symm=F, symbreaks=T, scale="none", trace="none", labRow=NA, labCol=NA)
	# do not show the row and column labels
	# trace argument should be put to "none", otherwise the trace is cyan and it will "eat" the heatmap
	# I asked in Biostar http://www.biostars.org/p/79444/#79457
	# look at the documentation for #heatmap.2 http://hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/gplots/html/heatmap.2.html
	# or you can use pheatmap
	pheatmap(m1, color=hmcols, breaks= bk, cluster_rows=TRUE, cluster_cols=FALSE, legend=FALSE, show_rownames=FALSE, show_colnames=FALSE)
	dev.off()

view raw heatmap_ChIP-seq.r hosted with ❤ by GitHub

The second gist:

	#The X axis is -3 kb to 3 kb around TSS.
	#It turns out that the heatmap.2 function use default hclust ( Hierachical Clustering) to cluster the #matrix.
	#alternatively, we can use K-means clustering to cluster the data and to see what's the pattern look like.
	km<- kmeans(m1,2) # determine how many cluster you want, I specify 2 here
	m1.kmeans<- cbind(m1, km$cluster) # combine the cluster with the matrix
	dim(m1.kmeans)
	# [1] 903 602
	# the last column is 602
	o<- order(m1.kmeans[,602]) # order the last column
	m1.keans<- m1.kmeans[o,] # order the matrix according to the order of the last column
	pheatmap( m1.kmeans[,1:601], cluster_rows = F, cluster_cols = F, col= hmcols, breaks = bk, legend=FALSE, show_rownames=FALSE, show_colnames=FALSE)
	# no need to cluster by pheatmap, so I set cluster_rows and cluster_cols to FALSE

view raw heatmap_ChIP-seq_2.r hosted with ❤ by GitHub

That's all!

I hope you have learned something after reading it:)
===============================
update on 09/17/13
arrange the rows in the heatmap by the coverage from strong to weak

	m.row.sum<- cbind(m1, rowSums(m1))
	o1<- rev(order(m.row.sum[,602]))
	m.row.sum<- m.row.sum[o1,]
	bk = unique(c(seq(-0.1,3, length=100),seq(3,10.35,length=100)))
	hmcols<- colorRampPalette(c("white","red"))(length(bk)-1)
	pheatmap( m.row.sum[,1:601], cluster_rows = F, cluster_cols = F, col= hmcols, breaks = bk, legend=FALSE, show_rownames=FALSE, show_colnames=FALSE)

view raw heatmap_ChIP-seq_3.r hosted with ❤ by GitHub

install GIMP 2.8 in Ubuntu 12.04

I was trying to make some figures, and google told me GIMP is a very popular software as an alternative to Photoshop.

Ubuntu 12.04 has GIMP 2.6 installed. To update to GIMP 2.8, follow the instructions here:
http://www.maketecheasier.com/upgrade-to-gimp-2-8-in-ubuntu/2012/05/04

and install the plugins:
http://www.webupd8.org/2013/02/gimp-284-released-install-it-in-ubuntu.html

It takes long time to make a good-looking figure.
I will write a post to demonstrate how to make a heatmap by R next time.

pdf organizer and reference manager

After being in the research field for a while, I found myself lost in a vast amount of pdf research papers.
I need a good pdf organizer.
especially when I start writing a paper, putting the reference  become a headache to me.
That's where Zotero and Mendeley come to help.
A comparison for these two:
https://www.youtube.com/watch?v=53QRA_SNSzE

http://www.zotero.org/
Zotero is open source, it  can organize other files like images and videos besides pdf.

http://www.mendeley.com
Mendeley is not open-source, it has a better support for pdf files

watch youtube videos on how to use it:
https://www.youtube.com/watch?v=xLtk6n8cFdk
https://www.youtube.com/watch?v=kSu4tO8Dwak
https://www.youtube.com/watch?v=VD1z0boSpQY
https://www.youtube.com/watch?v=XKAcQTea7E8
https://www.youtube.com/watch?v=zkrVbBSrK_w

I am currently using Mendeley.

Decorators in python

To help you understand Decorators in python, read the posts below:

http://simeonfranklin.com/blog/2012/jul/1/python-decorators-in-12-steps/

http://www.artima.com/weblogs/viewpost.jsp?thread=240808

and here:
http://pythonarticles.com/decorators_part1.html
http://pythonarticles.com/decorators_part2.html

Enjoy!

Two pipelines for NGS analysis

bcbio-nextgen from Brad Chapman

https://github.com/chapmanb/bcbio-nextgen
Documentation https://bcbio-nextgen.readthedocs.org/en/latest/index.html

Metaseq from Ryan Dale

https://github.com/daler/metaseq

Documentation http://pythonhosted.org/metaseq/

Very good R tutorial by Google developer

watch the videos here:
https://www.youtube.com/watch?v=iffR3fWv4xw&list=PLOU2XLYxmsIK9qQfztXeybpHvru-TrqAP

nice demonstration of *args and **kwargs in python function def

see a post here
http://freepythontips.wordpress.com/2013/08/04/args-and-kwargs-in-python-explained/

python for biologist

I just found this website:
http://pythonforbiologists.com/index.php/welcome-to-python-for-biologists/

it teaches you how to start programming with python for biologists.  I even bought the pdf ebook for $39. I had a look at the book, and I think it is worthy. I went through first 4 chapters really fast, they are really basic but teach biologists in a really friendly way. I am even thinking to hold a bioinformatics training course in the future :)  sounds ambitious, hmm..

Advanced topics can be found in the website:

Here is a link to a .zip file containing solutions to some of the exercises.

Section 1 : Large programs and complex data structures

Section 2 : Retrieving and sending information over the web

Section 3 : Large datasets and treelike data structures

Section 4 : Classes, modules and web applications

I've read the practical computing for biologists by Haddock Dunn http://www.amazon.com/Practical-Computing-Biologists-Steven-Haddock/dp/0878933913
it is a very good introductory book. Besides python, it also teaches you how to use linux command lines and some database basics. By the way, I wrote my first regular expression after reading the book.

well, I feel I just could get some practical problems ( text manipulation, Genomic Interval calculation using pybedtools, NGS by HTSeq) done by python though I began to play around with python in 04/24/2012. It takes time, but the investment definitely is rewarding!!

UCSC genome browser tutorial

see here openhelix  http://www.openhelix.com/cgi/tutorialInfo.cgi?id=76

I use UCSC genome browser so much, but I did not realize that it has much more functions I've missed before watching the video.

Related tutorials

This tutorial is a part of the tutorial group UCSC Tutorials. You might find the other tutorials in the group interesting:

• ENCODE Foundations: ENCyclopedia of DNA Elements
• UCSC Genome Browser: An Introduction: The UCSC Genome Browser Introduction
• UCSC Genome Browser: Custom Tracks and Table Browser: UCSC Genome Browser advanced topics
• UCSC Archaeal Genome Browser: Provides you with many research and analysis tools that can be used to examine the genomes of more than 50 microbial species from the domain archaea.
• ENCODE Data at UCSC: ENCODE Data at UCSC