HOMER

If you'd like to follow along on your own, we'll analyze the GRO-Seq data from IMR90 Fibroblast published by the Lis Lab (The first GRO-Seq publication).  To download the data, go to NCBI GEO and look up accession # GSE13518. 

When downloading published data sets from GEO, sometimes you have to carefully consider what data they have made available.  Often, the primary sequencing data is available, usually through the NCBI SRA (short read archive) or the Japanese/European equivalents.  Downloading this data and mapping it your self is the safest way to go, and often the ONLY way to go if they did not supply mapped reads in the GEO record.

However, mapping the data yourself is relatively time consuming, and many times authors post processed data containing the mapped reads.  In the case of the Core et al. study, they have provided BED alignment files for each of their replicates.  It is VERY IMPORTANT to take note of which genome they aligned their reads to.  For our purposes here, just download these files (GSM340901_lib1_aligned.bed.gz & GSM340902_lib2_aligned.bed.gz)  In fact, just "copy link location" in your web browser and use "wget <CTRL+V>" to download the files directly on the command line.  After downloading the files, unzip them


Also, later in the tutorial we'll talk about how GRO-Seq interacts with transcription factors.  Unfortunately, there isn't a lot of IMR90 transcription factor ChIP-Seq data.  However a study by Chicas et al. (2010) performed ChIP-Seq for Rb in IMR90 cells, which will work for our purposes.  You can download their mapped reads from GEO # GSE19898 and follow the ChIP-Seq tutorial to process the files.  Alternatively, I've supplied the peak file and UCSC genome browser file for pRb in senescent cells for use in the tutorial (pRb peak file, UCSC file hg18).

Pool the reads from both experiments into a single Tag Directory using makeTagDirectory (more details here).  In the case of technical replicates (i.e. runs from the same sequencing library), it is advisable to always "pool" the data.  If they are biological replicates, it is often a good idea to keep them separate for and take advantage of their variability to refine your analysis.  For some types of analysis, such as transcript identification, it is a good idea to create a single META-experiment that contains all of the GRO-Seq reads for a given cell type.  This will provide increased power for identifying transcripts de novo.  For our example:


This command will run out with something like this:

Every time you get a new data set, there are a couple things you want to keep your eye on.  One is the average tags per position, which describes how "clonal" the experiment was.  For RNA-Seq, you expect a higher level of clonality, but with GRO-Seq, reads should be found all over the genome and this value should be fairly low (i.e. less than 2), unless you sequenced a ton of reads.  Another parameter to check for GRO-Seq is the Same/Diff fold enrichment, which should be high.  If HOMER doesn't guess the sample is strand specific RNA-Seq, there may be problem with how the data was mapped or with how the file is formatted.  Finally, if the GC% distribution of the fragments is extremely different from the expected value (i.e. >10%), you might want to be very careful proceeding.

We can confirm the GC-distribution by opening the "tagGCcontent.txt" and "genomeGCcontent.txt" file in the IMR90-GroSeq directory and plotting together in EXCEL


It's also a good idea to check the nucleotide frequencies as a function of distance along the reads.  Opening "tagFreqUniq.txt" and graphing with EXCEL:


As with almost any sequencing technique, there are extreme nucleotide preferences at the 5' region where the reads at ligated to adapters.  There is a little bit more G than C and T than A (normally G/C and A/T should track with one another), but this imbalance reflects levels of nucleotides downstream of the TSS minus the GC-enrichment and probably reflects natural nucleotide imbalance found in transcribed regions (graph below made with annotatePeaks.pl tss hg18 -size 2000 -hist 1 -di > output.txt)


Also, examination of the autocorrelation data (open "tagAutocorrelation.txt" with EXCEL) yields a very typical GRO-Seq pattern.  With any strand-specific RNA sequencing method, read are predominately on the same strand relative to one another.  You'll notice a slight bump in the reads found on different strands - this is a manifestation of divergent transcription in the autocorrelation plot.

To visualize GRO-Seq experiments in the UCSC Genome Browser, we'll run the makeUCSCfile command (more info here).  Since GRO-Seq is strand specific, we need to specify options to ensure it is visualized on separate strands.  For our example:


To make the resulting file only 50 MB (default), HOMER will randomly remove 60% of the data points (replaced with their average).  You may want to try:

• Changing "-fsize <#>" to a very high number (i.e. 10e10) so that it does not remove any points).  However, there is a greater chance UCSC will time out when trying to upload the larger file, so you'll just have to try it out and see what happens.
  
• Make two separate files using "-strand +" and "-strand -" (save them to different file names).
• Make a bigWig file using the "-strand" option (covered here).
  

The other parameter to play around with is the fragment length ("-fragLengh <#>").  This will adjust how much coverage each tag generates in the output.

After you get the track you want, upload it to the UCSC Genome Browser.  First step - select the hg18 version of the human genome.  This is important - you always want to make sure you have the correct genome version!!!  Then click on "add custom tracks", and upload the "IMR90-GroSeq.ucsc.bedGraph.gz" file in the "IMR90-GroSeq/" directory.  Below is an example of the region at "chr14:73,956,837-74,768,357" (colors may differ since they are selected randomly)



GRO-Seq data should "coat" annotated gene bodies as you see above.  Usually there is a high number of reads right at the TSS where RNA polymerase II is frequently found "paused".  There will also be places outside of gene bodies with GRO-Seq read coverage, such as near transcription factor binding sites.  Below is an example of a distal pRb binding site (Not near a TSS, on chr12) with the typical GRO-Seq signature that is near transcription factor binding sites (pRb UCSC Track).

GRO-Seq reads have a special distribution near regions of transcription initiation (e.g. near TSS).  Transcription proceeds divergently from the TSS in each direction.  To create a histogram of GRO-Seq reads near the TSS, use the annotatePeaks.pl program in histogram mode (add "-hist <#>").


This will produce a histogram text file from -2000 to 2000 at 10 bp resolution (change with the "-hist <#>" parameter).  Graphing the 3rd and 4th columns (that separate reads by strand) with EXCEL:

One problem with GRO-Seq (or any RNA-Seq) data is that there are some positions that have a very high density of reads, such as snoRNAs/rRNA contamination.  These sites will appear as "spikes" in the histogram.  To smooth this out, you can add the option "-pc <#>" (e.g. "-pc 3") to limit the number or reads considered at each unique position to # (or 3, in this example).

GRO-Seq is also really cool because you can see RNA species that normally aren't visible using other methods.  For example, GRO-Seq readily detects eRNAs, or RNAs originating from distal enhancer regions.   To make a histogram of GRO-Seq data near a transcription factor, all we need to perform the histogram analysis around a set of peaks.  There is a trick though... we want make sure we are away from the TSS of any known genes.  To find the "distal" peaks, you can run annotatePeaks.pl -or- you can take advantage of the getDistalPeaks.pl program that automates that process.  Be default, getDistalPeaks.pl find peaks at least 3kb from an annotated TSS (change with "-d <#>").

Try the following (pRB peak file):


This will make a histogram around the distal pRb peaks.


However, many of the peaks are found within introns, so the level of background is relatively high.  getDistalPeaks.pl has additional options to isolate only intergenic peaks ("-intergenic").



One more optimization for this process is to remove transcription factor peaks near the TTS where RNA polymerase often continues for several kB past the end of the gene.  Removal of peaks from this region further reduce the background caused by read-through trancripts (add "-noTTS" option remove peaks within 10kb of the TTS).



Placed on top of one another you can get a better sense for how the "upstream" read-through signal is smaller using more refined lists:

To find transcripts directly from GRO-Seq, use the findPeaks command:

GRO-Seq analysis does not make use of an control tag directory

Basic Idea behind GRO-Seq Transcript identification

Using uniquely mappable regions to improve results

GRO-Seq analysis output