View source for ChIP-Seq Tutorial

From GenPlay, Einstein Genome Analyzer

You do not have permission to edit this page, for the following reason:
The action you have requested is limited to users in the group: contributor.
You can view and copy the source of this page.
'''Goal:''' The objective of this tutorial is to illustrate how GenPlay can be used to isolate peaks from the data generated from a ChIP-Seq experiment. Then, to generate a list of genes that have a peak in their promoter and finally to associate  the score of the peak summit with each promoter.

'''Prerequisite:''' GenPlay need to be installed on your computer. If you haven't installed GenPlay yet, please visit the [[Downloads]] page and follow the instruction to download and install GenPlay.

'''Note:''' The final result of this tutorial is available as a project that can be loaded from the [[Projects#ChIP-Seq Tutorial| Projects]] page of this website.

== Start a New Project ==
After starting GenPlay you will be prompted to select a name, a clade, a genome and an assembly for your project.
You can enter "ChIP-Seq Tutorial" for the name, select the mammal clade and human genome (figure 1).
[[image: ChIP-Seq_tutorial_project_frame.png|center|frame|Figure 1: New Project Window]]
Then, click on the tool box button on the assembly line. A new window will appear allowing you to select chromosomes. For this tutorial we will work only on the first chromosome so loading times are shorter and the amount of memory needed is smaller. You can select all the other chromosomes and then click on the unselect button (figure 2).
[[image: ChIP-Seq_tutorial_project_chromosomes.png|center|frame|Figure 2: Project Chromosomes]]

== Load the Files ==
We first need to download the [http://www.genplay.net/library/tutorials/ChIP-Seq/ChIP-Seq_tutorial_IP.zip ChIP-Seq IP file],  the [http://www.genplay.net/library/tutorials/ChIP-Seq/ChIP-Seq_tutorial_Input.zip ChIP-Seq input control sample file] and the  [http://www.genplay.net/library/tutorials/ChIP-Seq/Genes_RefSeq_hg19_09.20.2013.zip RefSeq gene annotation file] and to uncompress this files.

The input file should contain reads obtained by sequencing the sample before the immnuno-precipitation step. This control is necessary to normalize for sequence dependent variations in the efficiency of the library preparation and sequencing procedure. One major difficulty is that to be able to really normalize for these parameters a considerable depth of sequencing of the input is necessary since this sample is not enriched.  If the depth is too low to be able to completely normalize, it is preferable to subtract the input from the IP rather than divide (which tend to magnify the differences).
The IP (immuno-precipitated) sample should contain the result of the sequencing of the sample after the precipitation step.

We now need to load the IP data. To achieve that you will need to right click on the track handler of the 1st row. A popup menu should appear, select the ''Add Layer(s)'' option (Figure 3).
[[image:ChIP-Seq_tutorial_add_layer.png|center|thumb|800px|Figure 3: Add Layer(s) Menu]]

Then select the IP file that you downloaded earlier. You will be prompted to choose a layer type. Choose ''Variable or Fixed Window Layer''.

We now need to tell GenPlay how to load the data. You can keep the default name for the track or change it if you prefer. Then, you need to check ''Bin Data'' and choose a window size and a method of score caculation. The size of the window that you should choose depends on the number of reads that are available. The smaller the windows the higher the resolution. For this example, we will choose a window of 100 bp. The option for the score calculation are discussed in the Documentation.  For the type of files used in this example, you should choose ''Addition'' as the method for the score calculation. 

To obtain the narrowest peak possible, it is generally advisable to correct the strand bias that is caused by the fact that the cross-linked DNA fragments are sequenced from the end while the actual binding site might be anywhere within the immuno-precipitated fragments. (reviewed in Wilbanks EG et al.  Evaluation of algorithm performance in ChIP-seq peak detection. PLoS One. 2010 Jul 8;5(7):e11471.) 

To correct the strand bias we need to specify what is the average fragment length of our library and what is the length of the sequenced reads. In our case we know that the average fragment length is 300bp and the reads were 50bp. Check the ''Define Fragment Length'' box and enter these informations in the ''Fragment Length'' section of the window. At the end you should have something similar to figure 4.
[[image:ChIP-Seq_tutorial_new_layer_window.png|center|frame|Figure 4: New Layer Window]]

Click on Okay and wait for the track to load. That can take more than a minute.
Once the data is loaded, double click on the handle of the track 1 to adjust the display options. Check ''Show horizontal lines', uncheck ''Auto rescale'' and enter 100 for the ''Maximum Score'' (figure 5).
[[image:ChIP-Seq_tutorial_track_options.png|center|thumb|800px|Figure 5: Track Options]]

You can now look at your data by zooming using the mouse wheel and moving within the chromosome by pressing the left button of the mouse on the track and releasing it after moving the mouse on the left or on the right.
[[image:ChIP-Seq_tutorial_IP_layer_result.png|center|thumb|800px|Figure 6: IP Layer Loaded]]

<br style="clear: both" />

== Load Input ==
It's now time to load the input control sample track (the track resulting from sequencing the cross-linked chromatin before the immuno-precipitation). This track will help us remove the peaks (enrichment region) that are non specific, either because of problem with the method (cross-linking,antibody etc...) or because of sequencing artifact. These artifacts can be caused either by preferential sequencing of specific fragment by the instruments or by differences between the genome sequenced and the genome assembly used to align the reads, since any repeat in the sequenced genome that is present as a none
repeated region in the genome assembly will result in a peak. 

Identifying these outliers can be quite difficult. One useful method is to compare the IP  libraries with a control (input) library for the same sample. It's what we will to do in this tutorial. 

To load the input file, right click on the handler of an empty track, choose the ''Add Layer(s)'' option. Then select the Input file downloaded earlier and choose ''Variable or Fixed Window Layer''. On the next window set the parameters as shown on the figure 7.
[[image:ChIP-Seq_tutorial_input_layer_window.png|center|frame|Figure 7: Load Control Parameters]]

Once the layer is loaded you can modify the display options of the track by double clicking on the track handler the same way you did for the IP layer.

<br style="clear: both" />

== Normalize IP and Input ==
In order to be comparable, the IP and the input tracks need first to be normalized.

The normalize operation of GenPlay divides each score by the sum of all the scores of the track and multiply the result by a large constant specified by the user. The only purpose of the constant is to make the score more readable.

Let's start by normalizing the IP.  You need first to right click on the track handler. After that, select the ''Normalize'' option of the layer sub-menu as shown on figure 8. You can keep the default constant, just click OK (figure 9). After that operation the score are expressed as number of read per windows per 1 billion reads (if you kept the default constant).

Once it's done repeat the operation with the input track.

[[image:ChIP-Seq_tutorial_normalize1.png|center|thumb|800px|Figure 8: Normalize Menu]]
[[image:ChIP-Seq_tutorial_normalize2.png|center|frame|Figure 9: Normalize Dialog Box]]

<br style="clear: both" />

== Add Constant to Input ==
One of the specificities of GenPlay is that while computing an operation, the windows with a score of zero are not taken into account.  

This is helpful in most cases because we don't know if a zero window is caused by the presence of a repeat (which cannot be map), a region that is not sequenced in the genome assembly or if it is a mappable  window in which no read mapped. 

In the case of this tutorial though, since we want to compare the IP track to the input track, we want to be able to do this comparison even when a window of the input has a null score. Since the comparison will be a subtraction we can add a really small constant (so no window has a score of zero) without impacting significantly the result of the subtraction.

We chose to add a constant of 0.01 to all the windows of the input track. In order to do that, you first need to right click on the handler of the input track (figure 10). Then, choose the "Constant Operation" option of the lalyer sub-menu. Enter 0.01 in the next windows as shown on figure 11 and check the ''Apply to null windows'' option before clicking OK.

[[image:ChIP-Seq_tutorial_add_cst1.png|center|frame|Figure 10: Addition Constant Menu]]
[[image:ChIP-Seq_tutorial_add_cst2.png|center|frame|Figure 11: Constant Dialog Box]]

<br style="clear: both" />

== Compare Input to Input ==
We can now compare the IP to the input track in order to filter out the artifactual peaks. As we already decided, we want to subtract the input scores from the IP scores.

To do so, you need to right right click on the handler of the IP track and select the IP layer sub-menu. Then click on "Two Tracks Operation" as shown on figure 12.

In the first dialog window you need to choose the second track for the operation. Select the input track (figure 13). 

The operation to choose is subtraction (figure 14). You also need to choose in which empty track you want the result to be displayed. The figure 15 shows the result of the operation.

We choose to subtract rather than divide (as would be done in a PCR experiment) because in chip-Seq experiments the input sample is not enriched wile the IP might be enriched 10 to 20 time or more. As a result, to have comparable coverage the input would have to be sequenced 10 to 20 times deeper than the IP sample (which would be very expensive). This lower coverage for the input control causes a sampling error that can drastically affect the results if the data are normalised by calculating the ratio IP/IN.  The sampling error problem on the input  increase when the binning size decrease (when the window size used loading the track decrease). An alternative to subtracting the input from the IP track is to smooth the input track in order to attenuate the sampling error on that track. That can be done in GenPlay using either the moving average or the Gaussian smoothing operation. In practice to decrease the sampling error, we find that smoothing  over 5 to 10 kb is necessary. Finally, whatever the method use, some of the peak present in the IP track might not have been completely eliminated from the IP track, or some artifact might be present only in the IP track. If that is the case, one can also eliminate the outliers by simply using the filter to eliminate the highest peak. This can be done easily in GenPlay using either either absolute or relative values in the filter operation. We find that there are no optimal method to analyse all tracks and that it is best for the users to proceed by trial and error and exercise their judgment in performing these operations.
[[image:ChIP-Seq_tutorial_subtract1.png|center|thumb|800px|Figure 12: Two Tracks Operation Menu]]
[[image:ChIP-Seq_tutorial_subtract2.png|center|frame|Figure 13: Dialog to Choose the Second Track]]
[[image:ChIP-Seq_tutorial_subtract3.png|center|frame|Figure 14: Dialog to Choose the Operation Type]]
[[image:ChIP-Seq_tutorial_subtract4.png|center|thumb|800px|Figure 15: Result of the Subtraction]]

We can now change the display properties of the track by double clicking on the track handler of the layer resulting from the subtraction like we did for the two previous tracks. From the same window we can also change the layer name as shown in figure 16.
[[image:ChIP-Seq_tutorial_layer_name.png|center|thumb|800px|Figure 16: Change Layer Name]]

<br style="clear: both" />

== Filter Negative Values ==
Since we're just interesting in the enriched regions we can try to remove all the negative values in order to make the track easier to comprehend.

To filter the track you need to select the "Filter" option of the layer sub-menu (figure 17). On the left panel of the filter window select the threshold option. Set the parameters of the right panel as shown on figure 18. The figure 19 shows the result of the filter operation.

[[image:ChIP-Seq_tutorial_filter1.png|center|thumb|800px|Figure 17: Filter Operation Menu]]
[[image:ChIP-Seq_tutorial_filter2.png|center|frame|Figure 18: Filter Operation Window]]
[[image:ChIP-Seq_tutorial_filter3.png|center|thumb|800px|Figure 19: Filter Operation Result (the negative values of the track 5 have been removed)]]

We now need to remove the background noise and to keep only the islands (the peaks).

<br style="clear: both" />

== Isolate Peaks ==
The goal of this step is to remove the background noise from the track so just the peaks remain.

To do so, right click on the track handler, choose the "Input - Output" sub-menu and click on the "Find Peaks" option (figure 20). 
[[image:ChIP-Seq_tutorial_find_peaks1.png|center|thumb|800px|Figure 20: Find Peaks Operation]]

After the ''Find Peaks'' dialog opens, choose the "Island Finder" option on the left panel and set the parameters as shown on the figure 21.
[[image:ChIP-Seq_tutorial_find_peaks2.png|center|frame|Figure 21: Find Peaks Menu]]

The island finder is described in the documentation section of this website.
You'll notice that the selected output is "Peak Summits". This means that for each island, the only selected window will be the one with greatest score on the island. This generally corresponds to the window with the highest probability of containing the binding site. Other software can be used to identify the exact sequence of the binding site. In addition to choosing the summit as the output, one could add the Island score or the  "start value". These two values are a much more quantitative estimate of the level of enrichment of the signal than the summit value.  
The result should be similar to what is shown on figure 22.
[[image:ChIP-Seq_tutorial_find_peaks3.png|center|thumb|800px|Figure 22: Find Peaks Result (the track 4 contains only the summit of the peaks)]]

<br style="clear: both" />

== Extract Gene Promoters ==
First, we need to load a gene track. Right click on the track handler of the track with the summit and select "Add Layer(s)". Select the RefSeq file that we've already downloaded when prompted and then select ''Gene Annotation Layer''.

When it's done, right click on the track handler of the gene track and select "Extract Intervals" in the RefSeq sub-menu (Figure 23).
[[image:ChIP-Seq_tutorial_extract_promoters1.png|center|thumb|800px|Figure 23: Extract Intervals Menu]]

A dialog box will pop-up. We decide to define a promoter as a region that starts 500 bp before a gene start position and ends 500bp after. In order to do so, fill in the parameters as shown in figure 24.
[[image:ChIP-Seq_tutorial_extract_promoters2.png|center|frame|Figure 24: Extract Intervals Dialog]]

You'll finally be asked to select the result track position in the track list. The result track represents only the promoters of the genes of the reference track (figure 25).
[[image:ChIP-Seq_tutorial_extract_promoters3.png|center|thumb|800px|Figure 25: Gene Promoters (the track 4 contains only the gene promoters)]]

<br style="clear: both" />

== Score Promoters ==
Now that we have a track with the peaks and a track with the promoters we can score the promoters using the score of the peaks and export the result as a bed file.

To score the promoters, right click on the handler of the track with the promoters and select the "Score Exons" option of the gene layer sub-menu (figure 26).
[[image:ChIP-Seq_tutorial_score_exons1.png|center|thumb|800px|Figure 26: Score Exons Menu]]

You'll be prompted to choose the layer containing the scores. Select the layer with the peaks extracted (figure27).
[[image:ChIP-Seq_tutorial_score_exons2.png|center|frame|Figure 27: Choose Layer with Scores]]

Then select ''Maximum Coverage'' method of calculation (figure 28). Note that the color of the promoters represents the scores associated to the promoters (as described in the gene layer section of the documentation). The result is shown in figure 29.
[[image:ChIP-Seq_tutorial_score_exons3.png|center|frame|Figure 28: Method of Calculation]]
[[image:ChIP-Seq_tutorial_score_exons4.png|center|thumb|800px|Figure 29: Score Exons Result]]

<br style="clear: both" />

== Export Result ==
The last thing we need to do is to export the result of our analysis. Right click on the track handler of the track with the gene layer. Select "Save As" in the gene layer submenu (figure 30). Choose where you want to save the track and make sure that the file type is set to Bed file. 
[[image:ChIP-Seq_tutorial_save_layer.png|center|thumb|800px|Figure 30: Save Layer]]

Once it's done, you can open the file that you created with a text editor such as notepad. You'll notice that the result file contains the position  (field 1 to 3) of the promoters, the name of the genes (field 4), the strand of the gene (field 6) as well as the scores of the promoters (field 5). For more details about the result file you can refer to the File Type section of the documentation.
Return to ChIP-Seq Tutorial.
Retrieved from "http://genplay.net/wiki/index.php/ChIP-Seq_Tutorial"