get the peaks that shared in multiple samples

define intronic, exonic and intergenic regions using GRanges

You might want to use the bedtools for solutions. see Dave Tang's post.

I am going to use GRanges to do exactly the same thing. However, do bear in mind that R uses 1 based coordinate system.  bed files in Granges are 1 based.  bedtools uses 0 based coordinate system for bed files. It could be very confusing and one can make 1 base-off mistake easily.

when you have a bed file in 0 based coordinate system, if you read it into R by using import function in rtracklayer package, it will convert to 1 based internally. However, if you read the 0 based bed file into R by read.table and want to convert the dataframe to GRanges, do not forget to add 1 to the start before converting to GRanges.

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.

mysql to get all the disease associated SNPs

I want to find SNPs occur at the transcription binding sites which may potentially affect the binding of the transcription factor.

I need  to have a list of SNPs that are associated with diseases first.
people usually go to dbSNP http://www.ncbi.nlm.nih.gov/SNP/
but the UCSC snp138 table hg19 is much better accessible in terms of parsing.

See posts here: http://www.biostars.org/p/1288/
http://www.biostars.org/p/7073/
http://www.biostars.org/p/11701/

find SNPs associated with OMIM gene  (by Pierre )

" Inspired by Khader's comment. The following mysql query for the mysql anonymous server at UCSC answers the SNPs in the OMIM genes:

mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A  -D hg18 -e '
select
   concat(left(title1,30),"..."),
   omimId,
   S.name,
   S.func,
   G.chrom,
   S.chromStart,
   S.chromEnd
from
   omimGene as G,
   omimGeneMap as M,
   snp130 as S
where
  G.name=M.omimId and
  G.chrom=S.chrom and
  S.chromStart>=G.chromStart and
  S.chromEnd <= G.chromEnd
limit 10;'

Result:

+-----------------------------------+--------+------------+--------------------+-------+------------+----------+
| concat(left(title1,30),"...")     | omimId | name       | func               | chrom | chromStart | chromEnd |
+-----------------------------------+--------+------------+--------------------+-------+------------+----------+
| Nucleolar complex-associated p... | 610770 | rs72904505 | untranslated-3     | chr1  |     869480 |   869481 |
| Nucleolar complex-associated p... | 610770 | rs6605067  | untranslated-3     | chr1  |     869538 |   869539 |
| Nucleolar complex-associated p... | 610770 | rs2839     | untranslated-3     | chr1  |     869549 |   869550 |
| Nucleolar complex-associated p... | 610770 | rs3196153  | untranslated-3     | chr1  |     869586 |   869587 |
| Nucleolar complex-associated p... | 610770 | rs1133980  | untranslated-3     | chr1  |     869614 |   869615 |
| Nucleolar complex-associated p... | 610770 | rs28453979 | untranslated-3     | chr1  |     869781 |   869782 |
| Nucleolar complex-associated p... | 610770 | rs61551591 | intron,near-gene-3 | chr1  |     870079 |   870080 |
| Nucleolar complex-associated p... | 610770 | rs3748592  | intron,near-gene-3 | chr1  |     870100 |   870101 |
| Nucleolar complex-associated p... | 610770 | rs3748593  | intron,near-gene-3 | chr1  |     870252 |   870253 |
| Nucleolar complex-associated p... | 610770 | rs74047418 | missense           | chr1  |     870364 |   870365 |
+-----------------------------------+--------+------------+--------------------+-------+------------+----------+

however the table is not available any more in the UCSC databases.
instead you should do:
also  by Pierre 

1) Register an access to the FTP site of omim: http://omim.org/downloads and download mim2gene:

$ curl -s  "ftp://anonymous:xxxxxxx@xxxxx.edu/OMIM/mim2gene.txt" | head
# Mim Number    Type    Gene IDs    Approved Gene Symbols
100050  phenotype   -   -
100070  phenotype   100329167   -
100100  phenotype   -   -
100200  phenotype   -   -
100300  phenotype   100188340   -
100500  moved/removed   -   -
100600  phenotype   -   -
100640  gene    216 ALDH1A1
100650  gene/phenotype  217 ALDH2

get a list of the gene symbols:

~$ curl -s  "ftp://anonymous:xxxxx@xxxxxx.edu/OMIM/mim2gene.txt" |\
   egrep -v "#" | cut -d '  ' -f 4 | egrep -v '^\-$' |\
   sort | uniq > list1.txt

2) get your list of SNP associiated to the gene symbol. Something like:

mysql -N --user=genome --host=genome-mysql.cse.ucsc.edu -A  -D hg19 -e 'select  distinct
  G.geneSymbol,
  S.name
from snp132 as S,
kgXref as G,
knownGene as K where
    S.chrom=K.chrom and
    S.chromStart>=K.txStart and
    S.chromEnd<=K.txEnd and
    K.name=G.kgId 
    /* AND something to restrict the result to YOUR list of SNPs or gene */
' | sort -t '    ' -k1,1 > list2.txt

3) use unix join to join the two lists:

join -1 1 -2 1 list1.txt list2.txt

you should get a list with two columns: the OMIM gene and your SNP.

I want to get a list of SNPs that are clinical associated. description of the snp138 table
http://genome.ucsc.edu/cgi-bin/hgTables
at terminal:

mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -D hg19 -e ' SELECT * FROM snp138 AS s WHERE s.bitfields LIKE  "clin%" ' > clinic_associated_SNPs_hg19.txt

remember you are on the remote server of UCSC, you can not connect the database first, and do something like:

SELECT *
FROM snp138 AS s
WHERE s.bitfields LIKE "clin%"
INTO OUTFILE '/home/tommy/clinic_associated_SNPs_hg.txt'
FIELDS TERMINATED BY '\t'
LINES TERMINATED BY '\n';

Because we do not have the write privilege
after I got the txt file, I  changed it to a bed file:

tommy@tommy-ThinkPad-T420:~$ cat clinic_associated_SNPs_hg19.txt | cut -f2-7 | head
chrom chromStart chromEnd name score strand
chr1 985954 985955 rs199476396 0 +
chr1 1199488 1199489 rs207460006 0 +
chr1 1245103 1245104 rs144003672 0 +
chr1 1265153 1265154 rs307355 0 +
chr1 1265459 1265460 rs35744813 0 -
chr1 1469330 1469331 rs145324009 0 +
chr1 1635334 1635335 rs201004006 0 +
chr1 1689555 1689556 rs207460007 0 +
chr1 1959074 1959075 rs121434580 0 +

tommy@tommy-ThinkPad-T420:~$ wc -l clinic_associated_SNPs_hg19.txt 
103174 clinic_associated_SNPs_hg19.txt

There are total 103174 SNPs in the file.

Now, I can just use bedtools to intersect the SNPs bed file with the ChIP-seq bed file generated by MACS.

tommy@tommy-ThinkPad-T420:~$ bedtools intersect -a TF.bed -b SNPs_hg19.bed -wo 

Just keep in mind the difference between the 0 based and 1 based coordinates system. See a post here:

http://www.biostars.org/p/84686/

Alternatively, you can create a database to do this kind of intersection by using Join command.

http://gettinggeneticsdone.blogspot.com/2013/11/using-database-joins-to-compare-results.html

http://www.sitepoint.com/understanding-sql-joins-mysql-database/

other tools

regulomeDB
http://regulome.stanford.edu/

a quick google
 http://www.genomatix.de/online_help/help_gems/SNPInspector.html
http://viis.abdn.ac.uk/regsnp/Home.aspx

Finally a python package for Database https://dataset.readthedocs.org/en/latest/quickstart.html

My first play with GRO-seq data, from sam to bedgraph for visualization

GRO-seq (Global run on sequencing) is a technique to capture nascent RNAs in a global manner.  Different from RNA-seq which sequences the matured mRNA to infer gene expression level, GRO-seq sequences the newly generated short RNAs to infer RNApol2 binding or the production of enhancer RNAs.

It was first carried out by John T.Lis at Cornell University. See their science paper in 2008:
Nascent RNA Sequencing Reveals Widespread Pausing and Divergent Initiation at Human Promoters
http://www.ncbi.nlm.nih.gov/pubmed/19056941

Since then, many labs use this technique to study their own subjects.
Michael G. Rosenfeld from UCSD just published a nature paper in June, 

Functional roles of enhancer RNAs for oestrogen-dependent transcriptional activation

http://www.nature.com/nature/journal/v498/n7455/full/nature12210.html

I want to have a look at the data set since I work with MCF7 breast cancer cells a lot.

They put the data in GEO http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE45822

although they have a bigwig file there which one can direct upload to UCSC genome browser http://www.biostars.org/p/42844/ or local IGV to visualize, the data were mapped to hg18.

I want to have a hg19 version. So, I went to download the raw sequence sra file

wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR817/SRR817000/SRR817000.sra

then I used the sra toolkit to dump it to fastq file and mapped the data with bowtie 

in UF HPC terminal:

module load sra

fastq-dump SRR817000.sra

the resulting fastq file is around 9Gb.

PBS file for mapping:

[mtang@dev1 MCF7_GRO_seq]$ cat bowtie.pbs 

#!/bin/bash

#

#PBS -N bowtie

#PBS -M mtang@ufl.edu

#PBS -m abe

#PBS -o bowtie.test.out

#PBS -e bowtie.test.err

#PBS -l nodes=1:ppn=8

#PBS -l pmem=3000mb

#PBS -l walltime=10:30:00

#

cd $PBS_O_WORKDIR

# Load the module for bowtie

module load bowtie

# Run bowtie

bowtie -p 8 --best  /project/bio/bowtie/hg19  -q SRR817000.fastq  -S  SRR817000.sam

Essentially, the sequences are mapped back to genomic DNA as you see.

It took me around 30 mins to finish the mapping which is pretty fast compared with RNA-seq data. 

using Tophat, same configuration as the pbs file above, it takes me 50hrs to finish mapping. RNA-seq data are bigger (~25Gb fastq) and more complicated when mapping (figure out exon-exon junctions).

Now, I have the sam file after bowtie mapping. 

Next, I want to generate a bedGraph file for visualization.

There are several ways to do it:

1. use Homer http://biowhat.ucsd.edu/homer/

I strongly recommend you read the tutorials to understand the various sequencing techniques. Although the formatting is not that great, it is still pretty informative.

first make Tag directory 

makeTagDirectory MCF7_gro_seq_ctrl/  SRR817000.sam

then make the bedgraph file for plus and minus strand respectively:

makeUCSCfile MCF7_gro_seq_ctrl/  -strand + -o minus.bedGraph

makeUCSCfile MCF7_gro_seq_ctrl/ -strand - -o plus.bedGraph

2. use bedtools genomecov

bedtools genomecov -ibam SRR817000.sorted.bam -bg -strand + > bedtools_plus.bedGraph 

bedtools genomecov -ibam SRR817000.sorted.bam -bg -strand - > bedtools_minus.bedGraph 

the bam file needs to be sorted first with samtools

it takes much longer for bedtools compared with Homer (Homer generates the tagDirectory for later quick access). gzip the bedGraph file to .gz

all the files I have:

[mtang@dev1 MCF7_GRO_seq]$ ls -sh .

total 19G

 87M bedtools_plus.bedGraph.gz  4.0K fetchChromSizes     16M plus.bedGraph.gz  1.5G SRR817000.sorted.bam

4.0K bowtie.pbs                 4.0K hg19.chrom.sizes   1.9G SRR817000.bam     1.2G SRR817000.sra

3.1M bowtie.test.err            4.0K MCF7_gro_seq_ctrl  8.6G SRR817000.fastq

   0 bowtie.test.out             14M minus.bedGraph.gz  5.6G SRR817000.sam

3. Use HTSeq python package: http://www-huber.embl.de/users/anders/HTSeq/doc/tour.html#tour

I then imported the bedgraph files to IGV and had a look at my favourite gene VEGFA.

I published my first paper on it:) 

Restraint of angiogenesis by zinc finger transcription factor CTCF-dependent chromatin insulation

 http://www.ncbi.nlm.nih.gov/pubmed/21896759

The first track is the minus-strand bedgraph got from Homer

the second track is the plus-strand bedgraph got from Homer

the third track is the plus-strand bedgraph got from bedtools genomecov

I see clearly short RNAs spiking  at the promoter of VEGFa, suggesting that RNApol2 paused at promoter and does not undergo productive elongation. This gene is actively transcribed with full-length mRNA, but the promoter is still paused. ( read the 2008 science GRO-seq paper for more details of different groups of genes)

one thing I do noticed is that the bedgraph file from bedtools is bigger than the one from Homer 87M vs 16M, and if I zoom in I saw some small differences:

It looks like the bedgraph from bedtools exhibits a more detailed view of the data, but overall either bedgraph file would be fine for a quick inspection of interested genes.( I guess Homer uses some bigger bin for calculating or merge bins with small difference together to save space).

Finally, if you want to make a bigwig file for UCSC genome browser, you need to use bedGraphToBigWig  from here http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/

and follow http://www.biostars.org/p/42844/.

further readings: http://www.biostars.org/p/64495/#64680
http://davetang.org/muse/2012/03/15/visualising-rna-seq-like-data/
http://davetang.org/muse/2012/03/15/ucsc-genome-browser-custom-overlap-tracks/