HOMER

Software for motif discovery and ChIP-Seq analysis

Next-gen Sequencing Analysis: Creating a "Tag Directory" with makeTagDirectory

Input Alignment Files

• BED format
  
• SAM format
• BAM format (HOMER will use "samtools view file.BAM > file.SAM" to covert to a SAM file, so "samtools" must be available)
  
• default bowtie output format
  
• *.eland_result.txt or *_export.txt format from the Illumina pipeline

Specialized Input Files

As HOMER's capabilities expand, more and more file formats are being added.  Below is the current list:

• DNA methylation data: Files for mCpG data (coming soon).
  
• Hi-C summary format: Find out more about this in the Hi-C section.
  

Paired-End Reads

Creating Tag Directories

Filtering Reads/Selecting Uniquely Alignable Reads from SAM/BAM files

The makeTagDirectory program has several specialized options for processing SAM/BAM formatted files.  A quick introduction to SAM/BAM files can be found here.  Because of the various ways aligner create SAM/BAM files, it is often easiest to have HOMER filter the reads that are not uniquely alignable, etc.  Below is a list of the various options that are available to process SAM/BAM files:

-unique : selecting this option will cause HOMER to only keep uniquely alignable reads (this is also the default).  For HOMER to consider a read unique, it must NOT have the "secondary alignment" flag set, and it's MAPQ score must be greater than 10.  To change the MAPQ cutoff, use the "-mapq <#>" option.  If you set the option to "-mapq -1", HOMER will try to read AS:i:# and XS:i:# flags (if set) to figure out if the secondary alignment score is the same or less than the alignment.

-keepOne : selecting this will cause HOMER to keep all primary alignments, regardless of MAPQ score.  Alignments with the "secondary alignment" flag set will be discarded.

-keepAll : this will cause HOMER to keep all alignment in the SAM/BAM file.

-mis <#> : Maximum allowable mismatches allowed in the read.  Only works if the MD:Z: flag is set.

There are also ways to filter/select reads with specific properties:

-tbp <#> : only count one read per position per length of the read (i.e. if one read is 100 bp long, an another 75 bp long, but start at the same position, both with be kept).

-minlen <#> : minimum length of the read to keep (useful for small RNA processing, i.e. keeping 22-23 bp sequeces for miRNAs, 26-32 bp for piRNA, etc.)
-maxlen <#> : maximum length of read to keep (use in conjunction with -minlen <#>)

-filterReads <seq> <offset>  <keep | remove> : Advanced read filtering based on the presence of a specific oligo sequence in the genome (REQUIRES "-genome <genome>").  For example:

For example: -filterReads GATC -4 keep
This will only keep reads that contain a GATC sequence 4 bp upstream from the start of the read.

Interpreting the correct strand

By default, HOMER will interpret each read as it is presented in the alignment file.  This can differ from the expected results if you have paired-end sequencing from a strand-specific RNA-Seq experiment.  It can also result in reads aligning to the opposite strand if a 'first-strand' library synthesis protocol is used, such as methods that use the deoxy-UTP incorporation.  To control how HOMER interprets the strands you have two choices: Either use "-strand -" when running programs like annotatePeaks.pl when you would normally use "-strand +", or you can flip the strand during the makeTagDirectory step:

-flip : (will flip the strands of each read in the alignment file)
-sspe : (for Strand Specific Paired End sequencing - this will flip the 2nd read of the mate-pair to match the first read).

For example, if you're analyzing data from Illumina's newer strand specific RNA-Seq kits, you'll probably want to use "-flip" and also "-sspe" if you have a paired-end experiment.

What does makeTagDirectory do?

Basic Quality Control Analysis:

The following 4 basic quality control results are produced by default and are relatively inexpensive in terms of processing power to create while paring alignment files into Tag Directories, and do not require any extra information to produce.  These files are meant to be opened with a text editor or graphed using EXCEL or similar program.  More detailed descriptions of these files and how to interpret them are found in the sequencing technique-specific tutorials.  Also, if you want to rerun QC analysis on an existing Tag Directory, run the command using the "-update" option:

makeTagDirectory ExistingTagDirectory/ -update [additional options...]

Basic Tag information

tagInfo.txt - Contains basic configuration information, such as the total number of reads, the total number of unique positions with aligned reads (by genome and chromosome), and various other statistics.  One of the more important parameters is "fragmentLengthEstimate=##", which provides an estimate of the length of fragments used for sequencing.

Read Length Distribution

tagLengthDistribution.txt - File contains a histogram of read lengths used for alignment.

Clonal Tag Distribution

tagCountDistribution.txt - File contains a histogram of clonal read depth, showing the number of reads per unique position.  If an experiment is "over-sequenced", you start seeing the same reads over and over instead of unique reads.  Sometimes this is a sign there was not enough starting material for sequencing library preparation.  Below are examples of ideal and non-ideal results - in the case of the non-ideal experiment, you probably don't want to sequence that library anymore.

Autocorrelation Analysis

tagAutocorrelation.txt - The autocorrelation routine creates a distribution of distances between adjacent reads in the genome.  If reads are mapped to the same strand, they are added to the first column.  If adjacent reads map to different strands, they are added to the 2nd column.  The results from autocorrelation analysis are very useful for troubleshooting problems with the experiment, and are used to estimate the fragment length for ChIP-Seq and MNase-Seq.  The fragment length is estimated by finding the position were the "opposite strand" distribution is maximum.  HOMER will use this value as the fragment length unless overridden with the option "-fragLength <#>".

Different types of experiments (i.e. ChIP-Seq vs. DNase-Seq) produce different looking autocorrelation plots, and more detailed discussion of these differences can be found in the individual tutorials.  Below is an example from a successful ChIP-Seq experiment:

HOMER also uses the autocorrelation results to guess what type of experiment you conducted.  It computes 3 statistics:

• Same strand fold enrichment: Enrichment of reads on the same strand within 3x the estimated fragment length
• Diff strand fold enrichment: Enrichment of reads on different strands within 3x the estimated fragment length
• Same / Diff fold enrichment: Difference between enrichment of reads on the same strand or different strands

Below is a schematic to visualize how these are calculated (keep in mind that the background is calculated out to +/- 2kb):

Depending on the value of the autocorrelation quality statistics, HOMER will guess what your experiment is:

• If the Same/Diff Fold Enrichment is > 8 fold, it's a great chance the sample is strand-specific RNA.  If HOMER decides the sample is probably RNA due to the difference between same strand and different strand numbers, it will automatically set the estimated fragment length to 75 bp.  This is because it is difficult to estimate the fragment length for RNA/GRO-seq.  To manually set the fragment length, use "-fragLength <#>"
  
• If both the "Same Strand Fold Enrichment" and "Diff Strand Fold Enrichment" are both greater than 1.5 fold, there is good chance you're looking at a working ChIP-Seq experiment.

This is meant as immediate feedback, and not a definitive declaration of the quality or type of your data.  We've found it useful to help identify wrong annotations in sample IDs, for example, and helped give us an idea about how well an experiment worked. Good ChIP-Seq antibodies give "Same Strand Fold Enrichment" values well above 5.0 for transcription factors.  However, in may cases good results can be extracted out of experiments with lower enrichments, even those that yield a "Guessing sample is ChIP-Seq - may have low enrichment with lots of background" message.

Sequence Bias Analysis:

Invaluable information about a sequencing experiment can be found by examining the relationship between sequencing reads and the genomic sequence the came from.  Sequence bias analysis is not performed by default.  To perform this analysis, you must provide the genome and the "-checkGC" option.

makeTagDirectory <Output Directory Name> [options] -genome <genome> -checkGC <alignment file1> ...
i.e. makeTagDirectory Macrophage-PU.1-ChIP-Seq -genome mm9 -checkGC pu1.alignment.bed

This analysis will produce several output files in addition to the basic quality control analysis described above.

An important prerequisite for analyzing sequence bias is that the appropriate genome must by configured for use with HOMER.  In version v3.1, HOMER now handles custom/arbitrary genomes.  Instead of intalling/configuring a genome, you can specify the path to a file or directory containing the genomic sequence in FASTA format.  The genome can be in a single FASTA file, or you specify a directory where where each chromosome is in a separate file (named chrXXX.fa or chrXXX.fa.masked).  In either case, the FASTA headers must contain the chromosome names followed by white space, i.e. ">chr blahblahblah", not ">chr1-blahblahblah", or prefereably only ">chr1".

NOTE: If using a sequencing type other than ChIP-Seq or MNase-Seq, you may want to add "-fragLength <#>" since it may be difficult for HOMER to automatically determine the size of fragments used for sequencing.

Genomic Nucleotide Frequency relative to read positions

tagFreq.txt - Calculates the nucleotide and dinucleotide frequencies as a function of distance from the 5' end of all reads.
tagFreqUniq.txt - Same as tagFreq.txt, however individual genomic positions are only counted once.  If 10 reads mapped to the same position, those nucleotide counts would be added 10 times for "tagFreq.txt", but only once for "tagFreqUniq.txt". 

By default, HOMER calculates this frequency from -50 bp to +50 bp relative to the end of the estimated fragment (e.g. not +50bp relative to the 36 bp read, but +50 bp relative to the 200 bp ChIP-fragment).  This can be changed using the options "-freqStart <#>" and "-freqEnd <#>"

Below are some examples of this file graphed with EXCEL.  One is a pretty typical ChIP-Seq file, the other Chuck wouldn't say where it came from, but he's pretty sure there might be a problem with how they performed their experiment...  Quite a bit can be learned from looking at these types of plots (see individual tutorials of additional details)

Fragment GC% Distribution

tagGCcontent.txt - Distribution of fragment GC%-content.  GC% is calculated from the 5' end position of the read to end of the "fragment", not the end of the read.  For ChIP-Seq and MNase-Seq, the fragment length is automatically estimated with pretty high confidence.  If you are using another type of sequencing or want to set this manually, be sure to specify the desired fragment length manually ("-fragLength <#>").

genomeGCcontent.txt - Distribution of fragment GC%-content at each location in the genome (i.e. expected distribution)

These are very important files to check for any sequencing experiment.  Due to the numerous steps of library preparation involving amplification, size selection and gel extraction, it is very easy for the average GC%-content of the sequencing library to "shift".  This can be a disastrous problem.  Consider the following:

The problem with a GC% shifted sample is that even if the sample is random sequence, you will start to show "enrichment" at places with high GC-content in the genome, such as at CpG Islands.  This is unfortunate because most GC-rich areas are at transcription start sites, which might make you think the experiment worked, when in reality the sample was boiled instead of placed in the freezer.  See below to learn about GC-normalization.

Sequence Bias GC normalization:

HOMER has built in routines for normalizing GC-bias in a sequencing experiment.  In general, normalizing for GC-content is tricky.  First, you need to decide if it is worth normalizing at all.

When should I normalize for GC-content in my sequencing experiment?

First you must be able to identify that there is a GC-bias problem.  If the average GC% is within 5% of the expected genomic average, you're probably ok.  If you are comparing several experiments of similar type, and they all have the same approximate GC-bias, you're also probably better off just comparing the experiments as they are.  However, if one of your replicates is much more GC-rich than another, similar sample, you may want to consider normalizing it.

For some experiments, it tough to know what the expected GC-content should be.  For example, if sequencing H3K4me3 ChIP-Seq data, which is typically found near CpG Islands, you may expect the GC% to be higher.  However, even with H3K4me3 ChIP-Seq, most ChIP experiments are not very efficient, and most of the DNA being sequenced is background.  As a result, the average GC content should not be too far off the expected genomic average.  You can always "try" normalizing.

If the GC-bias is severe, you might be better off repeating the experiment.  Normalization essentially destroys information, so repeating the experiment (or at least the sequencing library preparation) may be a good move.  When re-preping the sample, try amplifying less and when extracting DNA from gels, dissolve the gels at room temperature.

Normalizing tags for GC-bias

To normalize for GC-bias, run the same makeTagDirectory command, but you must specify the genome (i.e. "-genome hg18") and a target normalization file ("-normGC <filename>").  To normalize to the expected genomic average, use "-normGC default".

makeTagDirectory <Tag Directory> -genome <genome> -normGC <GC profile file profile> -iterNorm 0.01 <alignment file 1> [alignment file 2] ...

i.e. makeTagDirectory Macrophage-PU.1-GC -genome mm9 -normGC default -iterNorm 0.01 aligned.reads.bed

The new GC normalization procedure attempts to be conservative by only normalizing by 'removing' reads from the experiment by down weighting the read counts associated with each position.  It is dangerous to mathematically "add reads" to the experiment - this is sort of like creating information which we don't really have, particularly for later if you are doing count-based differential expression.  HOMER will remove iteratively remove reads from GC ranges that are overrepresented until the overall distribution of GC is with the indicated percentage of the "-iterNorm <#>" parameter.  Of course, it doesn't actually remove the reads, just reduces them to fractional values.

The alternative/original GC normalization procedure is actually very simple.  HOMER calculates the distribution of fragment GC%-contents, then for each range of GC%, normalizes the tag counts to the expected distribution.  If your sample is more GC-rich that expected, reads in GC-rich regions will be reduced to a fractional value (limited by "-minNormRatio <#>").  Likewise, if reads are in AT-rich regions, their values will be increased (limited by "-maxNormRatio <#>").  HOMER is happy to deal with fractional tag values - there is no reason why a sequencing read can't be counted as a "half-a-read", for example.  The biggest problem problem with increasing the value of a read is that it is akin to creating information that doesn't exist, so HOMER puts a tight cap on that adjustment to only 2-fold (change with "-maxNormRatio <#>").  The resulting normalization information is recorded in the tagGCnormalization.txt file.  To enact this form of GC normalization, omit the "-iterNorm <#>" parameter.

Often, the expected genome GC-content is not appropriate.  To normalize to a custom GC-profile, provide a file after -normGC option ("-normGC <filename>").  This file should be formatted just like the tagGCcontent.txt files found in the directories.  In fact, the idea is that you can normalize one experiment to another by providing the experiments tagGCcontent.txt file as the argument for -normGC.  One potential strategy is to combine several experimental files across different experiments into a single tag directory, and then use the tagGCcontent.txt file from this directory to normalize each of the other experiment, normalizing each experiment to the average from all of the experiments.

The difference between the two normalization strategies can be seen in the resulting read counts of the experiment (the -iterNorm <#> method essentially removes reads until it can match the distribution):

Command line options of makeTagDirectory command:

Next: Creating UCSC Genome Browser visualization files