LBMC » NucleoMiner

% Generated by Sphinx.

% \def\sphinxdocclass{report}

\documentclass[letterpaper,10pt,english]{article}

\usepackage[utf8]{inputenc}

\DeclareUnicodeCharacter{00A0}{\nobreakspace}

\usepackage{cmap}

\usepackage[T1]{fontenc}

\usepackage{babel}

\usepackage{times}

\usepackage[Bjarne]{fncychap}

\usepackage{longtable}

\usepackage{multirow}

\usepackage{url}

\addto\captionsenglish{\renewcommand{\figurename}{Fig. }}

\addto\captionsenglish{\renewcommand{\tablename}{Table }}

\title{NucleoMiner2 Documentation}

\date{\today}

\author{Florent CHUFFART, Jean-Baptiste VEYRIERAS, Gael YVERT}

\newcommand{\sphinxlogo}{}

\makeindex

\def\PYGZus{\char`\_}

\begin{document}

\maketitle

\tableofcontents

\textbf{}\label{index::doc}

\section{Inferring Nucleosome Position and Extracting Read Counts}

\label{tuto:inferring-nucleosome-position-and-extracting-read-counts}

The second part of the tutorial uses R (\url{http://http://www.r-project.org}{http://http://www.r-project.org}).

It consists in 5 main steps: 

\begin{itemize}

\item {} 

Seting Up Configuration

\item {} 

Computing CURs

\item {} 

Extracting Nucleosome Maps from TemplateFilter Outputs

\item {} 

Building Count Tables from TemplateFilter Inputs

\item {} 

Analysing Count Tables using DESeq

\end{itemize}

\subsection{Seting Up Configuration}

This first step is in charge of:

\begin{itemize}

\item {} 

launching libraries used in the scripts

\item {} 

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

{\footnotesize

<<echo=TRUE, results=hide>>=

library(nucleominer)

@

<<echo=TRUE, results=verbatim>>=

package.version("nucleominer")

library(rjson)

json_conf_file = "nucleominer_config.json"

config = fromJSON(paste(readLines(json_conf_file), collapse=""))

all_samples = read.csv(config$CSV_SAMPLE_FILE, sep=";", head=TRUE, stringsAsFactors=FALSE)

head(all_samples)

markers = unique(all_samples$marker)

print(markers)

strains = unique(all_samples$strain)

print(strains)

combis = combn(strains, 2, simplify=FALSE)

print(combis)

@

}

\subsection{Computing CURs}

This second step is in charge of computing and caching Common Uinterrupted Regions (CURs). Caching means storing the information in the computer's memory.

\end{itemize}

Note that you can customize the function “translate”. This function allows you to use the alignments between genomes when performing various tasks.

\begin{itemize}

\item {} 

You may want to analyze data of a single strain (e.g. treatment/control, or only few mutations). In this case, the genome is identical across all samples and you do not need to define particular CURs (CURs are chromosomes). Simply use the default translate function which is neutral.

\item {} 

If you are analyzing data from two or more strains (as NucleoMiner2 was designed for), then you need to translate coordinates of one genome into the coordinates of another one. You must do this by aligning the two genomes, which will produce a .c2c file (see Appendice ``Generate .c2c Files'').  thenuse it to produce the list of regions and customise “translate”.

\end{itemize}

In our tutorial, we are in the second case and to perform all these steps run the following command line in your R console:

{\footnotesize

<<echo=TRUE, results=hide>>=

curs = compute_curs(diff_allowed=30, min_cur_width=4000, combis, config)

@

<<echo=TRUE, results=verbatim>>=

head(curs)

# Reduce first CUR length to accelerate demo, TODO: remove it before submission.

curs [1,3] = curs[1,3] - 110000

curs [1,5] = curs[1,5] - 110000

@

}

\subsection{Extracting Nucleosome Maps from TemplateFilter Outputs}

This third step is in charge of extracting Maps for well-positioned and sensitive nucleosomes. First of all, this script computes intra and inter-strain matches of nucleosome maps for each CUR. This step can be executed in parallel on many cores using the BoT library. Next, it collects results and produces maps of  well-positioned nucleosomes, sensitive nucleosomes and Unaligned Nucleosomal Regions .

The map of well-positioned nucleosomes for BY is collected in the result directory and is called \emph{BY\_wp.tab}. It is composed of following columns:

\begin{itemize}

\item {} 

chr, the number of the chromosome

\item {} 

lower\_bound, the lower bound of the nucleosome

\item {} 

upper\_bound, the upper bound of the nucleosome

\item {} 

cur\_index, index of the CUR

\item {} 

index\_nuc, the index of the nucleosome in the CUR

\item {} 

wp, 1 if it is a well positioned nucleosome, 0 otherwise

\item {} 

nb\_reads, the number of reads that support this nucleosome

\item {} 

nb\_nucs, the number of TemplateFilter nucleosome across replicates (= the number of replicates in which it is a well-positioned nucleosome)

\item {} 

llr\_1, for a well-positioned nucleosome, it is the LLR1 (log-likelihood ratio) between the first and the second TemplateFilter nucleosome on the chain.

\item {} 

llr\_2, for a well-positioned nucleosome, it is the LLR1 between the second and the third TemplateFilter nucleosome on the chain.

\item {} 

wp\_llr, for a well-positioned nucleosome, it is the LLR2 that compares consistency of the positioning over all TemplateFilter nucleosomes.

\item {} 

wp\_pval, for a well-positioned nucleosome, it is the p-value chi square test obtained from LLR2 (\emph{1-pchisq(2.LLR2, df=4)})

\item {} 

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

\end{itemize}

The sensitive map for BY is collected in the result directory and is called \emph{BY\_fuzzy.tab}. It is composed of following columns:

\begin{itemize}

\item {} 

chr, the number of the chromosome

\item {} 

lower\_bound, the lower bound of the nucleosome

\item {} 

upper\_bound, the upper bound of the nucleosome

\item {} 

cur\_index, index of the CUR

\end{itemize}

The map of common well-positioned nucleosomes aligned between the BY and RM strains is collected in the result directory and is called \emph{BY\_RM\_common\_wp.tab}. It is composed of following columns:

\begin{itemize}

\item {} 

cur\_index, the index of the CUR

\item {} 

index\_nuc\_BY, the index of the BY nucleosome in the CUR

\item {} 

index\_nuc\_RM, the index of the RM nucleosome in the CUR

\item {} 

llr\_score, , the LLR3 score that estimates conservation between the positions in BY and RM

\item {} 

common\_wp\_pval,  the p-value chi square test obtained from LLR3 (\emph{1-pchisq(2.LLR3, df=2)})

\item {} 

diff, the dyads shift between the positions in the two strains (in bp)

\end{itemize}

The common UNR map for BY and RM strains is collected in the result directory and is called \emph{BY\_RM\_common\_unr.tab}. It is composed of the following columns:

\begin{itemize}

\item {} 

cur\_index, the index of the CUR

\item {} 

index\_nuc\_BY, the index of the BY nucleosome in the CUR

\item {} 

index\_nuc\_RM,the index of the RM nucleosome in the CUR

\end{itemize}

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

{\footnotesize

<<echo=TRUE, results=hide>>=

build_maps(strains, combis, all_samples, curs, config)

@

<<echo=TRUE, results=verbatim>>=

BY_wp = read.table("outputs/BY_wp.tab", header=TRUE)

head(BY_wp)

RM_wp = read.table("outputs/RM_wp.tab", header=TRUE)

head(RM_wp)

BY_RM_common_wp = read.table("outputs/BY_RM_common_wp.tab", header=TRUE)

head(BY_RM_common_wp)

@

}

\subsection{Building Count Tables from TemplateFilter Inputs

}

\label{tuto:the-script-count-reads-r}

To associate a number of observations (read) to each nucleosome we run the script \emph{count\_reads.R}. It produces the files \emph{BY\_RM\_H3K14ac\_wp\_and\_nbreads.tab}, \emph{BY\_RM\_H3K14ac\_unr\_and\_nbreads.tab} \emph{BY\_RM\_Mnase\_Seq\_wp\_and\_nbreads.tab} and \emph{BY\_RM\_Mnase\_Seq\_unr\_and\_nbreads.tab}

for H3K14ac common well-positioned nucleosomes, H3K14ac UNRs, Mnase common well-positioned nucleosomes and Mnase UNRs respectively.

For example, the file \emph{BY\_RM\_H3K14ac\_unr\_and\_nbreads.tab} contains counted reads for well-positioned nucleosomes with the experimental condition ChIP H3K14ac. It is composed of the following columns:

\begin{itemize}

\item {} 

chr\_BY, the number of the chromosome for BY

\item {} 

lower\_bound\_BY, the lower bound of the nucleosome for BY

\item {} 

upper\_bound\_BY, the upper bound of the nucleosome  for BY

\item {} 

index\_nuc\_BY, the index of the BY nucleosome in the CUR for BY

\item {} 

chr\_RM, the number of the chromosome for RM

\item {} 

lower\_bound\_RM, the lower bound of the nucleosome for RM

\item {} 

upper\_bound\_RM, the upper bound of the nucleosome  for RM

\item {} 

index\_nuc\_RM,the index of the RM nucleosome in the CUR for RM

\item {} 

cur\_index, index of the CUR

\item {} 

BY\_H3K14ac\_36, the number of reads for the current nucleosome for the sample 36

\item {} 

BY\_H3K14ac\_37, \#reads for sample 37

\item {} 

BY\_H3K14ac\_53, \#reads for sample 53

\item {} 

RM\_H3K14ac\_38, \#reads for sample 38

\item {} 

RM\_H3K14ac\_39, \#reads for sample 39

\end{itemize}

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

{\footnotesize

<<echo=TRUE, results=hide>>=

for (form in c("wp", "unr")) {

	for (marker in markers) {

    build_count_table(marker, combis[[1]], form, curs, all_samples, config)

  }

}			

@

<<echo=TRUE, results=verbatim>>=

a_count_table = read.table("outputs/BY_RM_Mnase_Seq_wp_and_nbreads.tab", header=TRUE)

head(a_count_table)

@

}

\subsection{Analysing Count Tables using DESeq}

Finally, the script \emph{launch\_deseq.R} perform statistical analysis on each nucleosome using \emph{DESeq}. It produces files:

\begin{itemize}

\item {} 

results/current/BY\_RM\_H3K14ac\_wp\_snep.tab

\item {} 

results/current/BY\_RM\_H3K14ac\_unr\_snep.tab

\item {} 

results/current/BY\_RM\_H3K14ac\_wpunr\_snep.tab

\end{itemize}

These files are organised with the following columns (see file \emph{BY\_RM\_H3K14ac\_wp\_snep.tab} for an example):

\begin{itemize}

\item {} 

chr\_BY, the number of the chromosome for BY

\item {} 

lower\_bound\_BY, the lower bound of the nucleosome for BY

\item {} 

upper\_bound\_BY, the upper bound of the nucleosome  for BY

\item {} 

index\_nuc\_BY, the index of the BY nucleosome in the CUR for BY

\item {} 

chr\_RM, the number of the chromosome for RM

\item {} 

lower\_bound\_RM, the lower bound of the nucleosome for RM

\item {} 

upper\_bound\_RM, the upper bound of the nucleosome  for RM

\item {} 

index\_nuc\_RM,the index of the RM nucleosome in the CUR for RM

\item {} 

cur\_index, index of the CUR

\item {} 

form

\item {} 

BY\_Mnase\_Seq\_1, the number of reads for the current nucleosome for the sample 1

\end{itemize}

Next columns concern indicators for each sample:

\begin{itemize}

\item {} 

BY\_Mnase\_Seq\_2, \#reads for sample 2

\item {} 

BY\_Mnase\_Seq\_3, \#reads for sample 3

\item {} 

RM\_Mnase\_Seq\_4, \#reads for sample 4

\item {} 

RM\_Mnase\_Seq\_5, \#reads for sample 5

\item {} 

RM\_Mnase\_Seq\_6, \#reads for sample 6

\item {} 

BY\_H3K14ac\_36, \#reads for sample 36

\item {} 

BY\_H3K14ac\_37, \#reads for sample 37

\item {} 

BY\_H3K14ac\_53, \#reads for sample 53

\item {} 

RM\_H3K14ac\_38, \#reads for sample 38

\item {} 

RM\_H3K14ac\_39, \#reads for sample 39

\end{itemize}

The 7 last columns concern DESeq analysis:

\begin{itemize}

\item {} 

manip{[}a\_manip{]} strain{[}a\_strain{]} manip{[}a\_strain{]}:strain{[}a\_strain{]}, the manip (marker) effect, the strain effect and the snep effect. These are the coefficients of the fitted generalized linear model.

\item {} 

pvalsGLM, the pvalue resulting from the comparison of the GLM model considering the interaction term \emph{marker:strain} to the GLM model that does not consider it. This is the statsitcial significance of the interaction term and therefore the statistical significance of the SNEP.

\item {} 

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.0001.

\item {} 

mnase\_l2fc, this is the coefficients of the fitted generalized linear model based on MNASE samples only.

\item {} 

mnase\_l2fc\_pval, the pvalue resulting from the analysis of the simple model based on MNASE samples only. This is the statsitcial significance of the \emph{strain} term and therefore the statistical significance of the strain effect.

\end{itemize}

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

{\footnotesize

<<echo=TRUE, results=hide>>=

for (form in c("wp", "unr", "wpunr")) {

  sneps = analyse_count_table(markers[2], combis[[1]], form, all_samples, config=config)

}

@

<<echo=TRUE, results=verbatim>>=

a_sneps_table = read.table("outputs/BY_RM_H3K14ac_wp_snep.tab", header=TRUE)

head(a_sneps_table)

@

}

Here are the number of computed SNEPs for each forms.

\end{document}