hiddenDomains Manual, Tutorial and Troubleshooting Guide

hiddenDomains is a suite of programs used to identify significant enrichment of ChIP-seq reads that span large domains, like HK27me3. The input data can be in BAM format, or in a tab-delimited "reads per bin" format described below. The output is a BED formatted file the lists the enriched domains and their posterior probabilities.

hiddenDomains requires samtools, so if you don't already have that installed, you'll need to do that.

This manual and tutorial will teach you how to use the hiddenDomains suite of programs. In it we will analyze a mouse ChIP-seq dataset for H3k27me3

Quickstart Tutorial

1. Download and unpack hiddenDomains

   If you haven't already done this, download the latest version of hiddenDomains from the Sourceforge website: https://sourceforge.net/projects/hiddendomains/

   Now unpack it

   
   shell$ tar -xzvf hiddenDomains.VERSION.NUMBER.tar.gz
   

   Where VERSION and NUMBER represent the release that you downloaded.

2. Install R and required R packages (if necessary)

   Download and install R if you don't already have it.

   If you don't already have them, you will need to install two hidden Markov model libraries:

   • depmixS4
   • HiddenMarkov
   You can do this by starting R on the command line:

   
   shell$ R
   

   Or, if you prefer RStudio, you can use that.

   Regardless of how you started R, type the following commands to install the packages.

   
   > install.packages("depmixS4")
   > install.packages("HiddenMarkov")
   

3. Analyze the Data

   Once you have unpacked hiddenDomains, the next step is to analyze it.

   hiddenDomains comes with two example data files:

   • k27me3_chr6.bam, An H3K27me3 ChIP-seq dataset that we will use to identify enriched domains.
   • input_chr6.bam, An "Input" ChIP-seq dataset that we can use to limit the number of false positives by ruling out domains that are enriched regardless of the antibody.
   Both of these files only contain reads that mapped to chromosome 6 in order to limit the size of the files and to get you up and running as quickly as possible.

   To analyze the data with hiddenDomains.

   
   shell$ cd hiddenDomain.VERSION.NUMBER
   shell$ ./hiddenDomains -g ChromInfo_mm9.txt -b 200 -t k27me3_chr6.bam -c input_chr6.bam -o results
   

   hiddenDomains has three required parameters:

   • -g ChromInfo_FILE.txt :A tab delimited text file that lists the chromosome names (these have to match the chromosome names in the BAM files) and their sizes. hiddenDomains comes with ChromInfo.txt files for MM9, MM10, HG18 and HG19, and is easy to make a custom file for other genomes.

     To make your own ChromInfo.txt file, just create a tab delimited file with the chromosome names in the first column and their sizes in the second column. Here is an example ChromInfo.txt file for a genome with 3 chromosomes:

         chr1    197195432
         chr2    181748087
         chr3    159599783    
     

   • -t CHIP_SEQ.bam :The ChIP-seq dataset you want to analyze. The -t stands for "treatment". NOTE: While you can run hiddenDomains without a control (or "input") file, we strongly recommend using one. This will help eliminate false positives in the output.

   • -o OUTPUT_FILE_PREFIX :The prefix that you want all of the output filenames to have. hiddenDomains creates four or five output files. These are:
     1. OUTPUT_FILE_PREFIX_treatment_bins.txt :This file contains "binned" read counts.
     2. OUTPUT_FILE_PREFIX_control_bins.txt :This file contains "binned" read counts. hiddenDomains creates this only if you provide a control dataset.
     3. OUTPUT_FILE_PREFIX_domains.txt :This file contains the enriched domains and peaks, along with their posterior probabilities.
     4. OUTPUT_FILE_PREFIX_vis.bed :A BED file that is ready to be uploaded to the UCSC genome browser.
     5. OUTPUT_FILE_PREFIX_analysis.bed :A BED file that has all consecutive enriched bins merged into a single line, making it ideal for analysis with bedtools.

   In addition to these required parameters, hiddenDomains accepts many optional parameters (described in detail in the manual). The optional parameters that are used here are:

   • -b BIN_WIDTH :The width of the bins for analysis. The default size is 1000bp.
   • -c INPUT_SEQ.bam :The "input" dataset you want to analyze. The -c stands for "control".

   NOTES:

   1. Additional parameters and usage information can be found by entering hiddenDomains -h

   2. Depending on the computer, this command should run in less than 5 minutes with ample notices printed to the terminal while it works.

   3. If you get an error that hiddenDomains could not run, you either need to make it executable, with the command

      shell$ chmod 755 hiddenDomains

      and then re-run the command above, or run it directly with Perl using the following command

      shell$ perl ./hiddenDomains -g ChromInfo_mm9.txt -b 200 -t k27me3_chr6.bam -c input_chr6.bam -o results

   To verify that hiddenDomains has worked correctly, we have provided bedgraph files for the original data and you can upload these to the UCSC Genome Browser along with the BED files that you just created. You can then peruse the results at your liberty, however, there are two locations that we would like to point out that illustrate the value of having a control dataset.

   These are the files you need to upload to the UCSC Mouse mm9 Genome Browser:

   1. input_chr6.bedgraph.gz - this will display the normalized reads from the "input". The reads should appear uniformly distributed. You may want to set the track "Vertical viewing range" from 0 to 200. You do this by clicking on the grey bar on the left side of the track, setting the Vertical viewing range "min" and "max" values to 0 and 200, and then selecting "use vertical viewing range setting" form the "Data view scaling" pulldown menu.
   2. k27me3_chr6.bedgraph.gz - this will display the normalized reads - this will display the normalized reads from the K27me3 ChIP-seq experiment. The reads should appear piled up near genes that are repressed. You may want to set the track "Vertical viewing range" from 0 to 200.
   3. k27_hiddenDomain_no_control.bed
   4. k27_hiddenDomain.bed

   Once you finish uploading the files, you can point the Genome Browser to the Hox gene cluster at chr6:52,000,664-52,512,712. Here you see that both methods called the rather broad domain of H3k27me3 reads that cover this cluster of genes.

   

   chr6:103,414,566-103,717,194 highlights the importance of using a control dataset to limit false positives. In this example there is a location with a lot of reads in both the K27me3 dataset and the input dataset. Without using the input dataset as a control, hiddenDomains calls a peak at this location. However, when run with the input dataset as a control, hiddenDomains does not call this peak.

   

Manual

hiddenDomains can be run two different ways. In the tutorial, we used hiddenDomains to do three things:

1. Bin the reads

2. Look for domains

3. Generate BED files

However, for more control, you can do these things separately. The individual programs have additional options. For example, you can specify to only analyze certain chromosomes, or you can customize the track line in the "vis" BED file.

If you want to run the programs separately...

1. binReads.pl bins reads
2. hiddenDomains.R looks for domains
3. domainsToBed.pl and domainsMergeToBed.pl create BED files.

binReads.pl

Usage: binReads.pl [options] bamFile.bam > binned_reads.txt

binReads.pl takes a BAM formatted file of aligned reads and outputs the number of reads that overlap each bin, where a "bin" is an equally sized, non-overlapping portion of a chromosome.

Options

domainsMergeToBed.pl

Usage: domainsMergeToBed.pl [options] domainFile.txt > domainFile.bed

domainsMergeToBed.pl takes the output from hiddenDomains.R and merges all consecutive domains with a posterior greater than a threshold (default is 0; all domains are merged) into a single domain and generates track of enriched domains for the UCSC Genome Browser.

Options

    chr1    197195432
    chr2    181748087
    chr3    159599783    

    >= 0.9
    0.9 > posterior >= 0.8
    0.8 > posterior >= 0.7
    0.7 > posterior >= min posterior for significance

domainsToBed.pl

Usage: domainsToBed.pl [options] domainFile.txt > domainFile.bed

domainsToBed.pl takes the output from hiddenDomains.R and generates track of enriched domains and their posterior probabilities for the UCSC Genome Browser.

Options

    chr1    197195432
    chr2    181748087
    chr3    159599783    

hiddenDomains

Usage: hiddenDomains [options] -g ChromInfo -t TreatmentReads -o OutputPrefix

Options

    chr1    197195432
    chr2    181748087
    chr3    159599783    

    >= 0.9
    0.9 > posterior >= 0.8
    0.8 > posterior >= 0.7
    0.7 > posterior >= min posterior for significance

hiddenDomains.R

Usage: hiddenDomains(treat.bin.file="TREATMENT_bins.txt", control.bin.file="CONTROL_bins.txt", out.file.name="OUTFILE.txt", chr.names=["mouse"|"human"|c("chr1","chr2", "etc"]))

hiddenDomains.R takes the output from binReads.pl and attempts to identify enriched domains using Hidden Markov Models (HMMs).

Options

Troubleshooting Guide

• If the hiddenDomains.R script fails to find peaks/domains on a chromosome, but finds them on others...
  1. Consider tweaking the max.read.count and possibly the min.read.count parameters. I've found that setting max.read.count to 10 will fix things a lot of the time. A more rational method for figuring out the optimal parameter is to look at normalized bigWig tracks in the UCSC genome browser to get a sense of how tall the peaks really are and use that number.
  2. Try using a larger bin size (i.e. use -b 1000) for binReads.pl and domainsToBed.pl and re-running everything.