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.
This blog by Tommy Tang is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Monday, February 24, 2014
Tuesday, February 18, 2014
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
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...
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:
- 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.
- Remove any existing 'track' or 'browser' lines from your wig file so that it contains only data.
- Download the wigToBigWig program from the directory of binary utilities.
- 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.
- 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.) - Move the newly created bigWig file (myBigWig.bw) to a http, https, or ftp location.
- 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.bwautoScale <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. - Paste this custom track line into the text box in the custom track management page.
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. |
Thursday, February 13, 2014
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.
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.
Thursday, February 6, 2014
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.
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.
Tuesday, January 28, 2014
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!
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!
Thursday, January 23, 2014
MeDIP-seq and histone modification ChIP-seq analysis
I was reading a cell paper:
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!
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
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/S0092867413010878I 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
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:
[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.
Monday, January 20, 2014
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.
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.
Thursday, December 19, 2013
log2 fold gene expression change heatmap
I was looking at two different microarray data sets. In each experiment, the authors knocked down a transcription factor (TF) and examined the gene expression change respectively. I wanted to know whether there are any co-upregulated or co-downregulated genes after knocking down the TFs.
I analyzed the microarray data by bioconductor. See one of the examples I blogged before
http://crazyhottommy.blogspot.com/2013/12/geoquery-to-access-geo-datasets.html
Because those are two different experiments done in different platforms and in different labs, I can not put the raw expression data together. Instead, I got the log2 fold change (KD vs Control) for each data set.
I got a topgene list for each data set from the limma topTable. The third column is the gene name, the fourth column is the logFC. I cut these two columns out, delete the first line (the header) and get rid of the double quote around the gene names:
tommy@tommy-ThinkPad-T420:~/tommy$ cut -f3,4 topgenes_shTF1_FDR0.05_lfc1.txt | sed '1 d' | sed 's/"//g' | sort | uniq > shTF1_topgenes.txt
Inner join these two lists together produced a file with three columns. After join command, the file became white space delimited, I changed it back to tab delimited.
tommy@tommy-ThinkPad-T420:~/tommy$ join shTF1_topgenes.txt shTF2_topgenes.txt | tr ' ' '\t' > co_regulated.txt
I analyzed the microarray data by bioconductor. See one of the examples I blogged before
http://crazyhottommy.blogspot.com/2013/12/geoquery-to-access-geo-datasets.html
Because those are two different experiments done in different platforms and in different labs, I can not put the raw expression data together. Instead, I got the log2 fold change (KD vs Control) for each data set.
I got a topgene list for each data set from the limma topTable. The third column is the gene name, the fourth column is the logFC. I cut these two columns out, delete the first line (the header) and get rid of the double quote around the gene names:
tommy@tommy-ThinkPad-T420:~/tommy$ cut -f3,4 topgenes_shTF1_FDR0.05_lfc1.txt | sed '1 d' | sed 's/"//g' | sort | uniq > shTF1_topgenes.txt
tommy@tommy-ThinkPad-T420:~/tommy$ cut -f3,4 topgenes_shTF2.txt | sed '1 d' | sed 's/"//g' | sort | uniq > shTF2_topgenes.txt
Inner join these two lists together produced a file with three columns. After join command, the file became white space delimited, I changed it back to tab delimited.
tommy@tommy-ThinkPad-T420:~/tommy$ join shTF1_topgenes.txt shTF2_topgenes.txt | tr ' ' '\t' > co_regulated.txt
only keep the co-upregualted and co-downregulated genes:
tommy@tommy-ThinkPad-T420:~/tommy$ cat co_regulated.txt | awk '{if($2>0 && $3>0 || $2<0 && $3<0) print}' | head
Abhd5 1.22313483440383 0.718522081750001
Acta2 2.42294712761151 0.602606251249999
Acta2 2.7279547932607 0.602606251249999
Acta2 2.92237008906738 0.602606251249999
Akt3 1.62314377335472 0.717014153000002
Akt3 1.62314377335472 0.900214081749999
Akt3 1.62314377335472 1.10733403775
Akt3 2.21617654271768 0.717014153000002
Akt3 2.21617654271768 0.900214081749999
Akt3 2.21617654271768 1.10733403775
There are multiple probes for the same gene, only keep the first instance for each gene:
tommy@tommy-ThinkPad-T420:~/tommy$ cat co_up_or_down.txt | awk '!a[$1]++' > co_up_or_down_uniq.txt
Among these 116 genes, 85 are co-upregulated, 31 are co-down regulated.
Now, I can draw a heatmap based on the log2 fold changes.
The code above are heavily commented, you should be able to follow easily.
I once asked about the colour setting in the heatmap in Seqanwser http://seqanswers.com/forums/showthread.php?t=33023
Also see here http://stackoverflow.com/questions/17820143/how-to-change-heatmap-2-color-range-in-r
In the script, I mapped the colours to values:
So, the black and green coloured cells are down-regulated, the red cells are up-regulated.
After add the side bar:
In the dendrogram, I split the genes to two groups (co-up and co-down). In this particular example, I already pre-selected co-up and co-down regulated genes before I constructed the heatmap. If one did not do this and wants to find the clustered (heatmap.2 use hierarchical clustering as default, type ?hclust in R to see more) groups of gene members, one can find them by the code above.
The easiest way is to use cluster3.0 to cluster and Java tree viewer to visualize the result http://bonsai.hgc.jp/~mdehoon/software/cluster/software.htm
It is interactive, you can choose the sub-groups just by selecting them directly using mouse.
Heatmap is no mystery. It just uses colours to represent the values for visualization. However, you need to know what you are doing. Make sure you know the different distance calculations among rows http://stat.ethz.ch/R-manual/R-patched/library/stats/html/dist.html
and the different clustering methods (hierarchical, K-means etc)
http://stat.ethz.ch/R-manual/R-patched/library/stats/html/hclust.html
http://stat.ethz.ch/R-manual/R-devel/library/stats/html/kmeans.html
See also my previous post http://crazyhottommy.blogspot.com/2013/08/how-to-make-heatmap-based-on-chip-seq.html
I once asked about the colour setting in the heatmap in Seqanwser http://seqanswers.com/forums/showthread.php?t=33023
Also see here http://stackoverflow.com/questions/17820143/how-to-change-heatmap-2-color-range-in-r
In the script, I mapped the colours to values:
bk = unique(c(seq(-5.73,-0.5, length=100),seq(-0.5,0, length=100), seq(0,5,length=100)))hmcols<- colorRampPalette(c("green","black", "red"))(length(bk)-1)
So, the black and green coloured cells are down-regulated, the red cells are up-regulated.
The easiest way is to use cluster3.0 to cluster and Java tree viewer to visualize the result http://bonsai.hgc.jp/~mdehoon/software/cluster/software.htm
It is interactive, you can choose the sub-groups just by selecting them directly using mouse.
Heatmap is no mystery. It just uses colours to represent the values for visualization. However, you need to know what you are doing. Make sure you know the different distance calculations among rows http://stat.ethz.ch/R-manual/R-patched/library/stats/html/dist.html
and the different clustering methods (hierarchical, K-means etc)
http://stat.ethz.ch/R-manual/R-patched/library/stats/html/hclust.html
http://stat.ethz.ch/R-manual/R-devel/library/stats/html/kmeans.html
See also my previous post http://crazyhottommy.blogspot.com/2013/08/how-to-make-heatmap-based-on-chip-seq.html
Tuesday, December 17, 2013
ChIP-seq peaks overlapping significance test
The question is straightforward:
I have two ChIP-seq bed files (transcription factor binding) obtained by MACS peak calling, and I want to know whether these two transcription factors co-bind to the same sites with a probability greater than expected.
The answer is not trivial.
See a post here http://www.biostars.org/p/5484/
update on 10/30/2014: one can use an R package for this genometricorr
One can use bedtools shuffle http://bedtools.readthedocs.org/en/latest/content/tools/shuffle.html developed by Aaron Quinlan@UVA as suggested by the link above, but it is not automatic.
The coming bedtools2 seems to have this function implemented http://www.ashg.org/2013meeting/abstracts/fulltext/f130122205.htm
"In the context of biological discovery, the most exciting new functionality in bedtools2 is a comprehensive set of statistical measures for revealing associations between sets of genomic features (e.g., do these transcription factors co-associate more than expected by chance?). Here we present the new statistical tests in bedtools2, compare our tests to existing approaches such as the ENCODE Genome Structure Correction (GSC) metric, and provide needed insights into which metrics are most appropriate to common biological questions. We demonstrate typical misuses of these metrics and illustrate how our tests and associated visualization tools can reveal new biological insights. Given the speed and analytical flexibility of bedtools2, we anticipate that our new toolkit will be an invaluable resource for geneticists studying the impact of genetic variation and regulatory elements on human disease phenotypes."
I am looking forward to it.
Luckily, I found other two tools that automate this process and are very handy to use.
poverlap https://github.com/brentp/poverlap
and
Genomic Association Tester (GAT) http://www.cgat.org/~andreas/documentation/gat/contents.html
I will use GAT as a demonstration since it is very well documented and seems to be more powerful.
Following the tutorial here http://www.cgat.org/~andreas/documentation/gat/tutorialIntervalOverlap.html
I will need three bed files ready: a.bed file for TF1, b.bed file for TF2 and a workspace (genome) bed file that specifies the chromosome coordinates limits.
I have peak bed files for TFs, and I need to get the workspace bedfile.
See a post here http://www.biostars.org/p/16396/#16399
I downloaded the chromInfo.txt.gz file here ( I am working with mouse data)
http://hgdownload.soe.ucsc.edu/goldenPath/mm9/database/
Some manipulation to get the bed file:
tommy@tommy-ThinkPad-T420:~/mm9_genome_info$ cat chromInfo.txt
chr1 197195432 /gbdb/mm9/mm9.2bit
chr2 181748087 /gbdb/mm9/mm9.2bit
chr3 159599783 /gbdb/mm9/mm9.2bit
chr4 155630120 /gbdb/mm9/mm9.2bit
chr5 152537259 /gbdb/mm9/mm9.2bit
chr6 149517037 /gbdb/mm9/mm9.2bit
chr7 152524553 /gbdb/mm9/mm9.2bit
chr8 131738871 /gbdb/mm9/mm9.2bit
chr9 124076172 /gbdb/mm9/mm9.2bit
chr10 129993255 /gbdb/mm9/mm9.2bit
chr11 121843856 /gbdb/mm9/mm9.2bit
chr12 121257530 /gbdb/mm9/mm9.2bit
chr13 120284312 /gbdb/mm9/mm9.2bit
chr14 125194864 /gbdb/mm9/mm9.2bit
chr15 103494974 /gbdb/mm9/mm9.2bit
chr16 98319150 /gbdb/mm9/mm9.2bit
chr17 95272651 /gbdb/mm9/mm9.2bit
chr18 90772031 /gbdb/mm9/mm9.2bit
chr19 61342430 /gbdb/mm9/mm9.2bit
chrX 166650296 /gbdb/mm9/mm9.2bit
chrY 15902555 /gbdb/mm9/mm9.2bit
chrM 16299 /gbdb/mm9/mm9.2bit
chr13_random 400311 /gbdb/mm9/mm9.2bit
chr16_random 3994 /gbdb/mm9/mm9.2bit
chr17_random 628739 /gbdb/mm9/mm9.2bit
chr1_random 1231697 /gbdb/mm9/mm9.2bit
chr3_random 41899 /gbdb/mm9/mm9.2bit
chr4_random 160594 /gbdb/mm9/mm9.2bit
chr5_random 357350 /gbdb/mm9/mm9.2bit
chr7_random 362490 /gbdb/mm9/mm9.2bit
chr8_random 849593 /gbdb/mm9/mm9.2bit
chr9_random 449403 /gbdb/mm9/mm9.2bit
chrUn_random 5900358 /gbdb/mm9/mm9.2bit
chrX_random 1785075 /gbdb/mm9/mm9.2bit
chrY_random 58682461 /gbdb/mm9/mm9.2bit
Notice the weird chromosome names like chrUn and chr4_random. See a post here http://www.biostars.org/p/9577/
On HPC@UFL:
[mtang@dev1 mm9_genome_info]$ gat-run.py --segments=a.bed --annotations=b.bed --workspace=mm9_contigs.bed --ignore-segment-tracks --num-samples=1000 --log=gat.log > gat.tsv
a.bed file contains about 30,000 peaks and b.bed file contains 80,000 peaks. It took me 1 min to finish 1000 simulations:
For 1000 simulation, I got a minimal p value 0.001 indicating that these two TFs co-occupy the same genomic sites not randomly.
I have two ChIP-seq bed files (transcription factor binding) obtained by MACS peak calling, and I want to know whether these two transcription factors co-bind to the same sites with a probability greater than expected.
The answer is not trivial.
See a post here http://www.biostars.org/p/5484/
update on 10/30/2014: one can use an R package for this genometricorr
One can use bedtools shuffle http://bedtools.readthedocs.org/en/latest/content/tools/shuffle.html developed by Aaron Quinlan@UVA as suggested by the link above, but it is not automatic.
The coming bedtools2 seems to have this function implemented http://www.ashg.org/2013meeting/abstracts/fulltext/f130122205.htm
"In the context of biological discovery, the most exciting new functionality in bedtools2 is a comprehensive set of statistical measures for revealing associations between sets of genomic features (e.g., do these transcription factors co-associate more than expected by chance?). Here we present the new statistical tests in bedtools2, compare our tests to existing approaches such as the ENCODE Genome Structure Correction (GSC) metric, and provide needed insights into which metrics are most appropriate to common biological questions. We demonstrate typical misuses of these metrics and illustrate how our tests and associated visualization tools can reveal new biological insights. Given the speed and analytical flexibility of bedtools2, we anticipate that our new toolkit will be an invaluable resource for geneticists studying the impact of genetic variation and regulatory elements on human disease phenotypes."
I am looking forward to it.
Luckily, I found other two tools that automate this process and are very handy to use.
poverlap https://github.com/brentp/poverlap
and
Genomic Association Tester (GAT) http://www.cgat.org/~andreas/documentation/gat/contents.html
I will use GAT as a demonstration since it is very well documented and seems to be more powerful.
Following the tutorial here http://www.cgat.org/~andreas/documentation/gat/tutorialIntervalOverlap.html
I will need three bed files ready: a.bed file for TF1, b.bed file for TF2 and a workspace (genome) bed file that specifies the chromosome coordinates limits.
I have peak bed files for TFs, and I need to get the workspace bedfile.
See a post here http://www.biostars.org/p/16396/#16399
I downloaded the chromInfo.txt.gz file here ( I am working with mouse data)
http://hgdownload.soe.ucsc.edu/goldenPath/mm9/database/
Some manipulation to get the bed file:
tommy@tommy-ThinkPad-T420:~/mm9_genome_info$ cat chromInfo.txt
chr1 197195432 /gbdb/mm9/mm9.2bit
chr2 181748087 /gbdb/mm9/mm9.2bit
chr3 159599783 /gbdb/mm9/mm9.2bit
chr4 155630120 /gbdb/mm9/mm9.2bit
chr5 152537259 /gbdb/mm9/mm9.2bit
chr6 149517037 /gbdb/mm9/mm9.2bit
chr7 152524553 /gbdb/mm9/mm9.2bit
chr8 131738871 /gbdb/mm9/mm9.2bit
chr9 124076172 /gbdb/mm9/mm9.2bit
chr10 129993255 /gbdb/mm9/mm9.2bit
chr11 121843856 /gbdb/mm9/mm9.2bit
chr12 121257530 /gbdb/mm9/mm9.2bit
chr13 120284312 /gbdb/mm9/mm9.2bit
chr14 125194864 /gbdb/mm9/mm9.2bit
chr15 103494974 /gbdb/mm9/mm9.2bit
chr16 98319150 /gbdb/mm9/mm9.2bit
chr17 95272651 /gbdb/mm9/mm9.2bit
chr18 90772031 /gbdb/mm9/mm9.2bit
chr19 61342430 /gbdb/mm9/mm9.2bit
chrX 166650296 /gbdb/mm9/mm9.2bit
chrY 15902555 /gbdb/mm9/mm9.2bit
chrM 16299 /gbdb/mm9/mm9.2bit
chr13_random 400311 /gbdb/mm9/mm9.2bit
chr16_random 3994 /gbdb/mm9/mm9.2bit
chr17_random 628739 /gbdb/mm9/mm9.2bit
chr1_random 1231697 /gbdb/mm9/mm9.2bit
chr3_random 41899 /gbdb/mm9/mm9.2bit
chr4_random 160594 /gbdb/mm9/mm9.2bit
chr5_random 357350 /gbdb/mm9/mm9.2bit
chr7_random 362490 /gbdb/mm9/mm9.2bit
chr8_random 849593 /gbdb/mm9/mm9.2bit
chr9_random 449403 /gbdb/mm9/mm9.2bit
chrUn_random 5900358 /gbdb/mm9/mm9.2bit
chrX_random 1785075 /gbdb/mm9/mm9.2bit
chrY_random 58682461 /gbdb/mm9/mm9.2bit
tommy@tommy-ThinkPad-T420:~/mm9_genome_info$ cat chromInfo.txt | cut -f1-2 | awk '{print $1"\t"0"\t"$2}'
chr1 0 197195432
chr2 0 181748087
chr3 0 159599783
chr4 0 155630120
chr5 0 152537259
chr6 0 149517037
chr7 0 152524553
chr8 0 131738871
chr9 0 124076172
chr10 0 129993255
chr11 0 121843856
chr12 0 121257530
chr13 0 120284312
chr14 0 125194864
chr15 0 103494974
chr16 0 98319150
chr17 0 95272651
chr18 0 90772031
chr19 0 61342430
chrX 0 166650296
chrY 0 15902555
chrM 0 16299
chr13_random 0 400311
chr16_random 0 3994
chr17_random 0 628739
chr1_random 0 1231697
chr3_random 0 41899
chr4_random 0 160594
chr5_random 0 357350
chr7_random 0 362490
chr8_random 0 849593
chr9_random 0 449403
chrUn_random 0 5900358
chrX_random 0 1785075
chrY_random 0 58682461
tommy@tommy-ThinkPad-T420:~/Tet1/mm9_genome_info$ cat chromInfo.txt | cut -f1-2 | awk '{print $1"\t"0"\t"$2}' > mm9_contigs.bed
- The chr*_random sequences are unplaced sequence on those reference chromosomes.
- The chrUn_* sequences are unlocalized sequences where the corresponding reference chromosome has not been determined
On HPC@UFL:
[mtang@dev1 mm9_genome_info]$ gat-run.py --segments=a.bed --annotations=b.bed --workspace=mm9_contigs.bed --ignore-segment-tracks --num-samples=1000 --log=gat.log > gat.tsv
a.bed file contains about 30,000 peaks and b.bed file contains 80,000 peaks. It took me 1 min to finish 1000 simulations:
| observed | expected | CI95low | CI95high | stddev | fold | l2fold | pvalue | qvalue |
| 3998559 | 282473.979 | 267444 | 298581 | 9379.0988 | 14.1554 | 3.8233 | 1.0000e-03 | 1.0000e-03 |
For 1000 simulation, I got a minimal p value 0.001 indicating that these two TFs co-occupy the same genomic sites not randomly.
Wednesday, December 11, 2013
Tet1 enrichment at CpG rich-promoters
Continue with my last post http://crazyhottommy.blogspot.com/2013/12/geoquery-to-access-geo-datasets.html, I am going to analyze the Tet1 ChIP-seq data in the paper Dual functions of Tet1 in transcriptional regulation in mouse embryonic stem cells.
To reproduce fig1, I have to group the promoters into CpG-rich and CpG-poor promoters first.
Get the CpG island bed file and the refseq gene table first.
I got the file from UCSC table browser http://genome.ucsc.edu/cgi-bin/hgTables?command=start
choose Expression and regulation and CpG islands track, and download it.
similarly, get the refSeq gene table, send the output to Galaxy ( in order to get the -5kb promoter
regions, one needs to do some manipulation of the file, several options are there: bedtools, awk and Galaxy. The easiest way is through Galaxy if you are not familiar with linux command line)
Click get franks, and get the 5kb upstream coordinates, and cut columns that we need:
tommy@tommy-ThinkPad-T420:~/Tet1$ head 5kb_promoter_mm9.bed
chr1 134207701 134212701 NM_028778 +
chr1 134207701 134212701 NM_001195025 +
chr1 58752833 58757833 NM_175370 -
chr1 9289811 9294811 NM_027671 -
chr1 25886552 25891552 NM_175642 -
chr1 33726639 33731639 NM_008922 -
chr1 125942136 125947136 NM_199021 -
chr1 75503027 75508027 NM_178884 -
chr1 109057691 109062691 NM_198680 -
chr1 184490137 184495137 NM_130890 +
you can do the same thing with CpG islands with Galaxy, but I show you the linux command below:
tommy@tommy-ThinkPad-T420:~/Tet1$ head mm9_CpG_island.txt
#bin chrom chromStart chromEnd name length cpgNum gcNum perCpg perGc obsExp
107 chr1 36568608 36569851 CpG: 128 1243 128 766 20.6 61.6 1.09
151 chr1 82313020 82313491 CpG: 56 471 56 319 23.8 67.7 1.06
234 chr1 169213895 169215283 CpG: 145 1388 145 958 20.9 69 0.88
611 chr1 3521705 3521924 CpG: 27 219 27 167 24.7 76.3 0.86
612 chr1 3660700 3661155 CpG: 34 455 34 293 14.9 64.4 0.72
612 chr1 3661735 3662237 CpG: 45 502 45 348 17.9 69.3 0.75
619 chr1 4481782 4483754 CpG: 165 1972 165 1286 16.7 65.2 0.79
619 chr1 4487028 4487689 CpG: 47 661 47 393 14.2 59.5 0.81
619 chr1 4561722 4562156 CpG: 44 434 44 265 20.3 61.1 1.09
I only need the second, third and fourth column of the CpG file:
tommy@tommy-ThinkPad-T420:~/Tet1$ cat mm9_CpG_island.txt | cut -f2-4 > mm9_CpG.bed
tommy@tommy-ThinkPad-T420:~/Tet1$ head mm9_CpG.bed
chr1 36568608 36569851
chr1 82313020 82313491
chr1 169213895 169215283
chr1 3521705 3521924
chr1 3660700 3661155
chr1 3661735 3662237
chr1 4481782 4483754
chr1 4487028 4487689
chr1 4561722 4562156
chr1 4679265 4679478
To reproduce fig1a, I have all the files ready: a sorted Tet1 ChIP-seq bam file, CpG_rich_promoters.bed and CpG_poor_promoters.bed
To assign functional binding sites to their target genes, one may consider to use BETA developed in Shirely Liu's lab http://cistrome.org/BETA/. It integrates ChIP-seq data and gene expression data to infer the TF target genes.
You may be also interested in this paper:
To reproduce fig1, I have to group the promoters into CpG-rich and CpG-poor promoters first.
Get the CpG island bed file and the refseq gene table first.
I got the file from UCSC table browser http://genome.ucsc.edu/cgi-bin/hgTables?command=start
choose Expression and regulation and CpG islands track, and download it.
similarly, get the refSeq gene table, send the output to Galaxy ( in order to get the -5kb promoter
regions, one needs to do some manipulation of the file, several options are there: bedtools, awk and Galaxy. The easiest way is through Galaxy if you are not familiar with linux command line)
tommy@tommy-ThinkPad-T420:~/Tet1$ head 5kb_promoter_mm9.bed
chr1 134207701 134212701 NM_028778 +
chr1 134207701 134212701 NM_001195025 +
chr1 58752833 58757833 NM_175370 -
chr1 9289811 9294811 NM_027671 -
chr1 25886552 25891552 NM_175642 -
chr1 33726639 33731639 NM_008922 -
chr1 125942136 125947136 NM_199021 -
chr1 75503027 75508027 NM_178884 -
chr1 109057691 109062691 NM_198680 -
chr1 184490137 184495137 NM_130890 +
tommy@tommy-ThinkPad-T420:~/Tet1$ head mm9_CpG_island.txt
#bin chrom chromStart chromEnd name length cpgNum gcNum perCpg perGc obsExp
107 chr1 36568608 36569851 CpG: 128 1243 128 766 20.6 61.6 1.09
151 chr1 82313020 82313491 CpG: 56 471 56 319 23.8 67.7 1.06
234 chr1 169213895 169215283 CpG: 145 1388 145 958 20.9 69 0.88
611 chr1 3521705 3521924 CpG: 27 219 27 167 24.7 76.3 0.86
612 chr1 3660700 3661155 CpG: 34 455 34 293 14.9 64.4 0.72
612 chr1 3661735 3662237 CpG: 45 502 45 348 17.9 69.3 0.75
619 chr1 4481782 4483754 CpG: 165 1972 165 1286 16.7 65.2 0.79
619 chr1 4487028 4487689 CpG: 47 661 47 393 14.2 59.5 0.81
619 chr1 4561722 4562156 CpG: 44 434 44 265 20.3 61.1 1.09
tommy@tommy-ThinkPad-T420:~/Tet1$ cat mm9_CpG_island.txt | cut -f2-4 > mm9_CpG.bed
tommy@tommy-ThinkPad-T420:~/Tet1$ head mm9_CpG.bed
chr1 36568608 36569851
chr1 82313020 82313491
chr1 169213895 169215283
chr1 3521705 3521924
chr1 3660700 3661155
chr1 3661735 3662237
chr1 4481782 4483754
chr1 4487028 4487689
chr1 4561722 4562156
chr1 4679265 4679478
How many CpG island and refGenes are there?
tommy@tommy-ThinkPad-T420:~/Tet1$ wc -l 5kb_promoter_mm9.bed mm9_CpG.bed
32555 5kb_promoter_mm9.bed
16026 mm9_CpG.bed
48581 total
Now use bedtools (you can also use the intersect, substract tool in Galaxy) to get promoters that contain CpG island:
tommy@tommy-ThinkPad-T420:~/Tet1$ bedtools intersect -a 5kb_promoter_mm9.bed -b mm9_CpG.bed -wa -u | wc -l
16171
tommy@tommy-ThinkPad-T420:~/Tet1$ bedtools intersect -a 5kb_promoter_mm9.bed -b mm9_CpG.bed -wa -u > CpG_rich_promoters.bed
promoters without CpG island:
tommy@tommy-ThinkPad-T420:~/Tet1$ bedtools intersect -a 5kb_promoter_mm9.bed -b mm9_CpG.bed -v | wc -l
16384
tommy@tommy-ThinkPad-T420:~/Tet1$ bedtools intersect -a 5kb_promoter_mm9.bed -b mm9_CpG.bed -v > CpG_poor_promoters.bed
I have mapped the Tet1 ChIP-seq raw sequence data by bowtie at HPC@UFL. I got the sam file and converted it to sorted bam file with samtools. Then, I called peaks with MACS1.4, and got the peaks and bedgraph file that I can visualize in IGV genome browser (see my previous post http://crazyhottommy.blogspot.com/2013/11/chip-exo-data-analysis.html).
A partial reproduction of fig1b in the paper:
the first track is the IgG control, the second is the Tet1 bedgraph signal, the third is the Tet1 peaks called by MACS, and the last one is the mm9_CpG.bed track. I observed many CpG islands are overlapped with the Tet1 peaks. The peaks are very sharp indicating a successful ChIP-seq experiment.
There are total 30946 peaks that are called by MACS. The number is 35,564 in the paper, but they used an older version of MACS. Nevertheless, the numbers are very close.
tommy@tommy-ThinkPad-T420:~/Tet1$ wc -l Tet1_ChIP-seq_peaks.bed
30946 Tet1_ChIP-seq_peaks.bed
using CEAS http://liulab.dfci.harvard.edu/CEAS/usermanual.html from Shirely Liu's lab, I investigated the distribution of the peaks:
tommy@tommy-ThinkPad-T420:~/Desktop/Gene_Annotation_for_CEAS_package$ ceas -g mm9.refGene -b ~/Tet1/Tet1_ChIP-seq_peaks.bed
Most peaks are located at gene-rich regions, only 23.3% of them are located at distal intergenic regions. It will be very interesting to study the function of these distal Tet1 binding sites.
A significant portion of peaks (33.2%) are at introns, indicating that Tet1 plays a role in DNA methylation in the gene bodies.
Note that CEAS can also produce figures as in fig1a, but it requires a wig file as input. I will stick to my python script which only need a bam file ( to convert bam/sam files to wig file, see here
To reproduce fig1a, I have all the files ready: a sorted Tet1 ChIP-seq bam file, CpG_rich_promoters.bed and CpG_poor_promoters.bed
See the python code below, it uses HTSeq python package written by Simon Anders. http://www.embl.de/research/units/genome_biology/huber/members/?s_personId=CP-60002977. He is an expert in RNA-Seq analysis ( He wrote the famous DESeq package)
it generates a figure as shown below:
One can count the mapped the reads number by samtools flagstat first, and then plug in the number in the above script. It is pretty slow to count the mapped reads by the aligned_counts function above.
I tricked by counting the reads first by samtools and the whole script only took me 3 minutes to run which is pretty fast!
Different from fig1a in the original paper, I used normalized counts per million for Y axis. I do not know why the authors used the peak P value as Y axis. Anyway, the conclusion is still the same: Tet1 preferentially binds to CpG-rich promoters. A close inspection of the figure generated by ipython ( interactive visualization) reveals that Tet1 peaks at around 100bp downstream of TSS.
I also reproduced part of Fig3a where the authors show Tet1 binds to bivalent promoters:
Annotate the peaks by Homer:
tommy@tommy-ThinkPad-T420:~/Tet1$ annotatePeaks.pl Tet1_ChIP-seq_peaks.bed mm9 -noann > Tet1_targeted_genes.txt
Homer finds the closest TSS for each peak. If I set the closest distance limit to 10kb upstream and downstream of TSS, there are still around ~12,000 genes that are associated with these peaks (total 30946 peaks).
the 10th column is the distance to TSS, the 16th column is the gene name
tommy@tommy-ThinkPad-T420:~/Tet1$ cat Tet1_targeted_genes.txt | awk '{ if ($10<10000 && $10>-10000) print}'| cut -f16 | sort | uniq | wc -l
12202
it generates a figure as shown below:
I tricked by counting the reads first by samtools and the whole script only took me 3 minutes to run which is pretty fast!
Different from fig1a in the original paper, I used normalized counts per million for Y axis. I do not know why the authors used the peak P value as Y axis. Anyway, the conclusion is still the same: Tet1 preferentially binds to CpG-rich promoters. A close inspection of the figure generated by ipython ( interactive visualization) reveals that Tet1 peaks at around 100bp downstream of TSS.
I also reproduced part of Fig3a where the authors show Tet1 binds to bivalent promoters:
Annotate the peaks by Homer:
tommy@tommy-ThinkPad-T420:~/Tet1$ annotatePeaks.pl Tet1_ChIP-seq_peaks.bed mm9 -noann > Tet1_targeted_genes.txt
Homer finds the closest TSS for each peak. If I set the closest distance limit to 10kb upstream and downstream of TSS, there are still around ~12,000 genes that are associated with these peaks (total 30946 peaks).
the 10th column is the distance to TSS, the 16th column is the gene name
tommy@tommy-ThinkPad-T420:~/Tet1$ cat Tet1_targeted_genes.txt | awk '{ if ($10<10000 && $10>-10000) print}'| cut -f16 | sort | uniq | wc -l
12202
Summary:
1. It is a very successful ChIP-seq experiment with sharp peaks over IgG background. Majority of Tet1 binding sites are at gene-rich regions. 23.3% of Tet1 binding sites are in distal intergenic regions with unknow functions. Tet1 tends to bind to CpG-rich promoters.
2. If one assigns the nearest genes to the peaks, ~12,000 genes are associated with Tet1 binding sites. compare with the differentially expressed genes in my last post-only 1292 genes (~10%) are affected by Tet1 KD- the number of genes that are regulated by Tet1 is much smaller than the peak associated ones.
Several reasons: KD is not complete (~20% residue Tet1 remains according to the microarray data), the remaining Tet1 can still maintain expression of some genes. If we knock out Tet1 (by TALENT or CRISPR http://zlab.mit.edu/publications.html), it may have a bigger effect on gene expression. In addition, RNA-seq may reveal more differentially expressed genes compared with microarray.
Another reason is that binding does not infer function (especially for weak binding sites) and I think this is probably the mostly likely reason. Several papers have shown that changes of adjacent TF binding poorly correlates with gene expression change:
1. It is a very successful ChIP-seq experiment with sharp peaks over IgG background. Majority of Tet1 binding sites are at gene-rich regions. 23.3% of Tet1 binding sites are in distal intergenic regions with unknow functions. Tet1 tends to bind to CpG-rich promoters.
2. If one assigns the nearest genes to the peaks, ~12,000 genes are associated with Tet1 binding sites. compare with the differentially expressed genes in my last post-only 1292 genes (~10%) are affected by Tet1 KD- the number of genes that are regulated by Tet1 is much smaller than the peak associated ones.
Several reasons: KD is not complete (~20% residue Tet1 remains according to the microarray data), the remaining Tet1 can still maintain expression of some genes. If we knock out Tet1 (by TALENT or CRISPR http://zlab.mit.edu/publications.html), it may have a bigger effect on gene expression. In addition, RNA-seq may reveal more differentially expressed genes compared with microarray.
Another reason is that binding does not infer function (especially for weak binding sites) and I think this is probably the mostly likely reason. Several papers have shown that changes of adjacent TF binding poorly correlates with gene expression change:
Extensive Divergence of Transcription Factor Binding in Drosophila Embryos with Highly Conserved Gene Expression
http://www.plosgenetics.org/article/info%3Adoi%2F10.1371%2Fjournal.pgen.1003748
Transcription Factors Bind Thousands of Active and Inactive Regions in theDrosophila Blastoderm
The Functional Consequences of Variation in Transcription Factor Binding
http://arxiv.org/abs/1310.5166
" On average, 14.7% of genes bound by a factor were differentially expressed following the knockdown of that factor, suggesting that most interactions between TF and chromatin do not result in measurable changes in gene expression levels of putative target genes. "To assign functional binding sites to their target genes, one may consider to use BETA developed in Shirely Liu's lab http://cistrome.org/BETA/. It integrates ChIP-seq data and gene expression data to infer the TF target genes.
You may be also interested in this paper:
Assessing Computational Methods for Transcription Factor Target Gene Identification Based on ChIP-seq Data
http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003342#pcbi.1003342.s019
That's all for today! I will continue to analyze the DNA methylation data.
Subscribe to:
Posts (Atom)