-
Notifications
You must be signed in to change notification settings - Fork 22
Poster Example 1
Andy Pohl edited this page Nov 21, 2013
·
4 revisions
Install bwtool:
$ git clone https://github.com/andypohl/bwtool.git
$ cd bwtool
$ ./configure
$ make
$ make install
Obtain the ENCODE data from here.
$ curl -O ftp://ftp.sanger.ac.uk/pub/gencode/release_18/gencode.v18.annotation.gtf.gz
$ curl http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistoneMcf7H3k27acUcdSig.bigWig -o H3K27ac.bw
$ curl http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistoneMcf7H3k09me3UcdSig.bigWig -o H3K9me3.bw
$ curl http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistoneMcf7H3k27me3bUcdSig.bigWig -o H3K27me3.bw
$ curl http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistoneMcf7H3k36me3bUcdSig.bigWig -o H3K36me3.bw
$ curl http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistoneMcf7InputUcdSig.bigWig -o input.bw
Create a file from the GENCODE annotation containing merely the basic coordinates of the protein-coding genes i.e. chromosome, TSS, TTS, name, and strand. Proper bed files have a score as the fifth field and strand as the sixth field. Scores are not used by bwtool so can be set to zero:
$ zcat gencode.v18.annotation.gtf.gz \
| awk 'BEGIN{OFS="\t"}{if($3=="gene" && $20=="\"protein_coding\";"){print $1, $4-1, $5, $18, "0", $7}}' \
| sed 's/\"//g;s/;//' | sort -k1,1 -k2,2n > gencode_pc.bed
Run bwtool:
$ bwtool agg 5000:5000 -starts -long-form=TSS,H3K27ac,H3K27me3,H3K36me3,H3K9me3,input gencode_pc.bed \
H3K27ac.bw,H3K27me3.bw,H3K36me3.bw,H3K9me3.bw,input.bw plot.txt
Make a plot with R:
> library(ggplot2)
> library(extrafont)
> plots <- read.table(‘plot.txt’)
> colnames(plots)<-c('Feature','ChIP','Position','Signal')
> figure<-ggplot(plots, aes(x=Position,y=Signal,color=ChIP)) + geom_line(size=1.5) + theme_minimal() + scale_color_manual(values=rev(brewer.pal(5,"Blues"))) + theme(text=element_text(family="Helvetica")) + ylab("ChIP Signal")
> figure
> ggsave(“figure.png”)