HOMER

makeTagDirectory HiC-Unprocessed reads1-1.sam,reads1-2.sam reads2-1.sam,reads2-2.sam -tbp 1

Then experiment with different processing parameters (remember to include "-update" - this tells the program not to look for new data files but reprocess the existing tag directory)


This is a more efficient way to approach Hi-C processing if you're not sure what the data looks like yet.

Use makeTagDirectory to process paired-end sequencing into a tag directory for use with HOMER programs.  For Hi-C, you want to map the reads independently.  This is different than most "paired-end" alignments (such as for RNA-Seq) that constrain the relative position of the fragments.  For Hi-C these assumptions do not hold.  I'd recommend mapping the same way you would ChIP-Seq data, allowing only unique genomic matches. You may also want to check for restriction sites in your reads - sequence on either side of the restriction site might map to different locations, so it might be wise to trim the sequence after the restriction site to help with mapping (i.e. use "homerTools trim -3 AAGCTT reads.fq" to trim the reads before mapping)

To create a paired-end tag directory, you must run makeTagDirectory by entering the two halves of the mapped reads separated by a comma, without spaces.  For example:


HOMER matches reads between alignment files by matching their names, but to work with the original bowtie output files for Illumina sequencing (which often contains a 0 or 1 at the end to delineate the read source), it will need to remove the last character before trying to match.  In this case, add "-illuminaPE" to the end of the command.  If you download reads from the SRA they often have the same name for each paired end read - in these cases you don't want to remove the last character, so do NOT add "-illuminaPE" to the end of the command in that case.

Another typical option to add for Hi-C data is "-tbp 1", which will only consider read pairs with the exact same ends once.  Given the search space in Hi-C, you would never expect to see reads-pairs start from the exact same places twice in the data (at least in the case of sonication).  If you do, they are likely clonal and a result of over-sequencing your library or some other artifact.

On a side note, the HOMER tag directory created for paired-end reads is very similar to the directories created when running makeTagDirectory for single-end reads.  The paired end reads are stored in *.tags.tsv files (extended format to include the pair), only each read will get 2 entries in the appropriate *.tags.tsv so that the reads can be indexed on each chromosome in one file. A flag is also set in the tagInfo.txt file to inform HOMER to treat the experiment as a paired-end file.  Many of the programs in HOMER, such as makeUCSCfile, will treat paired-end tag directories as "single-end" directories to create output.  This can be useful if you want to visualize read coverage in the UCSC Genome Browser or want to find peaks in the Hi-C read coverage with findPeaks.

One important note: the tagInfo.txt file will report the "Total Tags" value as 2x the total number interactions.  This is because HOMER is storing each PE tag twice - the software later adjusts this number to reflect this redundancy, but keep this in mind when interpreting this number.

When running the makeTagDirectory program, several standard quality control analysis files will be created that resemble the analysis carried out for single-end reads.  You can also check for GC bias by specifying the genome (i.e. "-genome hg19") and "-checkGC" at the command line (analysis will be conducted by treating the samples as "single" reads).

The following two Hi-C/paired-end read specific QC files will also be produced by default:

petagLocalDistribution.txt
petagDistDistribution.txt

One major problem that can arise during Hi-C is that "background" ligation events may occur.  These are ligation events that occur after sonication [I think] as they seem to occur between random fragments of DNA with no regard to restriction site used for the Hi-C assay.  In theory, all reads should be in the vicinity of the restriction site used for the assay, and if they aren't, there is a good chance they are noise (more later on this).  To see the distribution of reads around your restriction site, use the option "-restrictionSite <seq>" (i.e. "-restrictionSite AAGCTT").  You MUST specify the genome so that HOMER can figure out where the restriction sites are located (i.e. "-genome mm9").  If the restriction site has a lot of star activity, you can also add "-rsmis <#>" (i.e. "-rsmis 1") to indicate how many mismatches are tolerated in the sites.  By default it only looks for perfect restriction sites.

Adding these options will produce a file named "petagRestrictionDistribution.<seq>.mis<#>.txt" (such as petagRestrictionDistribution.AAGCTT.mis0.txt).  This file will display the distribution of reads relative to the restriction site.  Graphing it with Excel will look something like this:

The Hi-C protocol can produce sequencing libraries that contain many read pairs that are not representative of true interactions.  This includes continuous genomic fragments, self-ligation or re-ligation products, and background non-specific ligation that occurs independent of a restriction site.

-removePEbg : This will remove paired-end reads that are likely continuous genomic fragments or re-ligation events.  If this option is specified, read pairs are removed if they are separated less than 1.5x the sequencing insert fragment length.  This distance is estimated from the distribution found in the file "petagLocalDistribution.txt" described above.  If you wish to specify your own size, use "-fragLength <#>" to set this length manually.

-removeSpikes <#size> <# fold> : This option will remove PE reads that originate from regions of unusually high tag density.  These regions arise for various reasons, for example, a region was duplicated in a given cell line.  These types of anomalies tend to create artifacts in the data.  Removing these from the analysis can help clean things up.  The size parameter is the window of tag density that will be examined to find anomalies, and #-fold is how many tags over the average is needed to determine if a locus should be removed from the dataset.  I'd recommend something like "-removeSpikes 10000 5" to get ride of the regions that are really bad.  For example, this will scan all 10kb regions in the genome, calculate the average read density, and remove reads from regions that contain more than 5x the average number of reads.

Filtering Reads based on Restriction Sites