Genomatix Software GmbH - RegionMiner ChIPSeq Workflow

If the read counts in clusters for two different conditions (here called "treatment" and "control" for simplicity) are to be compared, two statistical testing methods for evaluating differences in read abundance are available:

For defining enriched and depleted clusters between two conditions or samples, the following criteria are used (parameters set by the user):

Parameters

Input

Input file(s) with read positions from ChIPSeq ("Sample"/"Treatment")

Input data are accepted as a tab delimited file in BED / bigBed file format containing the input regions specified at least by chromosome number, start position and end position (in this order).
The maximum amount of regions and their maximum length can differ for various tasks. The limits are usually shown on top of the input pages.

Within this section you can either

choose from previously uploaded BED files
or add a new bed file to the list (by clicking "Add Bed file...")

When adding a new file, a new window will open, asking you to either

upload one or several BED files from your local computer
or import a BED file from the GMS (see more details)
or import a BED file from the GGA (see more details)

For the new BED files, you will have to select the correct organism, as the organism and the genome build are associated with the BED file for future use (the default is your latest choice in the current session).
Note that BED files critically depend on the underlying genome build, which can be changed by selecting a different ElDorado version on the top right of the page before uploading a BED file. You can see the list of genomes available in ElDorado.

Note that almost all browsers have a general upload limit of 2 GB, i.e. BED files bigger than this size should be zipped before uploading from your local computer. This restriction does not apply when using the direct import from the GGA/GMS.

Optionally you can specify a name for saving uploaded BED files on the server, otherwise the name of the uploaded file will be used. If several files are uploaded, the string given here will be used as prefix for each BED file name.

If any of the regions in the input file cannot be completely assigned to the selected genome (e.g. wrong chromosome numbering or wrong positions within a chromosome), an error message will appear and the regions will be skipped. If no valid region is found in an uploaded file, the complete file will be skipped.

After one or several BED files were uploaded successfully, and after closing the popup window, the list of available BED files will be automatically updated.

Uploaded BED files can be deleted from the project anytime via the project management.

Example input:

track description="sample treatment-control analysis with 3 treatments and 3 controls"
chr1       26519270        26519623
chr1       39723904        39724119
chr2       10841542        10841853
chr2       88937859        88938309

Control file(s)

Optional Control file(s)
for differential analysis

If one or several input control files (replicates) are available, they can be uploaded here. This is an optional field, so if no control files are available the checkbox should be left blank.

The upload options for control files are hidden by default. They appear only when checking the checkbox "Use second set of input files (control files) for differential analysis".

Workflow parameters

Read Classification

When checked, a read classification is done for each BED file from the input data. The number of input reads overlapping genomic elements like exons, introns, promoters and intergenic regions will be given in the result.

Peak Finding / Clustering

This step is mandatory, as all subsequent analysis steps rely on the results of the clustering.
For the peak finding / clustering one of three different algorithms can be selected:

Depending on the user selection, further parameters appear for the respective algorithm:

NGS Analyzer parameters
Window size for clustering:	The window size used for the clustering algorithm in basepairs. The default window size is 100 bp.
Minimum number of reads per cluster	A threshold for the clustering algorithm. By default (value -1), this number is calculated from the data set applying a Poisson distribution. Otherwise, values above 3 are allowed here.

MACS parameters
Tag size	The length of the input tags/reads. By default (-1), this value is determined from the input BED file, by reading the first 100 BED regions and calculating the average region length
p-value	P-value cutoff for peak detection. the default is 1e-5
Bandwidth	This value is used while building the shifting model. Default is: 300
mfold	The upper model fold value for MACS to select the regions with a high-confidence enrichment ratio against background to build a model. The lower model fold value is automatically set to 1. If no models are found, the no-model option is used by MACS automatically.
Redundancy threshold	The number of copies of identical reads allowed in a library. Values can be 'auto', 'all' or an integer 'auto': MACS calculates the maximum tags at the exact same location based on binomal distribution using 1e-5 as pvalue cutoff 'all': all input tags are kept integer: at most this number of tags will be kept at the same location

SICER parameters
Redundancy threshold	The number of copies of identical reads allowed in a library.
Window size	Resolution of SICER algorithm. For histone modifications, one can use 200 bp. Note from the SICER manual: The choice of window size and gap size has a large effect on outcome. In general, the broader the domain, the bigger the gap should be. For histone modifications H3K4me3, W=200 and (gap = 1 window) are suggested. For H3K27me3, W=200 and (gap = 3 windows) are suggested for first try. If even bigger gap size is found to work better, you might also want to try increasing the window size eg, window size = 1K, and gap = 3 windows)
Fragment size	Is for determination of the amount of shift from the beginning of a read to the center of the DNA fragment represented by the read. FRAGMENT_SIZE=150 means the shift is 75.
Gap size	Needs to be multiples of window size. Namely if the window size is 200, the gap size should be 0, 200, 400, 600, ...
FDR	The FDR is calculated using p-value adjusted for multiple testing, following the approach developed by Benjamini and Hochberg.
E-value	nr. of islands expected in random background, only if no control data supplied Note: E-value is not p-value. Suggestion for first try on histone modification data: E-value=100. If you find ~10000 islands using this evalue, an empirical estimate of FDR is 1E-2

Replicate Treatment

When using replicate data as input, first all input data is clustered separately. Then the clusters from the different replicates from the treatment group are either merged or intersected.
To this end, first a coverage profile is created. This helps to identify contiguous clusters from the different treatment cluster files. In particular, for each position within a contiguous region the number of the replicates (from the treatment group) having this position in common is calculated (i.e. how many replicates contain a cluster covering this position).
There are two parameters to control how the merging/intersection is done:

minimum number of replicates with cluster at a position:
This parameter can be used to set the required fraction of replicates in a cluster. If for example there are 3 replicates, and the parameter is set to 60%, it is sufficient that 2 replicates have overlapping clusters. (2/3 ≈ 67%).
minimum length of common region:
Only "candidate regions" which have a common region of the specified length will be kept.

The following example illustrates the effect of the parameters. Let's assume we do have three replicates (r1, r2, r3) with overlapping regions at different positions, as illustrated by the following graphic that also includes the coverage at each position:

r1 ----------     -----          -------------------
r2       -----------    ------     -----        -------
r3    ---------------              -----        -----------------------------
   11122233332222233211011111100011333331111111133332221111111111111111111111 <- coverage profile

With the minimum number of replicates set to 60% and a mininum region length of 10 we would get the following results ('+' denotes positions where the threshold for minimum number of replicates is met):

Merging clusters:
   --++++++++++++++++-- ------   --+++++--------+++++++----------------------
   M1                   M2       M3

Intersecting clusters:
     ++++++++++++++++              +++++        +++++++
     I1                            I2           I3

For the cluster M1 we have a length of 20 nucleotides, with 15 covered by at least 2 clusters, so we do have 67% coverage for 15 out of 30 nucleotides. These regions would therefore be merged into a cluster.
For M2 we do only have a single region from r2 which has no overlap and thus fails for both thresholds.
For M3 we have a length of 44, with at least 2 common regions for 12 nucleotides, satisfying both thresholds.
For I1 we have 16 positions where both thresholds are met, therefore I1 is a cluster.
I2 is only 5 positions long, below the minimum of 10 and I2 is rejected.
For I3 there are 7 positions, which is still below the minimum, so I3 is not considered as well.

Afterwards, the reads from all input BED files, treatment and control, are distributed on the new regions (found by merging resp. intersecting the clusters from the treatment as described above). This way, for each region we find the number of reads contained in this region, for each condition and replicate. Based on these count data the chosen statistical test is applied to evaluate the difference in read abundance between treatment and control.

Differential Analysis parameters

The differential analysis parameter section will only appear, if at least one control file was uploaded in the section above.
There are three available algorithms for calculating the differential expression/enrichment values:

DESeq (recommended for replicate data, but does work on non-replicates, too)
edgeR (recommended for replicate data)
the procedure introduced by Audic and Claverie (if no replicates are available)

For a short introduction to the different methods, see above in the Introduction to the Differential Analysis.

The thresholds that define a transcript as differentially expressed (or a region as enriched/depleted) can be set here. There are two criteria, that are combined (both must be satisfied for differential expression/enrichment):

a p-value threshold for the significance of observing the detected change

Note that the p-values calculated by the different methods (DESeq, edgeR, Audic-Claverie) can differ.

Also note, that setting the p-value to 1 allows skipping of this criterium.
a threshold for the log2 fold change of expression/enrichment level
A log2 ratio of 1 is a fold change of 2; a log2 ratio of 0.585 is a fold change of 1.5; e.g. if the log2 fold change of expression/enrichment is set to ≥ 1, the expression values must go up by at least 100% to appear in the differentially expressed transcripts/enriched regions list.

The log2 fold change thresholds can be set separately for up- and down-regulation (enrichment/depletion).

Note, that by setting the log2 fold change thresholds to 0, fold changes are ignored in the analysis.

Cluster Classification

When this option is checked, a cluster classification and statistics is done. The number of clusters (defined in the previous clustering step) overlapping genomic elements like exons, introns, promoters and intergenic regions will be given in the result.

Sequence Extraction

When this option is checked, the DNA sequence of all clusters (peak regions) defined in the clustering step is extracted in FASTA format. The resulting sequence file can be downloaded from the result and be used for further downstream analysis.

TFBS Overrepresentation

When this option is checked, the cluster sequences are checked for transcription factor binding sites (TFBS) and statistics on TFBSs together with overrepresentation values and Z-scores are generated.

The binding sites are based on the matrix families as available in MatBase and searched by MatInspector

For more details please see the RegionMiner Task Overrepresented TFs.

Definition of new TFBS

If this option is checked, a subset of cluster sequences is selected and then the Genomatix program CoreSearch is used to define new common motifs from the cluster sequences.
The program automatically sorts the cluster sequences by quality, and the number of "best sequences" used for CoreSearch can be set by the user. Here, the definition of "best" cluster sequences depends on the selected clustering algorithm:

for MACS the clusters are sorted by decreasing MACS score
for NGSAnalyzer the clusters are sorted by decreasing NE-value
if a Audic-Claverie evaluation was used, the clusters are sorted by increasing p-values

The sequences actually used for CoreSearch can also be downloaded from the result.

Output

Result

Here, you can edit the default name of the result file.

Email address

Here you can choose between two methods for receiving the results:

Show result directly in browser window
In this option the URL of the result is directly shown in your browser window.
Warning: Please use this option only for analyses which can be performed in a short time.
If the analysis takes longer than the timeout of the webserver, the connection will be terminated and you will receive an error message (e.g. "The document contained no data."). In this case, the results will not be available, please restart the analysis using the option below "Send the URL of the result to".
Send the URL of the result via email
In this option an email with the URL of the results will be sent to the user provided email address, when the analysis is finished.

The results will be available for a limited time on our server. For details of how long your results will be kept please see the result-email. After that period they will be deleted unless protected in the project management!

We recommend to use the email option for the ChIPSeq Workflow!

Please note, that a job started with the "show result in browser" option, will finish correctly even if the display in the window stops, and can then be found in the result management.

Output

1. Analysis Parameters

2. Main Results, depending on user selection

Sample Read Classification and Statistics

If the read statistics was selected as parameter, a table with the number of reads from the input overlapping genomic elements like exons, introns, promoters and intergenic regions is given.
Additionally a table with the distribution of the reads on the different chromosomes of the genome is shown. The content of this table is hidden by default, but can be shown by clicking the "Show details" link in the header.

If the read classification is included in the output, detailed annotation for each read can be downloaded as a tab-separated file. For a description of the format of this file, please see the cluster classification details.

Peak Finding / Cluster Generation

Here is a short overview of what happens at different input/parameter settings:

	no control data set	one sample control set	replicates for control data
one sample data set	input data is clustered no control / peak evaluation	input data is clustered peak evaluation by original algorithm (MACS/SICER) or Audic-Claverie for NGSanalyzer only significant enriched clusters are used for further analysis	input data is clustered the first control data set is used for peak evaluation by original algorithm (MACS/SICER) or Audic-Claverie for NGSanalyzer
replicates for sample data	each replicate data set is clustered separately clusters from replicates are merged no control / peak evaluation	each replicate data set is clustered separately clusters from replicates are merged reads from sample and control data are distributed on the merged clusters peak evaluation by differential analysis (DESeq/edgeR) only significant enriched clusters are used for further analysis	each replicate data set is clustered separately clusters from replicates are merged reads from sample and control data are distributed on the merged clusters peak evaluation by differential analysis (DESeq/edgeR) only significant enriched clusters are used for further analysis

Depending on the selected peak finding algorithm and the parameters, the output can look different in this section:

For all settings, the total number of clusters found by the program is shown and the resulting clusters can either be downloaded as a BED file or saved directly to the Genomatix project management, to be used with other RegionMiner tasks. Also, a link to the complete algorithm output is given, including details as described for NGS Analyzer or MACS or SICER respectively.

If "Peak Evaluation with Audic-Claverie Algorithm" was selected, the BED file contains only those clusters, which show a significant enrichment of reads, i.e. a subset of all significant clusters. "Significant" in this context means, that the Audic-Claverie p-value is at most as high as the cut-off specified on the input page. Additionally, a tab-separated file containing the p-values and other details for each input cluster can be downloaded (details for p-value file). This file contains all significant clusters, no matter if there is an enrichment or decrease of reads.

Note that the BED file format is zero-based and half-open, whereas numbering in the tab-separated p-value file is based at 1 and includes the end position.

cluster results output

If differential analysis by DEseq/edgeR was conducted, the output will contain a table with the numbers and download links for the analyzed regions, the significant regions (both enriched and depleted), the enriched regions and the depleted regions.

cluster results output

The analysis results are available in two formats, TAB-delimited (.tsv) and BED. Each single file can be dowloaded using the link in the table, but they are also available for download in one tar-archive.

Regions analyzed: 10.region_summary.tsv/bed
Regions with significant change: 11.significant_regions.tsv/bed
Enrichment: 12.enriched_regions.tsv/bed
Depletion: 13.depleted_regions.tsv/bed

The tsv-files have the following columns:

    1  Id: region id, an integer
    2  Contig: contig/chromosome accession number
    3  Chromosome: chromosome (with leading 'chr', e.g. 'chr1', 'chrX')
    4  Strand: strand of the region, '+', '-' or '0' if not strand-specific
    5  Start pos: region start position
    6  End pos: region end position
    7  Length: region length (end - start + 1)
    8  #Read total: sum of the read counts of all replicates
    9  p-value: p-value resulting from the hypothesis test of the selected test method
                for difference in read abundance (DESeq or edgeR)
   10  adj. p-value: Benjamini-Hochberg adjusted p-value (from column 9)
   11  log2(fold-change): logarithmic fold-change in read abundance in treatment and control
                          (> 0 means an enrichment in treatment,  < 0 means a decrease in treatment)
   12  Regulation: keyword indicating increase or decrease of reads in the region;
                   'up' if log2(fold-change) > 0, 'down' if log2(fold-change) < 0, 'no' otherwise
    + read count for each replicate
    + NE value for each replicate
   
i.e. if there are X replicates, there will be 12+2*X columns

The format of the BED files is as follows:

    1  chromosome: chromosome (with leading 'chr', e.g. 'chr1', 'chrX')
    2  start: start position of the region
    3  end: end position of the region
    4  id: region id
    5  (adj.) p-value: as in tsv-file column 9 or 10, depending on the parameter setting for
                       'p-value adjustment'
    6  strand: strand of the region, '+', '-' or '0' if not strand-specific

Plots

Note:The graphics are at first depicted as smaller icons, which can be enlarged by clicking on them.

93.fold_change_plot.png
graphics in PNG-format with a scatter plot of fold-change in the read abundance in treatment versus control (y-axis) against read abundance (x-axis). For edgeR, the measure of abundance is the concentration, for DESeq, it is the base mean, i.e. the values shown in the plot can be found in 92.output_replicate_analysis columns 4 and 5. Each data point resembles a region (i.e. merged/intersected clusters from the treatment group), those showing significant difference in read abundance ((adjusted) p-value below significance level) are colored in red. Regions having zero reads in the control group, i.e. no reads from any replicate for the control condition reside within the region, are not shown in this plot: DESeq assigns these regions a non-numeric base mean resp. fold-change ('NA' or 'Inf'), edgeR computes very large or small fold-changes resp. very small log-concentration values for these regions, which might cause problems with the scale of the graph. If you specified thresholds for the fold-change, they are shown in the plot as blue dashed lines.
94.fold_change_all_plot.png (only if test method is 'edgeR')
The same plot as 93.fold_change_plot.png, but containing all transcripts, also those with a total count of zero in (at least) one group
95.volcano_plot.png
graphics in PNG-format with a volcano plot of p-value or adjusted p-value (depending on the parameter setting for p-value adjustment) (y-axis) against fold-change in the read abundance of treatment versus control (x-axis). For edgeR, the measure of read abundance is the concentration, for DESeq, it is the base mean, i.e. the values shown in the plot can be found in 92.output_replicate_analysis columns 2 or 3, and 5, or in 10.region_summary.tsv columns 9 or 10, and 11. Each data point resembles a region found by merging/intersecting clusters from the treatment group. Regions having zero reads in the control group, i.e. no reads from any replicate for the control condition reside within the region, are not shown in this plot: DESeq assigns these regions a non-numeric fold-change ('NA' or 'Inf'), edgeR computes very large or small fold-change values for these regions, which might cause problems with the scale of the graph. If you specified thresholds for the fold-change, they are shown in the plot as blue dashed lines. To avoid problems with the plotting routine, p-values < 1e-311 are omitted from the volcano plot.
96.volcano_all_plot.png (only if test method is 'edgeR')
The same plot as 95.volcano_plot.png, but containing all transcripts, also those with a total count of zero in (at least) one group. To avoid problems with the plotting routine, p-values < 1e-311 are omitted from the volcano plot.

Cluster Classification

The cluster classification provides information for

the correlation of the found clusters with genomic elements, and the enrichment of the clusters for the genomic element (example: an enrichment value of 2 for promoters means, that the clusters are found twice as often overlapping promoters as expected from the fraction of promoters in the genome)
the distribution of the clusters on the individual chromosomes (can be displayed or hidden by clicking on the "Show/Hide Details" link)

Cluster Classification: Details for each Cluster

Download cluster classification

The file that can be downloaded here contains the classification for each cluster:

The file contains 7 columns (tab-separated) for each region:

  1 : read id
  2 : contig/chromosome accession number
  3 : chromosome
  4 : strand
  5 : start position of the read
  6 : end position of the read
  7 : genomic elements the read is associated with
      intergenic (intergenic region)
      exon
      intron
      partial (overlapping with exon)
      promoter
                An individual read is assigned to one of the four classes
                intergenic, exon, intron, partial and can be assigned to
                the class promoter in addition.

1428	NC_000001	chr1	0	1348237	1348455	intergenic
1429	NC_000001	chr1	0	2311642	2311953	promoter intron
1430	NC_000001	chr1	0	2450272	2450768	exon
1431	NC_000001	chr1	0	2469265	2469512	intergenic
1432	NC_000001	chr1	0	3556424	3556796	promoter partial
1433	NC_000001	chr1	0	3614623	3614831	exon

Sequence Extraction

If this step was selected in the parameter section, the DNA sequences for all clusters are extracted and written to a file in FASTA format with Genomatix annotation.

The number of sequences and the total number of basepairs is given, and the first few lines of the file are displayed. The complete file can be downloaded or saved to the project management for further analysis with other tasks.

TFBS Overrepresentation

When this option is checked, all cluster sequences are automatically checked for transcription factor binding sites (TFBS) (using MatInspector). The TFBS occurence is compared to expected values based on the background of occurrences of the TFBS

within the whole genomic sequence of the species
within all promoters as annotated by Genomatix for this species

The most overrepresented TFBS compared to genomic background as well as for promoter background is listed in the output. Additionally a link to the complete output table containing occurence, overrepresentation values and Z-scores for all analysed TBFSs is supplied ( table description) .

Definition of new TFBS

If this option is checked, a subset of cluster sequences is selected and then the Genomatix program CoreSearch is used to define new common motifs from those sequences.
The program automatically sorts all cluster sequences by quality, and the user-defined number of "best sequences" is used for CoreSearch. Here, the definition of "best" cluster sequences depends on the selected clustering algorithm:

for MACS the clusters are sorted by decreasing MACS score
for NGSAnalyzer the clusters are sorted by decreasing NE-value
if a Audic-Claverie evaluation was used, the clusters are sorted by increasing p-values

The motif(s) defined by CoreSearch are listed in the output together with their IUPAC representation and the re-value of the corresponding matrix. A link to the detailed CoreSearch output is given. From this page, the motif can be saved as a matrix for further downstream analysis, e.g. to search other sequences (e.g. all cluster sequences) for matches to this motif.

The sequences actually used for CoreSearch can be also downloaded from the result or saved directly to the Genomatix project management, to be used with other Genomatix tasks.

3. Download of Data Files

All result files that can be downloaded separately from the result page together with the statistics files (in text format) can be downloaded as an archive (tar-file).

RegionMiner subtask: Complete Workflow for ChIP-Seq Analysis (only available on GGA)

Introduction

Differential Analysis

Parameters

Output

1. Analysis Parameters

2. Main Results, depending on user selection

Sample Read Classification and Statistics

Peak Finding / Cluster Generation

Cluster Classification

Cluster Classification: Details for each Cluster

Sequence Extraction

TFBS Overrepresentation

Definition of new TFBS

3. Download of Data Files

RegionMiner subtask: Complete Workflow for ChIP-Seq Analysis
(only available on GGA)