qseq to fastq conversion

I was looking at some old ChIP-seq data with raw file in qseq format. I want to convert them to fastq file for mapping with bowtie.

A quick google I found:
http://www.biostars.org/p/6682/

the qseqtofastq C++ program http://www.dna.bio.keio.ac.jp/~krisp/qseq2fastq/ has to be compiled by scons, and I encountered some compiling problem. I gave it up and used this java program http://sourceforge.net/projects/snpeff/files/qseq2fastq.jar/download

Usage:

zcat myfile.qseq.tgz | java -jar qseqtofastq.jar -phred64 > myfile.fastq

it took me around 1 hour to finish the conversion of a 600MB tgz file on the computing cluster (single cpu 4GB ram).

I did have a problem at the end:

Exception in thread "main" java.lang.RuntimeException: java.lang.ArrayIndexOutOfBoundsException: 8
        at ca.mcgill.mcb.pcingola.Qseq2Fastq.main(Qseq2Fastq.java:51)
Caused by: java.lang.ArrayIndexOutOfBoundsException: 8
        at ca.mcgill.mcb.pcingola.Qseq2Fastq.main(Qseq2Fastq.java:43)

I counted the lines of the fastq file (divided by 4) and the original qseq file  , and they are equal.  So, I went ahead and mapped the fastq file with bowtie.

Several NGS bioinformatics training materials

1.  Next Generation sequencing wiki book http://en.wikibooks.org/wiki/Next_Generation_Sequencing_(NGS)

2. UCDavis training course https://training.bioinformatics.ucdavis.edu/documentation/

3. MSU training course by Titus Brown http://ged.msu.edu/angus/index.html

4. GOBLET http://www.mygoblet.org/training-portal

5.  more http://ged.msu.edu/angus/bioinformatics-courses.html

6. a list maintained by Stephen Turner in UVA http://stephenturner.us/edu.html

7. GVL practial protocols for ChIP-seq, RNA-seq, variant calling etc https://genome.edu.au/wiki/GVL

These will keep me busy for a while...I feel I do not have enough time to learn :)

Unix sort except the first line

Many times I have a file with a header  needs  to be sorted, but I do not want to sort the header.
I saw it on Twitter:

 Retweeted by Ming Tang

razor ‏@razZ0r  Mar 1

todays fav linux oneliner: command | (read -r; printf "%s\n" "$REPLY"; sort) > output sort everything except first line.

REPLY is the default variable for the read command.

one more general use of the "body" function
http://stackoverflow.com/questions/9281449/unix-skip-header-bash-function
http://unix.stackexchange.com/questions/11856/sort-but-keep-header-line-in-the-at-the-top

# print the header (the first line of input)
# and then run the specified command on the body (the rest of the input)
# use it in a pipeline, e.g. ps | body grep somepattern
body() {
    IFS= read -r header
    printf '%s\n' "$header"
    "$@"
}

IFS is the bash Shell Bourne variable: A list of characters that separate fields used by the shell to split text strings.

$#    Stores the number of command-line arguments that 
      were passed to the shell program.
$?    Stores the exit value of the last command that was 
      executed.
$0    Stores the first word of the entered command (the 
      name of the shell program).
$*    Stores all the arguments that were entered on the
      command line ($1 $2 ...).
"$@"  Stores all the arguments that were entered
      on the command line, individually quoted ("$1" "$2" ...).

Or an awk solution awk 'NR == 1; NR > 1 {print $0 | "sort -n"}'
It is very handy to use.

Test:
tommy@tommy-ThinkPad-T420:~$ cat body_sort.txt 
header
4
6
9
7
14
8

tommy@tommy-ThinkPad-T420:~$ cat body_sort.txt |( read -r; printf "%s\n" "$REPLY"; sort -n)

header

4

6

7

8

9

14

interactive visualization packages for R and python

I've been collecting potential packages for (interactive) visualization of Genomic data.
Here are the ones so far I've collected:

python:
Vispy http://vispy.org/gallery.html
pretty plotlib https://github.com/olgabot/prettyplotlib
D3py https://github.com/mikedewar/d3py
Yhat ggplot for python http://blog.yhathq.com/posts/ggplot-for-python.html
galry https://github.com/rossant/galry
Glue http://www.glueviz.org/en/latest/
Seaborn https://github.com/mwaskom/seaborn
Bokeh https://github.com/ContinuumIO/Bokeh

R:
Click me https://github.com/nachocab/clickme
plotly https://plot.ly/

others:
http://selection.datavisualization.ch/
Gpviz http://icbi.at/software/gpviz/gpviz.shtml
Jheatmap http://jheatmap.github.io/jheatmap/
Caleydo http://caleydo.github.io/caleydo-doc/3.1/index.html#!index.md

There are too many out there, I will need to focus on  one or two.

hosting bigwig by dropbox for UCSC visualization

First, you need to get an idea of what is a bigwig file: https://genome.ucsc.edu/goldenPath/help/bigWig.html

It is the binary form of wig file and allows UCSC genome browser to fetch only the data in the current window. I usually get wig file by MACS1.4 peak calling ChIP-seq data.
MACS2  now does not have the -w option any more. https://github.com/taoliu/MACS/

see a discussion in the google group:
https://groups.google.com/forum/#!searchin/macs-announcement/bedgraph$20ucsc$20track/macs-announcement/LBhAtmC-Zho/uZuxU8ZaqdEJ

MACS2 only creates a bedgraph file http://genome.ucsc.edu/goldenPath/help/bedgraph.html
and one can convert the bedgraph to bigwig:
https://github.com/taoliu/MACS/wiki/Build-Signal-Track
https://gist.github.com/taoliu/2469050
http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/

if you have a bam file, you can produce a bedgraph file by bedtools  genomeCoverageBed:
http://bedtools.readthedocs.org/en/latest/content/tools/genomecov.html
http://www.biostars.org/p/64495/#64680

you need to add a track line (--trackline option) for UCSC genome browser.
https://groups.google.com/forum/#!searchin/bedtools-discuss/bedtools$20bedgraph$20ucsc$20track/bedtools-discuss/3CibmlqIdWA/PW_bhgWQfVMJ

I followed the instructions https://genome.ucsc.edu/goldenPath/help/bigWig.html

To create a bigWig track from a wiggle file, follow these steps:

1. Create a wig format file following the directions here. Note that when converting a wig file to a bigWig file, you are limited to one track of data in your input file; you must create a separate wig file for each data track. Note that this is the file that is referred to asinput.wig in step 5 below.
2. Remove any existing 'track' or 'browser' lines from your wig file so that it contains only data.
3. Download the wigToBigWig program from the directory of binary utilities.
4. Use the fetchChromSizes script from the same directory to create the chrom.sizes file for the UCSC database you are working with (e.g. hg19). Note that this is the file that is referred to as chrom.sizes in step 5 below.
5. Create the bigWig file from your wig file using the wigToBigWig utility like so: wigToBigWig input.wig chrom.sizes myBigWig.bw
   (Note that the wigToBigWig program also accepts a gzipped wig input file.)
6. Move the newly created bigWig file (myBigWig.bw) to a http, https, or ftp location.
7. Construct a custom track using a single track line. The most basic version of the track line will look something like this:
   track type=bigWig name="My Big Wig" description="A Graph of Data from My Lab" bigDataUrl=http://myorg.edu/mylab/myBigWig.bw
   Optional values can include:

     autoScale         <on|off>                             # default is on
     alwaysZero        <on|off>                             # default is off
     gridDefault       <on|off>                             # default is off
     maxHeightPixels   <max:default:min>                    # default is 128:128:11
     graphType         <bar|points>                         # default is bar
     viewLimits        <lower:upper>                        # default is range found in data
     viewLimitsMax     <lower:upper>                        # suggested bounds of viewLimits, but not enforced
     yLineMark         <real-value>                         # default is 0.0
     yLineOnOff        <on|off>                             # default is off
     windowingFunction <mean+whiskers|maximum|mean|minimum> # default is maximum, mean+whiskers is recommended
     smoothingWindow   <off|[2-16]>                         # default is off
     transformFunc     <NONE|LOG>                           # default is NONE

   For further information on custom bigWig track settings, see the Track Database Definition Document. For further information on how bigWig settings are used in native Browser tracks, see the Configuring graph-based tracks page.
8. Paste this custom track line into the text box in the custom track management page.

I've encountered a coordinates problem mentioned here https://groups.google.com/forum/#!topic/macs-announcement/yPSPlKdTOwo when I tried to convert wig to bigwig.

the -clip option in the wigToBigWig program seems to resolve the problem, but it still gives warning messages. The best way is to use bedClip program here first http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/

Now, I need to use some public ftp or http to host the resulted bigwig file. I decided to use dropbox as I have used it for a while and I have 15G store space. The UCSC genome browser has problem to accept my link. I found a answer here:
http://bergmanlab.smith.man.ac.uk/?p=1989

"The first problem with the Share Link function is that the URL automatically generated by Dropbox cannot be read by the UCSC Genome Browser. For example, the link generated to the file “test.bed” in my Dropbox folder is “https://www.dropbox.com/s/7sjfbknsqhq6xfw/test.bed”, which gives an “Unrecognized format line 1″ error when pasted into the UCSC Browser.  This can easily be fixed if you just want to load a single custom track  to the UCSC Browser using Dropbox by simply replacing “www.dropbox” in the URL generated by Dropbox with “dl.dropboxusercontent”. In this example, the corrected path to the file would be “https://dl.dropboxusercontent.com/s/7sjfbknsqhq6xfw/test.bed”, which can be loaded by the UCSC Genome Browser automatically."

Finally I had the data visualized in UCSC!
lsd1 binding at the oct4 locus:

Next time, I will get a bedgraph file from the MACS1.4 or MACS2, and use the bdg2bw program https://gist.github.com/taoliu/2469050
" conversion to bedgraph is necessary to reduce the final bw size (up to 70%)"

=========================================
I got an email from dropbox, this was mentioned in the blog http://bergmanlab.smith.man.ac.uk/?p=1989
I guess I need to find another place to host my bigwig files...

Hi Ming, This email is an automated notification from Dropbox that your Public links have been temporarily suspended for generating excessive traffic. Your Dropbox will continue to function normally with the exception of Public links. For more information on suspended links, please visit the Help Center. If this is your first suspension, you may remove the suspension by visiting your account page.

how to get a genome-wide motif bed file

Someone was asking this question on Seqanswers http://seqanswers.com/forums/showthread.php?t=40762&highlight=genome+motif+bed

For motif analysis, the most popular program is MEME http://meme.nbcr.net/meme/
There are a bunch of tools in the suites including some useful ones for ChIP-seq
I also saw RAST http://rsat01.biologie.ens.fr/rsa-tools/index.html and oPOSSUM were mentioned http://opossum.cisreg.ca/oPOSSUM3/

I analyze ChIP-seq data a lot. Usually, one gets a bed file containing the positions of the peaks. However, many motif analysis programs require fasta file as input.

There are many ways to get fasta file based on coordinates.
See:
http://www.biostars.org/p/7481/
and my previous post http://crazyhottommy.blogspot.com/2013/04/batch-converting-coordinates-to.html

Going back to the question, I've found several ways to get a bed file containing the coordinates for the motifs.

1. Homer homepage
http://homer.salk.edu/homer/
at the very bottom, there are several links for human and mouse .The file contains all the known motif coordinates for the whole genome. Files are big (several Gb). one can get the specific motifs occurrences by grep.

2.  from UCSC software http://genome.ucsc.edu/ENCODE/analysisTools.html
ENCODE-motifs  http://compbio.mit.edu/encode-motifs/
at the bottom of the page:

only the human data are available.

3. Motif-map http://motifmap.ics.uci.edu/
click motif search

Several more organisms are supported

search the motif you want ( I use CTCF as an example)

Click save on the bottom right.

Click my motifs on the upper right.

export to bed or other format files.

several tools for Hi-C data, ChIP-seq and methylation data

See the links below:
http://omictools.com/differential-peak-calling/
http://omictools.com/chip-seq-and-beyond/3c-4c-5c-hi-c/hicup-s1473.html
http://omictools.com/analytical-pipelines8/
http://omictools.com/dmr/

There are too many tools out there for established high-throughput sequencing  based techniques.
1. you need to know what are out there. never try to re-invent the wheel
2. choose the right one to fit your own analysis. (only if none of them does what you want to do, write a script by yourself).  choose one with good documentation and is being actively developed.

manipulate GTF file with AWK

AWK is awesome! or I should say linux is awesome!
I had a post for linux one-liners here  http://crazyhottommy.blogspot.com/2013/10/useful-unix-one-liner-for-bioinformatics.html

Today, I  came across a blog series about GTF file manipulation using AWK. I feel it is cool, it teaches awk using a practical example. you will see the beauty of linux commands!
here are the links:
http://reasoniamhere.com/2013/09/16/awk-gtf-how-to-analyze-a-transcriptome-like-a-pro-part-1/
http://reasoniamhere.com/2013/09/17/awk-gtf-how-to-analyze-a-transcriptome-like-a-pro-part-2/
http://reasoniamhere.com/2013/09/18/awk-gtf-how-to-analyze-a-transcriptome-like-a-pro-part-3/

Enjoy!

MeDIP-seq and histone modification ChIP-seq analysis

I was reading a cell paper:

Dnmt3L Antagonizes DNA Methylation at Bivalent Promoters and Favors DNA Methylation at Gene Bodies in ESCs

http://www.sciencedirect.com.lp.hscl.ufl.edu/science/article/pii/S0092867413010878
I want to re-analyze the MeDIP-seq data and the histone modification ChIP-seq data.

I downloaded the MeDIP-seq data here: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE44644
and the H3k4me3 and H3k27me3 ChIP-seq data here: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE15519

The paper claimed that Dnmt3L maintains DNA hypomethylation at bivalent promoters ( with both H3k4me3 and H3k27me3 marks).  The authors used a Dnmt3L deficient mouse ES cell line and did MeDIP-seq on it. They found that many bivalent promoters gain DNA methylation in Dnmt3L negative cells. In contrast, DNA methylation loss mainly occurs in gene body.

I will work on the histone modification ChIP-seq data first.
I took a look at many different H3k4me3 and H3k27me3 ChIP-seq data sets in mES cells. H3k4me3 ChIP-seq data quality is usually very good with very specific and sharp peaks after MACS peak calling. However, H3k27me3 data are not that great. I did regular H3k4me3 ChIP, the antibody is very good and the enrichment of input is usually above 10% in the active loci I tested.

Similar to H3k36me3, an active mark in gene body, H3k27me3 (repressive) mark usually form broader but less-enriched regions. I played around with the MACS argument setting.

I will skip the process of transforming sra format data to fastq and mapping them to the reference genome. I have several posts for that already:
http://crazyhottommy.blogspot.com/2013/09/the-backtick-in-linux-command.html
bowtie for ChIp-seq http://crazyhottommy.blogspot.com/2013/05/sra-file-and-my-first-bowtie-run-on-uf.html
Tophat for RNA-seq http://crazyhottommy.blogspot.com/2013/10/rna-seq-analysis-samtools-sort-and.html

For H3k4me3 data, I only downloaded one, and called peaks:
[mtang@dev1 H3k4me3]$ macs -t SRR038984.sam -n initial_test -g mm 
INFO  @ Thu, 23 Jan 2014 13:37:47: 
# ARGUMENTS LIST:
# name = initial_test
# format = AUTO
# ChIP-seq file = SRR038984.sam
# control file = None
# effective genome size = 1.87e+09
# band width = 300
# model fold = 10,30
# pvalue cutoff = 1.00e-05
# Large dataset will be scaled towards smaller dataset.
# Range for calculating regional lambda is: 10000 bps

INFO  @ Thu, 23 Jan 2014 13:37:47: #1 read tag files... 
INFO  @ Thu, 23 Jan 2014 13:37:47: #1 read treatment tags... 
INFO  @ Thu, 23 Jan 2014 13:37:47: Detected format is: SAM 
INFO  @ Thu, 23 Jan 2014 13:37:56:  1000000 
INFO  @ Thu, 23 Jan 2014 13:38:05:  2000000 
INFO  @ Thu, 23 Jan 2014 13:38:13:  3000000 
INFO  @ Thu, 23 Jan 2014 13:38:21:  4000000 
INFO  @ Thu, 23 Jan 2014 13:38:29:  5000000 
INFO  @ Thu, 23 Jan 2014 13:38:38:  6000000 
INFO  @ Thu, 23 Jan 2014 13:38:46:  7000000 
WARNING @ Thu, 23 Jan 2014 13:38:56: NO records for chromosome chr16_random, minus strand! 
INFO  @ Thu, 23 Jan 2014 13:38:56: #1 tag size is determined as 36 bps 
INFO  @ Thu, 23 Jan 2014 13:38:56: #1 tag size = 36 
INFO  @ Thu, 23 Jan 2014 13:38:56: #1  total tags in treatment: 5000614 
INFO  @ Thu, 23 Jan 2014 13:38:56: #1 user defined the maximum tags... 
INFO  @ Thu, 23 Jan 2014 13:38:56: #1 filter out redundant tags at the same location and the same strand by allowing at most 1 tag(s) 
INFO  @ Thu, 23 Jan 2014 13:38:57: #1  tags after filtering in treatment: 3982997 
INFO  @ Thu, 23 Jan 2014 13:38:57: #1  Redundant rate of treatment: 0.20 
INFO  @ Thu, 23 Jan 2014 13:38:57: #1 finished! 
INFO  @ Thu, 23 Jan 2014 13:38:57: #2 Build Peak Model... 
INFO  @ Thu, 23 Jan 2014 13:39:15: #2 number of paired peaks: 22337 
INFO  @ Thu, 23 Jan 2014 13:39:16: #2 finished! 
INFO  @ Thu, 23 Jan 2014 13:39:16: #2 predicted fragment length is 164 bps 
INFO  @ Thu, 23 Jan 2014 13:39:16: #2.2 Generate R script for model : initial_test_model.r 
INFO  @ Thu, 23 Jan 2014 13:39:16: #3 Call peaks... 
INFO  @ Thu, 23 Jan 2014 13:39:16: #3 shift treatment data 
INFO  @ Thu, 23 Jan 2014 13:39:17: #3 merge +/- strand of treatment data 
INFO  @ Thu, 23 Jan 2014 13:39:18: #3 call peak candidates 
INFO  @ Thu, 23 Jan 2014 13:40:14: #3 use self to calculate local lambda and  filter peak candidates... 
INFO  @ Thu, 23 Jan 2014 13:40:27: #3 Finally, 32919 peaks are called! 
INFO  @ Thu, 23 Jan 2014 13:40:27: #4 Write output xls file... initial_test_peaks.xls 
INFO  @ Thu, 23 Jan 2014 13:40:27: #4 Write peak bed file... initial_test_peaks.bed 
INFO  @ Thu, 23 Jan 2014 13:40:27: #4 Write summits bed file... initial_test_summits.bed 
INFO  @ Thu, 23 Jan 2014 13:40:27: #5 Done! Check the output files! 

MACS successfully built the model and identified 32919 peaks.

For H3k27me3 data, there are two replicates, I merged them with samtools merge first

[mtang@dev1 H3k27me3]$ samtools merge merged_H3k27me3.bam SRR038975.bam SRR038976.bam 

then at HPC@UFL:

[mtang@dev1 H3k27me3]$ macs -t merged_H3k27me3.bam  -n initial_test -g mm 

it gave me only 4941 peaks.

[mtang@dev1 H3k27me3]$ wc -l merged_H3k27me3_peaks.bed 
4941 merged_H3k27me3_peaks.bed

I changed the arguments:
[mtang@dev1 H3k27me3]$ macs -t merge.bam -n no_model  -g mm --nomodel --shiftsize 73 --pvalue 1e-3 

note that I specify no peaking model building, shift size 73 (  fragment size 73x2=146 bp, a nuclesome size), and a less stringent pvalue 1e-3 ( the default is 1e-5)

This time I got 24522 peaks.

[mtang@dev1 H3k27me3]$ wc -l no_model_peaks.bed 

24522 no_model_peaks.bed

In my next post, I will use those two peak bed files generated by MACS to find the bivalent promoters.

=================================================================

I will switch gear and analyze the MeDIP-seq data. MeDIP-seq data are essentially similar to ChIP-seq data. Although there are specific software to analyze it, MACS  works well with this kind of data.

see a list of programs here http://seqanswers.com/wiki/Software/list

I have  MeDIP_shGFP.sam and MeDIP-shDnmt3L.sam files ready. 

To visualize the MeDIP data in IGV, I have the bedgraph files generated by MACS respectively. I also want to have a bedgraph file that with MeDIP-shGFP signal subtracted from MeDIP-shDnmt3L signal.

MACS2 has a bdgcpm program for that purpose

https://github.com/taoliu/MACS/wiki/Build-Signal-Track#wiki-Run_MACS2_bdgcmp_to_generate_foldenrichment_and_logLR_track

macs2 bdgcmp -t shDnmt3L_treat_afterfiting_all.bdg -c shGFP_treat_afterfiting_all.bdg -o subtracted.bdg -m subtract

Now, I have almost every file for visualization. Just to reproduce fig4E in the paper:

The first track is the MeDIP-seq signal for shGFP; the second is the MeDIP-seq signal for shDnmt3L; the third is the shDnmt3L signal subtracted with shGFP signal; the fourth track is the H3k4me3 ChIP-seq signal; the fifth is the H3k4me3 peaks called by MACS; the sixth is the H3k27me3 ChIP-seq signal; the last is the H3k27me3 peaks called by MACS.

The next question is to find the differentially changed peaks after ablating Dnmt3L. I have a post discussing this topic before http://crazyhottommy.blogspot.com/2013/10/compare-chip-seq-data-for-different.html

MACS2 has some new implements for this kind of purpose, but it requires a IgG or Input control file for each condition. https://github.com/taoliu/MACS/wiki/Call-differential-binding-events

I only have two bam files, so I will use http://code.google.com/p/diffreps/

diffReps needs  bed files to be the input. I converted the bam files to bed files with bedtools first:

http://code.google.com/p/bedtools/wiki/Usage#bamToBed

[mtang@dev1 MeDIP-seq]$ bamToBed -i SRR764931_shGFP.bam > SRR764931_shGFP.bed

[mtang@dev1 MeDIP-seq]$ bamToBed -i SRR764932_shDnmt3L.bam > SRR764932_shDnmt3L.bed

since I do not have biological replicates, I used G-test 

[mtang@dev1 MeDIP-seq]$ diffReps.pl -tr SRR764932_shDnmt3L.bed -co SRR764931_shGFP.bed -gn mm9 -re diff.nb.txt -me gt

It generates three files:

diff.nb.txt

diff.nb.txt.annotated

diff.nb.txt.hotspot

for more information of these three files, see the diffReps link above.

ChIP-seq analysis programs

This is the first blog in the year of 2014.  I am graduating in the coming August, there are so many things to do and I do not have much time to blog. However, I do not want it to be inactive.

I have extensive experience with ChIP-seq data analysis, and I want to list all the programs that I usually use.

With raw Fastq data, one needs to align it to the reference genome by an aligner. I usually use Bowtie http://bowtie-bio.sourceforge.net/index.shtml, others like BWA (developed by Heng Li, a legendary bioinformatician), MAQ, SOAP etc are also very popular (http://massgenomics.org/short-read-aligners)

Many times, I download data from SRA, sra toolkit http://eutils.ncbi.nih.gov/Traces/sra/?view=software is used to convert the sra format to fastq, and then one can align it with bowtie.

After mapping with Bowtie, one get a sam file http://genome.sph.umich.edu/wiki/SAM. samtools (also developed by Heng Li) is used to convert it to bam file, a binary form of sam file.

MACS https://github.com/taoliu/MACS/ is the most widely used peak calling program developed by Tao Li, previously in shirley Liu's lab at Harvard. MACS2 now can detect differentially changed peaks (see the link above). It also generates a bedgraph file to be visualized by IGV http://www.broadinstitute.org/igv/.

For motif enrichment analysis:
Homer  http://homer.salk.edu/homer/ngs/index.html
pscan-ChIP http://159.149.160.51/pscan_chip_dev/
MEME http://meme.nbcr.net/meme/
GREAT for GO analysis http://bejerano.stanford.edu/great/public/html/
Homer, CEAS http://liulab.dfci.harvard.edu/CEAS/ or PAVIS http://manticore.niehs.nih.gov:8080/pavis/ to annotate peaks

you might be interested in my previous blogs also:
http://crazyhottommy.blogspot.com/2013/12/chip-seq-peaks-overlapping-significance.html
http://crazyhottommy.blogspot.com/2013/10/compare-chip-seq-data-for-different.html
http://crazyhottommy.blogspot.com/2013/08/how-to-make-heatmap-based-on-chip-seq.html
http://crazyhottommy.blogspot.com/2013/04/how-to-make-tss-plot-using-rna-seq-and.html

2014 will be an exciting year for me! Follow me and many others I follow on twitter https://twitter.com/tangming2005.