LBMC » NucleoMiner

Tutorial

========

This tutorial describes steps allowing to perform quantitave analysis of 

nucleosomal epigenome. We assume that files are organised around a given 

hierarchie and that all command lines are launched from project's root.

This tutorial is divided into t=wo main parts. First one consists in the python 

script `wf.py` that aligns and convert Illumina reads. Second one is the R 

script `main.r` that extracts information (nucleosome position and indicators) 

from the dataset.

Python and R Common Configuration File

--------------------------------------

First of all we define in one place some configuration variables that will be launch by python and R scripts. This file is **configurator.py**. The execution of this python script dumps variables into the **nucleo_miner_config.json** that will be launch by both kind of scriopts (R and puython).

To do this launch at the root of your project the following command line:

.. code:: bash

  python src/current/configurator.py

$$$ other python script to describe:

- libcoverage.py

- wf.py

Dataset and Configuration Variables

-----------------------------------

We want to compare nucleosomes of 3 yeast strains: 

- BY

- RM

- YJM

For each strain we perform Mnase-Seq and ChIP-Seq using the 5 following 

markers: 

- H3K4me1

- H3K4me3

- H3K9ac

- H3K14ac

- H4K12ac

In order to simplify the design of experiment, we considere Mnase as a marker. 

For each couple `(strain, marker)` we perform 3 replicates. So, theoritically 

we should have `3 * (1 + 5) * 3 = 54` samples. In practice we only obtain 2 

replicates for `(YJM, H3K4me1)`. Each one of the 53 samples is indentify by a 

uniq identifier. The file `CSV_SAMPLE_FILE` sums up this information.

.. autodata:: configurator.CSV_SAMPLE_FILE

    :noindex: 

We use a convention to link sample and Illumina fastq outputs. Illumina output 

files of the sample `ID` will be stored in the directory 

`ILLUMINA_OUTPUTFILE_PREFIX` + `ID`. For example, sample 41 outputs will be 

stored in the directory `data/2012-09-05/FASTQ/Sample_Yvert_Bq41/`.

.. autodata:: configurator.ILLUMINA_OUTPUTFILE_PREFIX

    :noindex: 

For BY (resp. RM and YJM) we use following reference genome 

`saccharomyces_cerevisiae_BY_S288c_chromosomes.fasta`

(resp. `saccharomyces_cerevisiae_rm11-1a_1_supercontigs.fasta` and 

`saccharomyces_cerevisiae_YJM_789_screencontig.fasta`).

The index `FASTA_REFERENCE_GENOME_FILES` stores this information.

.. autodata:: configurator.FASTA_REFERENCE_GENOME_FILES

    :noindex: 

Each chromosome/contig is identify in the fasta file by an obscure identifier. 

For example, BY chromosome I is identify by `gi|144228165|ref|NC_001133.7|` when 

TemplateFilter is waiting for an integer. So, we translate it. The index 

`FASTA_INDEXES` stores this translation.

.. autodata:: configurator.FASTA_INDEXES

    :noindex: 

From a pragamatical point of view we discard some part of the genome (repeated 

sequence etc...). The list of the black listed area is explicitely detailled in 

`AREA_BLACK_LIST`.

.. autodata:: configurator.AREA_BLACK_LIST

    :noindex: 

For BY-RM (resp. BY-YJM and RM-YJM) genome sequence alignment we use previously 

compute .c2c file `data/2012-03_primarydata/BY_RM_gxcomp.c2c` (resp. 

`BY_YJM_GComp_All.c2c` and `RM_YJM_gxcomp.c2c`). For more information about 

.c2c files, please read section 5 of the manual of `NucleoMiner`, the old 

version of `NucleoMiner2` 

(http://www.ens-lyon.fr/LBMC/gisv/NucleoMiner_Manual/manual.pdf).

.. autodata:: configurator.C2C_FILES

    :noindex: 

`nucleominer` uses specific directory to work in, these are described in 

`INDEX_DIR`, `ALIGN_DIR` and `LOG_DIR`.

Finally, `nucleominer` use external ressources, the path to these resspources 

are describe in `BOWTIE_BUILD_BIN`, `BOWTIE2_BIN`, `SAMTOOLS_BIN`, 

`BEDTOOLS_BIN` and `TF_BIN` and `TF_TEMPLATES_FILE`.

All paths, prefixes and indexes could be change in the 

`src/current/nucleominer_config.json` file.

.. autodata:: wf.json_conf_file

    :noindex: 

Preprocessing Illumina Fastq Reads for Each Sample

--------------------------------------------------

This preprocessing step consists in the 4 main steps embed in the `wf.py` and 

described bellow. As a preamble, this script computes `samples` `samples_mnase` 

and `strains` that will be used along the 4 steps.

.. autodata:: wf.samples

    :noindex: 

.. autodata:: wf.samples_mnase

    :noindex: 

.. autodata:: wf.strains

    :noindex: 

Creating Bowtie Index from each Reference Genome

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

For each strain, we need to create bowtie index. Bowtie index of a strain is a 

tree view of the genemoe reference for this strain. It will be used by 

bowtie to align reads. This step is performed by the following part of the 

`wf.py` script:

.. literalinclude:: ../../../snep/src/current/wf.py

   :start-after: # _STARTOF_ step_1

   :end-before: # _ENDOF_ step_1

   :language: python

The following table sum up involved file sizes and process durations concerning 

this step.

======  ======================  ======================  ================

strain  fasta genome file size  bowtie index file size  process duration

======  ======================  ======================  ================

BY      12 Mo                          25 Mo                    11 s.

RM      12 Mo                          24 Mo                    9 s.

YJM     12 Mo                          25 Mo                    11 s.

======  ======================  ======================  ================

Aligning Reads to Reference Genome 

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Next, we launch bowtie to align reads to the reference genome. It produces a 

`.sam` file that we convert into a `.bed` file. Binaries for `bowtie`, `samtools` 

and `bedtools` are wrapped using python `subprocess` class. This step is 

performed by the followinw part of the `wf.py` script:

.. literalinclude:: ../../../snep/src/current/wf.py

   :start-after: # _STARTOF_ step_2

   :end-before: # _ENDOF_ step_2

   :language: python

Convert Aligned Reads for TemplateFilter

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

TemplateFilter use particular input format for reads, so we convert `.bed` file. 

TemplateFilter expect reads as following: `chr coord strand #read` where:

- chr is the number of the chromosome;

- coord is the coordinate of the reads;

- strand is `F` for forward and `R` for reverse;

- #reads the number of reads for this position.

Each entry is *tab*-separated.

**WARNING** for reverse strand bowtie returns the position of left first 

nucleotid when TemplateFilter is waiting for right one. So this step takes it 

into account and add lenght of reads (in our case 50) to reverse reads 

coordinate.

This step is performed by the followinw part of the `wf.py` script:

.. literalinclude:: ../../../snep/src/current/wf.py

   :start-after: # _STARTOF_ step_3

   :end-before: # _ENDOF_ step_3

   :language: python

The following table sum up number of reads, involved file sizes and process durations concerning 

the two last steps. In our case, aligment process have been multuthreaded over over 3 cores.

==  ==============  =========================  ======  ================  ==================  ================  

id  Illumina reads  aligned and filtred reads  ratio   `.bed` file size  TF input file size  process duration

==  ==============  =========================  ======  ================  ==================  ================

1   16436138        10199695                   62,06%  1064 Mo           60  Mo              383   s.

2   16911132        12512727                   73,99%  1298 Mo           64  Mo              437   s.

3   15946902        12340426                   77,38%  1280 Mo           65  Mo              423   s.

4   13765584        10381903                   75,42%  931  Mo           59  Mo              352   s.

5   15168268        11502855                   75,83%  1031 Mo           64  Mo              386   s.

6   18850820        14024905                   74,40%  1254 Mo           69  Mo              482   s.

7   15591124        12126623                   77,78%  1163 Mo           72  Mo              405   s.

8   15659905        12475664                   79,67%  1194 Mo           71  Mo              416   s.

9   14668641        10960565                   74,72%  1052 Mo           70  Mo              375   s.

10  14339179        10454451                   72,91%  1049 Mo           51  Mo              363   s.

11  18019895        13688774                   75,96%  1378 Mo           59  Mo              474   s.

12  13746796        10810022                   78,64%  1084 Mo           54  Mo              360   s.

13  15205065        11766016                   77,38%  990  Mo           54  Mo              381   s.

14  17803097        13838883                   77,73%  1154 Mo           60  Mo              452   s.

15  15434564        12307878                   79,74%  1032 Mo           57  Mo              394   s.

16  16802587        12725665                   75,74%  1221 Mo           48  Mo              438   s.

17  16058417        12513734                   77,93%  1192 Mo           63  Mo              422   s.

18  16154482        13204331                   81,74%  1277 Mo           52  Mo              430   s.

19  21013924        17102120                   81,38%  1646 Mo           59  Mo              555   s.

20  17213114        14433357                   83,85%  1389 Mo           53  Mo              459   s.

21  17360907        14733001                   84,86%  1203 Mo           55  Mo              450   s.

22  18136816        15389581                   84,85%  1257 Mo           53  Mo              469   s.

23  14763678        12173025                   82,45%  1140 Mo           56  Mo              393   s.

24  15541709        12890345                   82,94%  1057 Mo           48  Mo              398   s.

25  16433215        13094314                   79,68%  1241 Mo           57  Mo              433   s.

26  17370850        14264136                   82,12%  1347 Mo           51  Mo              466   s.

27  14613512        8654495                    59,22%  887  Mo           56  Mo              339   s.

28  15248545        11367589                   74,55%  1166 Mo           67  Mo              405   s.

29  14316809        10767926                   75,21%  1103 Mo           63  Mo              379   s.

30  15178058        12265794                   80,81%  1030 Mo           66  Mo              390   s.

31  14968579        11876186                   79,34%  1009 Mo           63  Mo              387   s.

32  16912705        13550508                   80,12%  1143 Mo           70  Mo              442   s.

33  16782154        12755111                   76,00%  1227 Mo           65  Mo              438   s.

34  16741443        13168071                   78,66%  1260 Mo           71  Mo              442   s.

35  13096171        10367041                   79,16%  992  Mo           62  Mo              350   s.

36  17715118        14092985                   79,55%  1404 Mo           68  Mo              483   s.

37  17288466        7402082                    42,82%  741  Mo           48  Mo              339   s.

38  16116394        13178457                   81,77%  1101 Mo           63  Mo              420   s.

39  14241106        10537228                   73,99%  880  Mo           57  Mo              348   s.

40  13784738        10598464                   76,89%  1005 Mo           64  Mo              358   s.

41  12438007        9620975                    77,35%  911  Mo           60  Mo              326   s.

42  13853959        11031238                   79,63%  1045 Mo           64  Mo              365   s.

43  12036162        6654780                    55,29%  684  Mo           46  Mo              268   s.

44  13873129        10251074                   73,89%  1048 Mo           61  Mo              365   s.

45  19817751        14904502                   75,21%  1520 Mo           72  Mo              528   s.

46  13368959        10818619                   80,92%  912  Mo           63  Mo              350   s.

47  7566467         6139001                    81,13%  520  Mo           44  Mo              201   s.

48  32586928        21191363                   65,03%  1816 Mo           82  Mo              766   s.

49  30733184        18791373                   61,14%  1801 Mo           89  Mo              721   s.

50  41287616        30383875                   73,59%  2911 Mo           112 Mo              1065  s.

51  40439965        31177914                   77,10%  2981 Mo           117 Mo              1070  s.

53  40876476        33780065                   82,64%  3316 Mo           103 Mo              1165  s.

55  52424414        47117107                   89,88%  3811 Mo           119 Mo              1477  s.

==  ==============  =========================  ======  ================  ==================  ================  

For some reasons (manipulation efficiency, e.g. PCR...), we remove samples 33, 45, 48 and 55.

Run TemplateFilter on Mnase Samples

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Finally, for each sample we perfome TemplateFilter analysis. 

**WARNING** TemplateFilter returns a list of nucleosomes. Each nucleosome is 

define by its center and its width. An odd width leads us to considere non 

interger lower and upper bound.

**WARNING** TemplateFilter is not design to deal with replicate. So we choose to 

keep a maximum of nucleosome and filter in a second time using the benefit of 

replicate. To do that we set a low correlation threshold parameter (`0.5`) and a 

particularly high value of overlaping (`300%`).

This step is performed by the followinw part of the `wf.py` script:

.. literalinclude:: ../../../snep/src/current/wf.py

   :start-after: # _STARTOF_ step_4

   :end-before: # _ENDOF_ step_4

   :language: python

==  ======  ==========  =============  ================

id  strain  found nucs  nuc file size  process duration

==  ======  ==========  =============  ================

1    BY     96214       68 Mo          1022 s.                     

2    BY     91694       65 Mo          1038 s.                      

3    BY     91205       65 Mo          1036 s.                       

4    RM     88076       62 Mo          984 s.                      

5    RM     90141       64 Mo          967 s.                      

6    RM     87517       62 Mo          980 s.                      

7    YJM    88945       64 Mo          566 s.                      

8    YJM    88689       64 Mo          570 s.                      

9    YJM    88128       63 Mo          565 s.                    

==  ======  ==========  =============  ================

Inferring Nucleosome Position and Extracting Read Counts

--------------------------------------------------------

The second part of the tutorial uses `R` (http://http://www.r-project.org). It consists in a set of R scripts taht will be sourced in an R console launched at the root of your project. the R srcipts are:

  - headers.R

  - extract_maps.R

  - compare_common_wp.R

  - split_samples.R

  - count_reads.R

  - get_size_factors  

  - launch_deseq.R

The Script headers.R

^^^^^^^^^^^^^^^^^^^^

The script header.R is included in each other scripts. It is in charge of: 

  - launching libraries used in thes scripts

  - launching configuration (design, strain, marker...)

  - computing and caching CURs

In your R console, run the following command line:

.. code:: bash

  R CMD BATCH src/current/header.R

The Script extract_maps.R

^^^^^^^^^^^^^^^^^^^^^^^^^

This script is in charge of extracting Maps for well positioned and fuzzy nucleosomes. First of all, this script computed intra and inter strain nucleosome maps for each CUR. This step is executed in parallel on many cores using the BoT library. Next, it collects results and produces well positioned, fuzzy and UNR maps.

The well-positioned map for BY is collected in the result directory and is called **BY_wp.tab**. It is composed of following columns:

 - chr, the number of the chromosome 

 - lower_bound, the lower bound of the nucleosome

 - upper_bound, the upper bound of the nucleosome 

 - cur_index, index of the CUR

 - index_nuc, the index of the nucleosome in the CUR

 - wp, 1 if it is a well positioned nucleosome, 0 else

 - nb_reads, the number of reads that supports this nucleosome

 - nb_nucs, the number of TemplateFilter nucleosome across replicates (= the number of replicates if it is a well-positioned nucleosome)

 - llr_1, for a well-positioned nucleosome, it is the LLR1 between the first and the second TemplateFilter nucleosome.

 - llr_2, for a well-positioned nucleosome, it is the LLR1 between the second and the first TemplateFilter nucleosome.

 - wp_llr, for a well-positioned nucleosome, it is the LLR2 overall TemplateFilter nucleosomes.

 - wp_pval, for a well-positioned nucleosome, it is the p-value chi square test obtained with the LLR2 (**1-pchisq(2.LLR2, df=4)**)

 - dyad_shift, for a well-positioned nucleosome, it is shift between the two extreme TemplateFilter nucleosome dyad positions. 

The fuzzy map for BY is collected in the result directory and is called **BY_fuzzy.tab**. It is composed of following columns:

 - chr, the number of the chromosome 

 - lower_bound, the lower bound of the nucleosome

 - upper_bound, the upper bound of the nucleosome 

 - cur_index, index of the CUR

The common well-position map for BY and RM strains is collected in the result directory and is called **BY_RM_common_wp.tab**. It is composed of following columns:

 - cur_index, the index of the CUR

 - index_nuc_BY, the index of the BY nucleosome in the CUR

 - index_nuc_RM,the index of the RM nucleosome in the CUR

 - llr_score, the LLR3 score between th eBy and RM nucleosomes

 - common_wp_pval,  the p-value chi square test obtained with the LLR3 (**1-pchisq(2.LLR3, df=2)**)

The common UNR map for BY and RM strains is collected in the result directory and is called **BY_RM_common_unr.tab**. It is composed of following columns:

 - cur_index, the index of the CUR

 - index_nuc_BY, the index of the BY nucleosome in the CUR

 - index_nuc_RM,the index of the RM nucleosome in the CUR

To execute this script, run the following command line in your R console:

.. code:: bash

  source("src/current/extract_maps.R")

The Script compare_common_wp.R

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

This script is used to compare inter strain distances between common well-positioned nucleosomes. 

For example, it compute the file **BY_RM_common_wp_diff.tab** that contains dyad shifts between two well-positioned nucleosomes. It is composed of following columns:

 - cur_index, the index of the CUR

 - index_nuc_BY, the index of the BY nucleosome in the CUR

 - index_nuc_RM,the index of the RM nucleosome in the CUR

 - llr_score, the LLR3 score between th eBy and RM nucleosomes

 - common_wp_pval,  the p-value chi square test obtained with the LLR3 (**1-pchisq(2.LLR3, df=2)**)

 - diff, the dyad shifts between two well-positioned nucleosomes

It also translates well-positioned nucleosome maps from a strain to an other strain and stores it into a table. 

For example, the file **results/2014-04/RM_wp_tr_2_BY.tab** contains RM well-positioned nucleosome translated into the BY genome referential. It is composed of following columns:

 - strain_ref, the reference genome (in which positioned are defined)

 - begin, the translated lower bound of the nucleosome

 - end, the translated upper bound of the nucleosome

 - chr, the number of chromosome for the reference genome (in which positioned are defined)

 - length, the length of the nucleosome (could be negative)

 - cur_index, the index of the CUR

 - index_nuc, the index of the nucleosome in the CUR

To execute this script, run the following command line in your R console:

.. code:: bash

  R CMD BATCH src/current/compare_common_wp.R

Split and Compress Samples According CURs

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. code:: bash

  R CMD BATCH src/current/split_samples.R

Count Reads for Each Nucleosome

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. code:: bash

  R CMD BATCH src/current/count_reads.R

Get Size Factors Using DESeq

^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. code:: bash

  R CMD BATCH src/current/get_size_factors.R

Performing DESeq Analysis

^^^^^^^^^^^^^^^^^^^^^^^^^

.. code:: bash

  R CMD BATCH src/current/launch_deseq.R

Results

-------

Output Files Organisation

^^^^^^^^^^^^^^^^^^^^^^^^^

Previous steps produce following 45 files. 

Each filename is under the form 

.. code:: bash

  results/current/[combi]_[marker]_[form]_snep.tab 

Where combi is in {BY_RM, BY_YJM, RM_YJM} for each strain combination, marker is 

in {H3K4me1, H3K4me3, H3K9ac, H3K14ac, H4K12ac} for each post translational 

histone modification and form is in {wp, fuzzy, wpfuzzy} considering well 

positioned nucleosomes, fuzzy nucleosomes or both for SNEP computation.

chr_BY lower_bound_BY upper_bound_BY index_nuc_BY chr_RM lower_bound_RM upper_bound_RM index_nuc_RM roi_index form BY_Mnase_Seq_1 BY_Mnase_Seq_2 BY_Mnase_Seq_3 RM_Mnase_Seq_4 RM_Mnase_Seq_5 RM_Mnase_Seq_6 BY_H3K14ac_36 

BY_H3K14ac_37 BY_H3K14ac_53 RM_H3K14ac_38 RM_H3K14ac_39 pvalsGLM 

For each file, there is 1 line per nucleosome and each line is composed of many columns divided into 3 main topics:

  - nuc information

  - number opf reads for each sample

  - DESeq analysis results.

For exemple for the file *BY_RM_H3K14ac_wp_snep.tab* informations are: 

  - chr_BY, the BY chr involved

  - lower_bound_BY, the lower bound of the BY nuc

  - upper_bound_BY, the upper_bound of the BY nuc

  - index_nuc_BY, the index of the nuc in the entire list of BY nucs

  - chr_RM, lower_bound_RM, upper_bound_RM, index_nuc_RM 

	are the same information for the RM strain

  - roi_index, the index of the region of interrest involved.

Next cols concern indicators for each sample. They are labeled [strain]_[marker]_[sample_id] and each value represents the number of reads for the current nuc for the sample *sample_id*. 

The 5 final columns concern DESeq analysis:

  - manip[a_manip] strain[a_strain] manip[a_strain]:strain[a_strain], the manip (marker) effect, the strain effect and the snep effect.  

  - pvalsGLM, the pvalue resulting of the comparison of the GLM model considering or the interaction term *marker:strain* 

  - snep_index, a boolean set to TRUE if the *pvalueGLM* value is under the threshold computed with FDR function with a rate set to 0.01%. 

It also produces the file that explicts size factor for each involved sample in differents strain combination and nucleosomal region type:

TODO: include this file... /home/filleton/analyses/snepcatalog/data/2013-10-09/current/README.txt

.. code:: bash

  results/current/size_factors.tab

Number of SNEPs

^^^^^^^^^^^^^^^

Here are the number of computed for each forms.

.. code:: bash

  [1] "wp"

         #nucs H3K4me1 H3K4me3 H3K9ac H3K14ac H4K12ac

  BY-RM  30234     520     798     83    3566      26

  BY-YJM 31298     303     619    102     103     128

  RM-YJM 29863     129     340     46    3177      18

  [1] "fuzzy"

         #nucs H3K4me1 H3K4me3 H3K9ac H3K14ac H4K12ac

  BY-RM  10748     294     308    101    1681      42

  BY-YJM 10669     122     176    124      93      87

  RM-YJM 11478      54     112     41    1389      20

  [1] "wpfuzzy"

         #nucs H3K4me1 H3K4me3 H3K9ac H3K14ac H4K12ac

  BY-RM  40982     770    1136    183    5404      73

  BY-YJM 41967     439     804    214     198     199

  RM-YJM 41341     184     468     87    4687      37

TODO: 

  - Print/study intra/inter strain LODs.

  - Check the normality of sample using Shapiro–Wilk (Hypothesis for computing LODs)