Homec4science
Diffusion PhD_thesis (master) Browse

R9524/main/ch_group_projects.tex
ef1db7ca3a32master

Last Change
Raw File

ch_group_projects.tex
View Options

\cleardoublepage
\chapter{Published laboratory projects}
\markboth{Published laboratory projects}{Published laboratory projects}
\addcontentsline{toc}{chapter}{Published laboratory projects}
This section contains the description of projects that have been accomplished as part of a laboratory-wide effort. Resources such as softwares and databases which have been subjected to a publication in peer-reviewed article are listed here.
\section{Mass Genome Annotation repository}
As part of an effort to enhance reproducibility and data re-use in genomics, my laboratory developed and maintain the Mass Genome Annotation repository (MGA). This repository is a resource that has been designed to store publicly available sequencing data released on GEO [REFERENCE] or ArrayExpress [REFERENCE] in a standardized manner with a high quality data annotation (metadata). This repository has been already been described in \citep{dreos_mga_2018}.
My laboratory developed and maintain several resources, among which several databases. One of them is Mass Genome Annotation repository (MGA). The MGA is a resource that has been designed to store publicly available sequencing data in a highly standardized manner to facilitate their re-use by the scientific community. My involvement in this project was related to data curration and annotation.
\subsection{Introduction}
owaday, in order to publish a study in a peer-reviewed article, the authors are requested to submit their sequencing data to primary repository, such as GEO [REFERENCE] or ArrayExpress [REFERENCE]. However, in both cases, the raw sequencing data are deposited and the metadata are often incomplete, difficult to understand if not wrong.
Because working with such data requires quality control checks and mapping, it is a substantial investment of work and time which is often a barrier impeding the re-use of the data and the reproducibility of published studies.
\subsection{MGA content and organization}
\begin{figure}[!htbp]
\begin{center}
\includegraphics[scale=0.8]{images/ch_group_projects/mga_figure1.jpeg}
\captionof{figure}{\textbf{Content of the MGA repository} \textbf{A} Proportion of samples in the database grouped by type. \textbf{B} Proportion of samples grouped by organism. Assemblies belonging to the same organism are merged together. \textbf{C} Samples numbers stratified by type and organism. Dot areas are proportional to the total number of samples in that category. The corresponding numbers can be found in a weakly updated table posted on the MGA home page at \url{http://ccg.vital-it.ch/mga}.
Figure and legend taken and adapted from \citep{dreos_mga_2018}.}
\label{lab_projects_mga_stats}
\end{center}
\end{figure}
The MGA contains 7 different types of data : ChIP-seq reads, ChIP-seq peaks, transcription profiling data, nucleosome/chromatin accessibility data, DNA methylation data, genome annotation data (TSSs, intro-exons boundaries) and sequence derived data (sequence conservation scores [REFERENCESSSS], genetic variants, predicted TF binding sites computed from binding models). The MGA contains more than 24'000 thousand samples in 15 different species (Figure \ref{lab_projects_mga_stats}).
All the data available on the MGA are stored in Simple Genome Annotation (SGA, [REFERENCE]). SGA is a single coordinate format. In essence, all data are represented as a single coordinate along the genome. However this does not compromise the accuracy of the data. Furthermore, if needed, SGA data can easily be converted to BED format data. Finally, different data type are represented differently : ChIP-seq peaks and paired-end sequenced fragments are represented by the middle position, reads are represented by their 5' end, TSSs are single base coordinates anyway, and so one. On the advantage to minimize the amount of data stored and to speed up their treatment by in-house developed tools (ChIP-cor, ChIP-extract, SSA-server tools, and so one).
TODO put ENCODE, Fantom5, RoadMap data somewhere
In order to enhance original results reproducibility as much as possible, processed data, if available on the primary repositories, are always preferred. Otherwise, the raw sequencing data are downloaded and processed using a general pipeline comprising i) read mapping using Bowtie [REFERENCE] or Bowtie 2 [REFERENCE], ii) conversion to bed [REFERENCE] and a final conversion to SGA [REFERENCE]. As on GEO, data (called samples) for a given study/article belong to a same serie. Finally, the metadata are documented. A full description of the data, their biological significance, their processing is available in HTML format. Two additional machine readable text files are available with i) the sample information ii) the serie information.
Importantly, it should be highlighted that the MGA is fully interconnected in-house developed analysis tools hosted on the ChIP-seq [REFERENCE] and SSA [REFERENCE] servers. All the data hosted can thus be analysed using any of these in-house developed tools.
\subsection{Conclusions}
The MGA repository is an important asset to the scientific community. It allows anybody to undertake quickly a wide range of data analyses, together with the ChIP-seq and SSA servers. Additionally, all these data are readily available and can be downloaded in MGA or BED format. Using them is so convenient that I used MGA hosted data in almost all my projects.
\section{Eukaryotic Promoter Database}
\section{PWMScan}
This section describes PWMScan, a software that has been developed to predict TF binding sites given a binding specificity model. This work has been conducted in close collaboration with Giovana Ambrosini, a senior scientist of the laboratory.
My contribution to this worked is substantial. Among other things, I contributed to the development of the matrix\_scan program, set up the entire benchmark procedure, run it and finally developed the necessary backend Python code for an efficient parallelization of matrix\_scan.
The content of this section is taken and adapted - with the author permission - from the original article \citep{ambrosini_pwmscan:_2018}. I produced all the figures presented in this section.
% Because TF regulate gene expression by binding to short specific DNA sequences of 5–20 bp, it is of prime importance to be able to predict the binding sites of a TF given its binding specificity.
% To fulfill this task, I participated in the development a genome scanner that uses a PWM as TF specificity model to scan a genome and returns candidate binding sites.
% A web-based version of this software has been made available to scan server-resident genomes. Additionally, the user can provide a PWM or select one among the different hosted public databases such as JASPAR \citep{khan_jaspar_2018}, HOCOMOCO \citep{kulakovskiy_hocomoco:_2018} and many others.
\subsection{Introduction}
% article Figure 1
% \begin{figure}[!htbp]
% \begin{center}
% \includegraphics[scale=0.4]{images/ch_group_projects/pwmscan_fig1.png}
% \captionof{figure}{Screen shot of the PWMScan graphical user interface.}
% \label{lab_projects_gui}
% \end{center}
% \end{figure}
% flowchart present in supplemental material
\begin{figure}
\begin{center}
\includegraphics[scale=0.1]{images/ch_group_projects/pwmscan_flowchart.png}
\captionof{figure}{\textbf{PWMScan workflow :} the input is composed of a PWM and a score threshold specifying the minimum score for a sequence to achieved to be considered as a match. Letter probability matrices or count matrices are also accepted and are converted into PWMs. The score threshold can also be given as a p-value or a percentage of the maximum score, in which case it is converted into a threshold score. Based on the length of the PWM, Bowtie or pwm\_scan can be used to find the matches on the genome. If Bowtie is used, the set of k-mers achieving a better score than the threshold score is computed using branch-and-bound algorithm (mba) and mapped on the genome. On the other hand, if matrix\_scan is used, the PWM is used to score every possible sub-sequence in the genome. The regions corresponding to the sequences achieving a score at least as good as the threshold score are then returned under BED format.
Figure and legend taken and adapted from \citep{ambrosini_pwmscan:_2018}.}
\label{lab_projects_pwmscan_pipeline}
\end{center}
\end{figure}
Knowing where transcription factors (TFs) bind to the genome is the key to understanding gene regulation. The binding specificity of a TF is commonly represented by a numerical matrix, either as a position weight matrix (PWM), a position frequency matrix (PFM) or a letter probability matrix (LPM). The three representations are information-wise equivalent and inter-convertible. A PWM contains weights for each base at each motif position. By summing up weights at corresponding positions, a binding score can be computed for any base sequence of the same length as the PWM. Large collections of TF specificity matrices are nowadays available frpublishedom public libraries such as JASPAR \citep{khan_jaspar_2018} or HOCOMOCO \citep{kulakovskiy_hocomoco:_2018}.
PWMScan is a web-server for rapid scanning of large genomes for high-scoring matches to a user-supplied or server-resident PWM. Compared to other web-based PWM scanning tools, PWMScan is unique in that it scans server-resident whole genomes rather than user-uploaded DNA sequences. Other key features are: (i) menu-driven access to genomes of >30 model organisms; (ii) menu-driven access to >300 public PWM libraries; (iii) support of various PWM representations and formats; (iv) cut-off values can be specified as match scores or P-values; (v) output in BEDdetail format with match scores and P-values; (vi) links to UCSC genome browser for visualization of results; and (vii) action buttons to transfer match lists to analysis tools.
A short description of the PWMScan server follows. Technical details about algorithms, programs and data are provided under Supplementary Material in the original article \citep{ambrosini_pwmscan:_2018}.
\subsection{Data and methods}
Genome sequences were downloaded from the NCBI in FASTA format. Indexed versions for rapid scanning were generated for Bowtie \citep{langmead_ultrafast_2009}. The motif databases offered by PWMScan have been downloaded from the MEME Suite website \citep{bailey_meme_2009}. LPMs have been converted to integer PWMs (see Supplementary Material in \citep{ambrosini_pwmscan:_2018}).
% The input form of the PWMScan server is shown in Figure \ref{lab_projects_pwmscan_gui}. The user chooses a genome assembly from a menu. Optionally, a BED file may be uploaded to restrict the search to genomic regions of particular interest, e.g. open chromatin regions. The right side of the form offers several ways to specify a DNA motif. PWMs from a server-resident database are chosen from a pull-down menu. Alternatively, matrices can be entered into a text area or uploaded. Accepted motif types are: PFMs, LPMs, real or integer PWMs and IUPAC consensus sequences. PFMs can be entered in several formats, including TRANSFAC and JASPAR.
All motif types have to be converted into integer PWMs for input to the genome search engines (see Supplementary Material in \citep{ambrosini_pwmscan:_2018}). Default conversion parameters are proposed and can be changed by the user.
For instance, real PWMs can be rescaled on input by a multiplication factor to ensure sufficient resolution after integer conversion. IUPAC consensus sequences are converted into binary matrices consisting of 0 and 1.
For all matrix formats, the cut-off value can be specified as PWM score, as P-value or as percentage of the score range (0\% = minimal score and 100\% = maximal score). For IUPAC consensus sequences, the cut-off value is specified as a maximal number of mismatches allowed.
The P-value of a PWM score x is defined as the probability that a random k-mer sequence of the length of the PWM has a binding score $\geq$ x given the base composition of the genome.
The whole genome scan takes as input an integer PWM and a corresponding cut-off. The output is a list of sequence regions that match the PWM with a match score higher or equal to the cut-off value. Depending on the length of the PWM and the cut-off, one of the following search strategies is chosen: (i) Bowtie, a fast memory-efficient short read aligner using indexed genomes and (ii) matrix\_scan, a C program developed by our group using a conventional search algorithm, as shown in Figure \ref{lab_projects_pwmscan_pipeline}.
The first strategy is more efficient for short PWMs and high cut-off values. It requires as a first step the generation of a list of all k-mers that match the PWM with the given cut-off. The list of k-mers is then mapped to the genome using Bowtie.
The second strategy takes genome sequences in FASTA format as input. Individual chromosomes are processed in parallel and distributed to multiple cores by a Python script. We empirically found that this approach becomes more efficient if the number of k-mers exceeds 10 5 sequences. matrix\_scan was benchmarked for speed together with five other matrix scanners and was found to be the fastest (see Supplementary Material in \citep{ambrosini_pwmscan:_2018}).
The basic search step outputs a list of PWM matches, including the genomic coordinates, the DNA sequence and the match score. Post-processing of this list involves computation of the corresponding P-values, addition of the matrix name and, optionally, elimination of overlapping matches. The final match list is provided in BEDdetail format. The output page further shows the total number of PWM matches and a sequence logo reflecting the letter-probabilities of the input matrix. Action buttons are provided for: (i) sending the match lists to analysis tools of the ChIP-Seq and SSA servers \citep{ambrosini_chip-seq_2016}, (ii) extracting DNA sequences around the matches, (iii) sending the output to the UCSC genome browser for visualization and (iv) liftover of the match coordinates to other assemblies of the same or related species.
PWMScan is meant to support many types of genomic data analysis and designed to be interoperable with other tools from our group and elsewhere. An example of a typical workflow involving ChIP-seq data is presented in Supplementary Material in \citep{ambrosini_pwmscan:_2018}.
PWMScan is also available as a command-line software package from SourceForge, including a master script scheduling all computational steps running during a web job.
\subsection{Benchmark}
% supplemental figure 3 in article
\begin{figure}
\begin{center}
\includegraphics[scale=0.2]{images/ch_group_projects/pwmscan_figure_s1.png}
\captionof{figure}{\textbf{Benchmark :} PWMScan speed performances were measured and compared with 6 other well known genome scanners. In all cases, the h19 genome sequence was scanned with a 19bp CTCF matrix and a 11bp STAT1 matrix, 10 times. The run times are represented as boxplots. For PWMScan, both pwm\_scan and Bowtie strategies were run.
Figure and legend taken and adapted from \citep{ambrosini_pwmscan:_2018}.}
\label{lab_projects_pwmscan_benchmark}
\end{center}
\end{figure}
\begin{table}
\begin{center}
\begin{tabular}{| c | c | c | c | c | c | c |}
\hline
Program & \multicolumn{3}{|c|}{CTCF} & \multicolumn{3}{|c|}{STAT1} \\
\hline
& \thead{Nb of\\matches} & \thead{Perc.\\overlap} & \thead{Score corr.} & \thead{Nb of\\matches} & \thead{Perc.\\overlap} & \thead{Score\\corr}. \\
\hline
matrix\_scan & 70607 & 100 & 1 & 138531 & 100 & 1 \\
\hline
Patser & 70607 & 100 & 1 & 138531 & 100 & 1 \\
\hline
PossumSearch & 70607 & 100 & 1 & 138531 & 100 & 1 \\
\hline
FIMO & 70131 & 99.16 & 0.9989 & 134901 & 100 & 1 \\
\hline
RSAT & 70701 & 99.85 & 0.9999 & 134901 & 100 & 1 \\
\hline
HOMER & 70701 & 99.95 & 0.9999 & 134901 & 100 & 1 \\
\hline
STORM & 70704 & 99.83 & 0.9999 & 134901 & 100 & 0.9999 \\
\hline
\end{tabular}
\caption{\textbf{Motif scanning software comparison}. The performances of matrix\_scan were assessed by comparing how many of the regions listed by matrix\_scan were also returned by other programs and if the region scores were comparable. For the percentage of overlap with the match list returned by matrix\_scan, the shorter of the two lists always serves as the reference (100\%). For the score correlations with matrix\_scan scores, the Spearman correlation was used.
Table and legend taken and adapted from \citep{ambrosini_pwmscan:_2018}.}
\label{lab_projects_pwmscan_benchmark_table}
\end{center}
\end{table}
% article Table 1
% \begin{center}
% \captionof{table}{Speed is expressed in seconds. The benchmarking tests have been run on a Linux/CentOS7/x86\_64 workstation with 48 CPU-cores and 256 GB of DRAM.
% The performance measurements for matrix\_scan were achieved using 1 (left value) and 10 (right value) CPU-cores in parallel.}
% \begin{tabular}{ c c c }
% PWM/P-value & Bowtie speed & Matrix\_scan speed \\
% \hline
% STAT1(len = 11 bp)/$10^{-5}$ & 3 & 30/5 \\
% STAT1(len = 11 bp)/$10^{-4}$ & 8 & 40/8 \\
% STAT1(len = 11 bp)/$10^{-3}$ & 60 & 65/30 \\
% CTCF(len = 19 bp)/$10^{-5}$ & 12 & 40/6 \\
% CTCF(len = 19 bp)/$10^{-4}$ & 90 & 50/10 \\
% CTCF(len = 19 bp)/$10^{-3}$ & 720 & 90/35 \\
% \label{lab_projects_pwmscan_table}
% \end{tabular}
% \end{center}
The problem of pattern-matching using a matrix is considered an important problem and many algorithms and program have been developed. The challenge comes from the fact that this problem should be solved in a time and memory efficient manner but also because it poses a more fundamental problem : most of the prediction turn to be false positives.
The most straight forward approach to this problem is the sliding window approach. This way of doing is implemented by matrix\_scan, Patser \citep{hertz_identification_1990}, PossuMSearch \citep{beckstette_fast_2006}, RSAT matrix-scan (later on simply referred to as RSAT, \citep{turatsinze_using_2008}), HOMER \citep{heinz_simple_2010} or FIMO \citep{grant_fimo:_2011}. PossuMSearch \citep{beckstette_fast_2006} and matrix\_scan implements "lookahead search", allowing to eliminate sub-sequences that do not match the matrix faster than a regular scoring method. PossuMSearch also propose an alternative search strategy, based on the use of enhanced suffix arrays (ESA). This is a two step method. First the genome is indexed, to stored it using ESA and dumped to file. Then, the ESA are searched for matches using a matrix, allowing a sub-linear search time.
These programs also generally compute a p-value associated with each score $s$. The p-value is the probability that a random sequence get at least a score as good as $s$. FIMO also perform an expensive False Discovery Rate (FDR) estimate in order to correct the pvalues for multiple testing. On its side, STORM \citep{schones_statistical_2007} uses a database of regulatory sequences to count the expected number of occurrences of kmers to compute a p-value. STORM creates a so called ($g$,$k$) table in order to store the number of occurrences of all possible kmers of size $k$ with a maximum gap of size $g$ present in a sequence database.
Because both statistical and run-time efficiency aspects of algorithms have been subject of intense work, we decided to assessed the performances of PWMScan both in term of speed than in terms of scoring similarity with other well established methods.
Also, for completness, it is worth mentioning that other programs have been released, such as MotifScanner \citep{aerts_toucan:_2003}, MotifViz \citep{fu_motifviz:_2004} or TRED \citep{zhao_tred:_2005}. However they could not be included in the benchmark for diverse reasons. MotifScanner could only be used through a GUI client connected to a (offline?) server, making it unsable for batch analyses. MotifViz relied on deprecated 32bits libraries and TRED input sequence size was limited to 10kb.
The runtimes of PWMScan and of other competitor programs were measured by scanning the human genome (UCSC assembly hg19) with two different PWMs from JASPAR, STAT1 with length 11 bp and CTCF with length 19 bp, and different cut-off values expressed as P-values. Results are shown in Figure \ref{lab_projects_pwmscan_benchmark} and and Table \ref{lab_projects_pwmscan_benchmark_table}. PoSSuMSearch was used with its sliding window strategy because to use of ESA, on an already indexed genome, turned out to be slower (not shown). Also, note that for longer motifs and higher P-values, the Bowtie-based approach becomes inefficient, whereas matrix\_scan remains reasonably fast (not shown).
\subsection{Conclusions}
We developed PWMScan, a program that searches matches inside a genome using a PWM. We implemented two alternative search strategies : i) a traditional sliding window approach and ii) an innovative matching sub-sequences mapping using Bowtie. Both strategies are competitive in terms of speed and score similarity compared to other existing programs. Furthermore, our window sliding approach turned out to be fastest of all.
Additionally, PWMScan is available as a standalone distribution as well as through a dedicated web-server with access to dozens of species genomes and thousand of PWMs from the most popular databases.
To conclude, PWMScan is an asset for the research community because of its simplicity of use and of its high performances.
\clearpage
\section{SPar-K}
This section describes SPar-K (Signal Partitioning with K-means), a modification of the K-means algorithm to cluster genomic regions based on their chromatin organization, defined by by their sequencing profiles.
My contribution to this work is major. I developed the algorithm, implemented it and run all the benchmark and other analyses.
The content of this section is taken an adapted from the original article \citep{groux_spar-k:_2019}. I produced all the figures presented in this section.
\subsection{Introduction}
The sequence-specific binding of transcription factors (TFs) to target sites and the associated changes in chromatin state are considered key events in gene regulatory processes that turn on or off genes. The advent of powerful assays to characterize chromatin features genome-wide (e.g. ChIP-seq, DNase-seq, MNase-seq) and the large-scale application of these technologies to many cell types by consortia such as ENCODE \citep{consortium_integrated_2012} have created a wealth of data. To gain knowledge from such data, it is common to analyze positional correlations between chromatin features, e.g. the position of nucleosomes (revealed by MNase-seq) relative to transcription factor binding regions (mapped by ChIP-seq). However, as noted in \citep{kundaje_ubiquitous_2012}, chromatin patterns tend to be heterogeneous and often asymmetric. The latter is a problem if the orientation of a given functional region is not defined by experiments. Furthermore, limited mapping precision leading to moderate misalignment between functionally equivalent regions can also obscure a chromatin pattern.\\
Several methods and software tools have been developed for discovering chromatin pattern by clustering and/or realignment of signal profiles for genomic regions, including ChromaSig, CATCHProfiles, CAGT and ChIP-Partitioning, see \citep{nair_probabilistic_2014} and references therein. All of these programs have some limitations. Some do not realign, others are restricted to count data or lack an efficient implementation. To fill this gap, we developed SPar-K (Signal Partitioning with K-means).
% figure 1 from article
\begin{figure}
\centerline{\includegraphics[scale=0.4]{images/ch_group_projects/spark_figure1.pdf}}
\caption{Partitioning of DNaseI hypersensitivity profiles around SP1 binding sites in K562 cells. The optimal number of clusters was determined using the elbow method (Figure \ref{fig_s07}). \textbf{A.} Input data based on peak summits provided by ENCODE. \textbf{B.} Same regions clustered, re-aligned and oriented by SPar-K. Clusters 1, 2 and 3 are indicated by colored bars in red, blue, and green, respectively. \textbf{C.} MNase-seq read densities for the same regions, ordered, aligned and oriented as in B. \textbf{D.} Predicted SP1 binding motifs for the same regions, ordered, aligned and oriented as in B. \textbf{E.} Proportion of binding sites within each cluster having a confirmed promoter-associated TSS within +/- 300bp. \textbf{F.} Aggregations profiles for DNase-seq (red), MNase-seq (blue), promoter TSS (green) and CAGE-seq data (violet) for cluster 2 (aligned and oriented as in B). \textbf{G.} Motifs found by MEME-ChIP and Tomtom in the narrow footprints of each cluster. (*) known SP1 interactor, (c) central enrichment. Cluster 2 left and right refer to the left and right footprints seen in \textbf{B}.
Figure and legend taken and adapted from \citep{groux_spar-k:_2019}.}
\label{lab_projects_spark_dnase}
\end{figure}
\subsection{Methods}
SPar-K is a modified K-means algorithm which iteratively partitions, aligns and orients a set of genomic regions based on their signal density profiles obtained by a chromatin profiling assay (or in some other way). SPar-K is implemented as a stand-alone, fully multithreaded, C++ program. A brief description of its overall design follows. Details can be found in the Supplementary Material where all the algorithms are described.\\
Input data are the signal over $N$ regions of length $L$ stored as a $N$ rows and $L$ columns matrix. The signal resolution may be at single-base or at a larger bin size. The regions are typically defined by relative positions to an anchor point, e.g. a ChIP-seq peak summit.\\
The user sets via command-line options the number of clusters to find $K$, the shifting freedom $S$, whether flipping is allowed, the maximum number of iterations, whether outlier smoothing should be performed, the seeding strategy, and optionally a seed for the random number generator. The shifting freedom is specified as the maximal number of positions two regions may be shifted relative to each other ($S$), or pattern width ($W$=$L$-$S-1$). The partitioning is initialized by choosing $K$ rows from the input matrix according to the specified seeding strategy. These regions define the initial class centers (signal profiles). During iterative optimization, each region is compared to the current class centers in all possible alignments (shifting states) and orientations and assigned to the closest one in terms of correlation distance. The class centers are then updated by aggregation of the optimally aligned regions of its members. The process stops after a convergence criterion or the maximal number of iterations is reached.\\
The program returns a table listing for each region the cluster assignment, the shift state and the orientation. The software distribution includes R scripts for visualizing the data as heatmaps as shown in Figure \ref{lab_projects_spark_dnase}.
% supplemental figure 2 from article
\begin{figure}
\centerline{\includegraphics[scale=0.4]{images/ch_group_projects/spark_supplemental_figure2.pdf}}
\caption{\textbf{Clustering accuracy using random seeding :} to compare the clustering accuracies of the different methods, several simulated dataset containing 3 classes, different coverages (10, 50 and 100 reads per region indicated as "cov10", "cov50" and "cov100") and noise proportions (no noise, 10\% noise, 50\% noise and 90\% noise indicated as "0.0", "0.1", "0.5" and "0.9") were generated. Each dataset was clustered 50 times with each method. The Adjusted Rand Index (ARI) was computed for each partition. The ARI values are displayed as boxplots. SPar-K and ChIPPartitioning were run allowing flipping and shifting. The ARI was measured on each of the resulting data partitions. For SPar-K, "smooth" indicates outlier smoothing. For the regular K-means, "eucl." and "corr." refer to the euclidean and correlation distances. "R" stands for "random" and indicates the ARI values obtained when comparing the true cluster labels with a randomly shuffled version of it, 100 times.
Figure and legend taken and adapted from \citep{groux_spar-k:_2019}.}
\label{lab_projects_spark_ari}
\end{figure}
% supplemental figure 4 from article
\begin{figure}
\centerline{\includegraphics[scale=0.4]{images/ch_group_projects/spark_supplemental_figure4.pdf}}
\caption{\textbf{Median SSE :} for the simulated ChIP-seq dataset containing 3 classes, with coverage 100 and no noise, partitioned into 2 to 5 clusters. To judge whether the elbow method could be used to estimate the optimal number of clusters, this dataset was partitioned with SPar-K, allowing flip and shifting, into 2 to 5 clusters, 50 times for each set of parameters. For each number of clusters, the median SSE is shown, +/- 1 standard deviation (bars). \textbf{A} Seeding done at random, \textbf{B} seeding done at random and outlier smoothing \textbf{C} seeding done with the K-means++ method \textbf{D} seeding done with the K-means++ method and outlier smoothing. In all cases, the optimal number of clusters seemed to be 3 (which was the expected value).
Figure and legend taken and adapted from \citep{groux_spar-k:_2019}.}
\label{lab_projects_spark_sse}
\end{figure}
% supplemental figure 4 from article
\begin{figure}
\centerline{\includegraphics[scale=0.4]{images/ch_group_projects/spark_supplemental_figure5.pdf}}
\caption{\textbf{Running times :} to compare the run times of each program, the synthetic dataset with coverage 100 and no noise was partitioned 20 times with each program. The run times (wall clock) in second were measured. For all SPar-K and the regular K-means, the partitions were initialized using a random and K-means++ (indicated as "k++"). For ChIPPartitioning, only a random seeding was used. The partitions were then optimized for 30 iterations at most. For SPar-K and ChIPPartitioning, a shifting of 71 bins and flipping were allowed. For SPar-K, only one thread was used and "smooth" indicates outlier smoothing. For the regular K-means, "eucl." and "corr." refer to the euclidean and correlation distances.
Figure and legend taken and adapted from \citep{groux_spar-k:_2019}.}
\label{lab_projects_spark_time}
\end{figure}
\subsection{Results}
First we compared SPar-K, regular K-means and ChIP-partitioning on synthetic datasets (see Supplemental Figure 1 in \citep{groux_spar-k:_2019}) exhibiting properties that are plausible for ChIP-seq profiles for genomic regions: count vectors, 3 classes with different shapes, limited misalignment, random orientation and various amounts of background noise. Performance was assessed by the Adjusted Rand index (see Figure \ref{lab_projects_spark_ari} and Supplemental Figure 3 in \citep{groux_spar-k:_2019}) and the optimal number of classes was estimated by the elbow method (Figure \ref{lab_projects_spark_sse}). As expected, regular K-means performed poorly. On the contrary, SPar-K was equally accurate as ChIP-Partitioning except for the lowest coverage class. Considering speed, Spar-K outperformed ChIP-partitioning by a factor of at least 20 (Figure \ref{lab_projects_spark_time}).\\
We applied SPar-K with $K=3$ to DNaseI accessibility profiles (2bp resolution) around 7'206 ChIP-seq SP1-binding peaks (+/-300bp relative to peak summit) in K562 cells (Figure \ref{lab_projects_spark_dnase}A). The results revealed the presence of clear footprints in all the clusters (Figure \ref{lab_projects_spark_dnase}B). To validate these footprints, we checked whether they are consistent with to location of nucleosomes (Figure \ref{lab_projects_spark_dnase}C) and SP1 binding motifs (Figure \ref{lab_projects_spark_dnase}D), which was indeed the case. De novo motif analysis of the narrow footprints seen in Figure \ref{lab_projects_spark_dnase}B with MEME-ChIP and Tomtom \citep{bailey_meme_2009} discovered SP1-related, NFYA/B and GATA motifs (Figure \ref{lab_projects_spark_dnase}G) the latter two reportedly being interaction partners of SP1. Taken together, these results suggest that SPar-K is able to precisely refocus initially misaligned DNaseI profiles around SP1 binding sites.\\
The partitioning of SP1 binding regions reveals distinct chromatin landscapes. Cluster 1 (red) groups binding sites lying between two closely spaced nucleosomes. Cluster 2 (blue) shows strong asymmetry suggestive of promoter regions, an interpretation supported by the presence of promoter-associated transcription starts sites (TSSs) and CAGE tags (Figure \ref{lab_projects_spark_dnase}E and F). Finally, the symmetrical cluster 3 (green) contains binding sites located on large nucleosome-free regions reminiscent of enhancer regions.\\
As a second example, we ran the same type of analysis on nucleosome profiles around CTCF binding sites (see Supplemental Figure 8 in \cite{groux_spar-k:_2019}). Overall, the results confirm observations published by \citep{kundaje_ubiquitous_2012}. Strong nucleosome arrays became visible in all classes after realignment, with three out of four showing strong asymmetry in addition.
\subsection{Conclusion}
SPar-K is a useful partitioning method for moderately misaligned and randomly oriented chromatin regions. Compared to existing methods, it is competitive in terms of accuracy, superior in speed, applicable to a wider range of input signals (not restricted to count data) and easier to use.
c4science · Help