chapter_3.tex

\chapter{Реализация и использование}\label{implementation}
В этой главе будут подробно описаны реализованные методы анализа в проекте \emph{phantasus}, их реализация на стороне сервера и на стороне клиента, будет рассказано о способах запуска приложения и другие инфраструктурные подробности.
\section{Реализованные методы анализа экспрессии}
В ходе работы над проектом были реализованы следующие функции и методы анализа экспрессии:\begin{itemize}
\item \texttt{loadGEO} --- загрузка и визуализация данных из \emph{Gene Expression Omnibus};
\item \texttt{pcaPlot} --- реализация метода главных компонент и визуализация результата;
\item \texttt{kmeans} --- реализация кластеризации методом \emph{kmeans};
\item \texttt{limmaAnalysis} --- анализ дифференциальной экспрессии для сравнения образцов. \end{itemize}

В этом разделе будут один за другим описаны все вышеперечисленные методы.

\subsection{Загрузка данных из GEO}
В разделе~\ref{geointro} обзора были описаны форматы данных в репозитории \emph{Gene Expression Omnibus}.
В \emph{phantasus} загрузка данных из \emph{GEO} осуществляется следующим образом:
\begin{enumerate}
\item функция \texttt{loadGEO} принимает на вход идентификатор \emph{GEO};
\item в зависимости от его вида (\emph{GSE} или \emph{GDS}) запускаются дополнительные функции (\texttt{getGSE} на листинге~\ref{getGSE} и \texttt{getGDS} на листинге~\ref{getGDS});
\item в каждой из двух функций с помощью \texttt{GEOquery::getGEO} загружаются данные с аннотацией (или подгружаются из кэша, если он указан или если их уже загружали);
\item результат обрабатывается, создается \texttt{ExpressionSet} и отправляется в глобальные переменные;
\item в файл записываются сериализованные в \emph{ProtoBuf} данные, в том же формате, что и при создании \texttt{ExpressionSet} из внешних данных (смотри раздел~\ref{createESsection}), которые после считает и обработает клиент.
\end{enumerate}
\begin{lstlisting}[float=!h,caption={Загрузка данных типа GSE из Gene Expression Omnibus},label={getGSE},language=R]
getGSE <- function(name, destdir = tempdir()) {
  es <- getGEO(name, AnnotGPL = T, destdir = destdir)[[1]]
  featureData(es) <- featureData(es)[,grepl("symbol", fvarLabels(es), ignore.case = T)]
  phenoData(es) <- phenoData(es)[,grepl("characteristics", varLabels(es), ignore.case = T)
                                  | (varLabels(es) %in% c("title", "id", "geo_accession"))]
  chr <- varLabels(es)[grepl("characteristics", varLabels(es), ignore.case = T)]
  take <- function(x, n) {
    sapply(x, function(x) { x[[n]] })
  }
  rename <- function(prevName, x) {
    splitted <- strsplit(x, ": ")
    sumlength <- sum(sapply(as.vector(splitted), length))
    if (sumlength != 2 * length(x)) {
       return(list(name = prevName, x = x))
    }
    splittedFirst <- unique(take(splitted, 1))
    if (length(splittedFirst) == 1) {
       res = list(name = splittedFirst[1], x = take(splitted, 2))
    }
    else {
       res = list(name = prevName, x = x)
    }
    res
  }
  renamed <- lapply(chr, function(x) { rename(x, as.vector(pData(es)[,x])) })
  phenoData(es) <- phenoData(es)[, !(varLabels(es) %in% chr)]
  pData(es)[,take(renamed,1)] <- take(renamed,2)
  es
}
\end{lstlisting}

\begin{lstlisting}[float=!h,caption={Загрузка данных типа GDS из Gene Expression Omnibus},label={getGDS},language=R]
 getGDS <- function(name, destdir = tempdir()) {
  l <- getGEO(name, destdir = destdir)
  table <- slot(l, 'dataTable') # extracting all useful information on dataset
  data <- Table(table)  # extracting table ID_REF | IDENTIFIER/SAMPLE | SAMPLE1 | ...
  columnsMeta <- Columns(table) # phenoData
  sampleNames <- as.vector(columnsMeta[["sample"]])
  rownames <- as.vector(data[["ID_REF"]])
  symbol <- as.vector(data[["IDENTIFIER"]])
  data <- data[,sampleNames] # expression data
  exprs <- as.matrix(data)
  row.names(exprs) <- rownames
  row.names(columnsMeta) <- sampleNames
  # columnsMeta <- columnsMeta[,!(colnames(columnsMeta) %in% c('sample'))] 
  pData <- AnnotatedDataFrame(data.frame(columnsMeta, check.names = F))
  fData <- data.frame(matrix(symbol, nrow(exprs), 1));
  colnames(fData) <- "symbol"
  fData <- AnnotatedDataFrame(fData)
  featureNames(fData) <- rownames
  ExpressionSet(assayData = exprs, phenoData = pData, featureData = fData)
}
\end{lstlisting}

Основной код загрузки и дополнительные утилиты к нему можно увидеть на листинге~\ref{loadGEO}.

\begin{lstlisting}[float=!h,caption={Загрузка данных из Gene Expression Omnibus},label={loadGEO},language=R]
loadGEO <- function(name, type = NA) {
  es <- getES(name, type, destdir = "/var/phantasus/cache")
  assign("es", es, envir = parent.frame())
  data <- as.matrix(exprs(es)); colnames(data) <- NULL; row.names(data) <- NULL

  pdata <- as.matrix(pData(es)); colnames(pdata) <- NULL; row.names(pdata) <- NULL

  participants <- colnames(es)
  rownames <- rownames(es)

  fdata <- as.matrix(fData(es))
  colnames(fdata) <- NULL
  row.names(fdata) <- NULL

  res <- list(data = data, pdata = pdata,
              fdata = fdata, rownames = rownames,
              colMetaNames = varLabels(phenoData(es)),
              rowMetaNames = varLabels(featureData(es)))

  f <- tempfile(pattern = "gse", tmpdir = getwd(), fileext = ".bin")
  writeBin(protolite::serialize_pb(res), f)
  f
}
getES <- function(name, type = NA, destdir = tempdir()) {
  if (is.na(type)) {
     type = substr(name, 1, 3)
  }
  if (type == 'GSE') {
    es <- getGSE(name, destdir)
  }
  else if (type == 'GDS') {
    es <- getGDS(name, destdir)
  }
  else {
    stop("Incorrect name or type of the dataset")
  } 
  es
}
\end{lstlisting}

Такая реализация загрузки данных из \emph{GEO} позволяет избежать проблемы \emph{cross-origin request}, которая служила препятствием для загрузки из \emph{GEO} в исходном приложении \emph{Morpheus}. 


\subsection{Метод главных компонент и визуализация его результата}
Данный инструмент предназначен для построения графиков в соответствие с методом главных компонент.
В качестве аргументов на вход к инструменту подается:
\begin{itemize}
\item номера образцов для сравнения;
\item категориальная аннотация для различения точек по цвету (если не указана, то стандартный цвет);
\item числовая аннотация для различения точек по размеру (если не указана, то стандартный размер);
\item аннотация для подписей к точкам (если не указана, то без подписи);
\item функция замены \texttt{NA} в данных при вычислении матрицы \emph{PCA} (\texttt{mean} или \texttt{median}).
\end{itemize}

\begin{figure}[h]
  \caption{Инструмент \texttt{PcaPlotTool} и отрисованный график по данным GSE14308}
  \includegraphics[scale=0.4]{plotexample.png}
  \label{plotexample}
\end{figure}

Далее по алгоритму, описанному в разделе~\ref{functioncallalgo}, на \emph{OpenCPU}-сервер отправляется \emph{RPC}-вызов с аргументами: ключ сессии, содержащий актуальный \texttt{ExpressionSet}, и функция замены \texttt{NA}.

Данные аргументы приходят на вход к функции \texttt{pcaPlot}, реализованной в \emph{R}-пакете \emph{phantasus}, код которой можно увидеть на листинге~\ref{pcaPlot}. Предварительно все \texttt{NA}-значения заменяются в соответствии с переданной функцией, если этого не сделать, то дальнейшие вычисления будут невозможны. Матрица экспрессии из входного \emph{ExpressionSet} передается в стандартную функцию \texttt{prcomp} из \emph{R}-пакета \emph{stats}~\cite{stats}, которая и вычисляет результирующую матрицу.
\begin{lstlisting}[float=!h,caption={Вычисление матрицы главных компонент},label={pcaPlot},language=R]
  pcaPlot <- function(es, rows=c(), columns = c(), replacena = "mean") {
  rows <- getIndicesVector(rows, nrow(exprs(es)))
  columns <- getIndicesVector(columns, ncol(exprs(es)))
  data <- data.frame(exprs(es))[rows, columns]

  ind <- which(is.na(data), arr.ind = T)
  if (nrow(ind) > 0) {
    data[ind] <- apply(data, 1, replacena, na.rm = T)[ind[,1]]
  }
  ind1 <- which(!is.nan(as.matrix(data)), arr.ind = T)
  left.rows <- unique(ind1[,"row"])
  data <- data[left.rows,]
  data <- t(data)

  pca <- prcomp(data)
  explained <- (pca$sdev)^2 / sum(pca$sdev^2)

  xs <- sprintf("PC%s", seq_along(explained))
  xlabs <- sprintf("%s (%.1f%%)", xs, explained * 100)

  pca.res <- as.matrix(pca$x); colnames(pca.res) <- NULL; row.names(pca.res) <- NULL
  return(jsonlite::toJSON(list(pca = t(pca.res), xlabs=xlabs)))
}
\end{lstlisting}

На клиент в \emph{JSON}-формате приходит вычисленная матрица \emph{PCA}.

После, по дополнительным аргументам и вычисленной матрице, строится интерактивный график с помощью \emph{plotly.js}, пример которого можно увидеть на рисунке~\ref{plotexample}.

\subsection{Кластеризация методом kmeans}
Этот инструмент осуществляет разбиение генов на указанное пользователем число кластеров по алгоритму \emph{kmeans}.

\begin{figure}[h]
  \caption{Графический интерфейс инструмента \texttt{KmeansTool}}
  \includegraphics[scale=0.4]{kmeanstool.png}
   \label{kmeanstool}
\end{figure}

На клиенте в инструменте \texttt{KmeansTool}, который показан на рисунке~\ref{kmeanstool}, считываются следующие аргументы:
\begin{itemize}
\item количество кластеров, на которые нужно разбить данные;
\item функция замены \emph{NA} в данных.
\end{itemize}

Данные аргументы и ключ сессии актуального \texttt{ExpressionSet} отправляются на сервер в соответствующую функцию \texttt{kmeans}, код которой можно увидеть на листинге~\ref{kmeans}. Также как и в \texttt{pcaPlot}, здесь сначала заменяются \texttt{NA} на соответствующие данной функции значения, так как в этом случае \texttt{NA}-значения мешают вычислить среднее среди векторов. После данные отправляются в стандартную функцию \texttt{kmeans} пакета \emph{stats}~\cite{stats}. Результат возвращается как список соответствия каждого гена определенному кластеру, который приходит на клиент в \emph{JSON}-формате и отрисовывается как новая цветовая аннотация к строкам. Пример такой аннотации можно увидеть на рисунке~\ref{kmeansexample}.

\begin{lstlisting}[float=!h,caption={Кластеризация методом kmeans},label={kmeans},language=R]
kmeans <- function(es, columns = c(), rows = c(), k, replacena = "mean") {
  assertthat::assert_that(k > 0)

  rows <- getIndicesVector(rows, nrow(exprs(es)))
  columns <- getIndicesVector(columns, ncol(exprs(es)))
  data <- replacenas(data.frame(exprs(es))[rows, columns], replacena)

  data <- t(scale(t(data)))
  while (sum(is.na(data)) > 0) {
    data <- replacenas()
    data <- t(scale(t(data)))
  }

  km <- stats::kmeans(data, k)
  res <- data.frame(row.names = row.names(exprs(es)))
  res[["cluster"]] <- NA
  res[names(km$cluster), "cluster"] <- as.vector(km$cluster)
  return(toJSON(as.vector(km$cluster)))
}
\end{lstlisting}

\begin{figure}[h]
  \caption{Результат работы инструмена \texttt{KmeansTool} на данных GSE14308}
  \includegraphics[scale=0.4]{kmeansexample.png}
  \label{kmeansexample}
\end{figure}

\subsection{Анализ дифференциальной экспрессии}
Инструмент предназначен для анализа дифференциальной экспрессии: экспрессия генов сравнивается в двух группах образцов, и вычисляются несколько статистических характеристик, показывающих, насколько случайны различия этих групп.

\begin{figure}[h]
  \caption{Графический интерфейс инструмента \texttt{LimmaTool}}
  \includegraphics[scale=0.4]{limmatool.png}
  \label{limmatool}
\end{figure}

На клиенте, в инструменте, показанном на рисунке~\ref{limmatool}, осуществляется получение следующих аргументов:
\begin{itemize}
\item какие аннотации образцов участвуют в сравнении;
\item какая комбинация значений указанных выше аннотаций обозначает класс A;
\item аналогично для класса B.
\end{itemize}

Далее происходит подготовка аргументов к отправке на сервер: образцы разбиваются по выбранным аннотациям на три группы: A, B и не участвующие в сравнении.
Список соответствия образцов классам и ключ сессии, содержащей актуальный \texttt{ExpressionSet}, отправляются на сервер.

Аргументы приходят в функцию \texttt{limmaAnalysis}, код которой можно увидеть на листинге~\ref{limmaAnalysis}. Прежде чем использовать функцию \texttt{de} (\emph{differential expression}), реализованную в пакете \emph{limma}~\cite{limma}, которая помогает увидеть, насколько случайны различия между образцами, гены которых находятся в разных условиях, необходимо дополнить аннотацию образцов списком, содержащим идентификаторы классов сравнения.

Построенный дизайн сравнения передается в функцию \texttt{de}, которая возвращает матрицу статистических характеристик к каждому гену.
Эта матрица далее сериализуется в \emph{ProtoBuf} и записывается в файл.

\begin{figure}[h]
  \caption{Результат работы инструмента \texttt{LimmaTool} на данных GSE14308}
  \includegraphics[scale=0.4]{limmaexample.png}
  \label{limmaexample}
\end{figure}

Клиент, получив ключ временной \emph{OpenCPU}-сессии, читает файл с сериализованной в \emph{ProtoBuf} матрицей результатов, которые с помощью \emph{protobuf.js} разбираются и после отрисовываются в виде аннотации к строкам. Результат работы можно увидеть на рисунке~\ref{limmaexample}.

\begin{lstlisting}[float=!h,caption={Реализация анализа дифференциальной экспрессии в R-пакете phantasus},label={limmaAnalysis},language=R]
limmaAnalysis <- function(es, rows = c(), columns = c(), fieldValues) {
  assertthat::assert_that(length(columns) == length(fieldValues) || length(columns) == 0)
  rows <- getIndicesVector(rows, nrow(exprs(es)))
  columns <- getIndicesVector(columns, ncol(exprs(es)))
  fieldName <- "Comparison"
  fieldValues <- replace(fieldValues, fieldValues == '', NA)
  new.pdata <- pData(es)[columns,]
  new.pdata[[fieldName]] <- as.factor(fieldValues)
  new.pdata <- new.pdata[!is.na(new.pdata[[fieldName]]),]
  new.sampleNames <- row.names(new.pdata)
  es.copy <- es[rows, new.sampleNames]
  pData(es.copy) <- new.pdata
  fData(es.copy) <- data.frame(row.names=rownames(es.copy))
  es.design <- model.matrix(~0 + Comparison, data = pData(es.copy))
  colnames(es.design) <- gsub(pattern = fieldName,
                              replacement = '',
                              x = make.names(colnames(es.design)))
  fit <- lmFit(es.copy, es.design)
  fit2 <- contrasts.fit(fit, makeContrasts(B - A,
                                           levels=es.design))
  fit2 <- eBayes(fit2)
  de <- topTable(fit2, adjust.method="BH", number=Inf)
  de <- de[row.names(fData(es.copy)),]
  f <- tempfile(pattern = "de", tmpdir = getwd(), fileext = ".bin")
  writeBin(protolite::serialize_pb(as.list(de)), f)
  f
}
\end{lstlisting}

\section{Инфраструктура проекта phantasus}
\subsection{Структура git-репозитория}
Как следует из главы об архитектуре проекта, проект состоит из двух составляющих:
\begin{itemize}
\item \emph{phantasus.js} --- \emph{fork} репозитория \emph{morpheus.js}~\cite{morpheus};
\item \emph{phantasus} --- репозиторий для \emph{R}-пакета.
\end{itemize}
Внутри репозитория \emph{phantasus} находится подмодуль для репозитория \emph{phantasus.js}.
Соответственно, чтобы загрузить целиком весь проект достаточно вызвать команду из листинга~\ref{clonerepo}.
\begin{lstlisting}[float=!h,language=bash,label={clonerepo},caption={Клонирование репозитория проекта phantasus}]
  git clone --recursive https://github.com/ctlab/phantasus.git
\end{lstlisting}

\subsection{Запуск приложения}
На данный момент существует два варианта запуска приложения \emph{phantasus}. В данном разделе будут описаны оба.
\subsubsection{Единый R-пакет phantasus}
Так как проект теперь существует в виде  единого \emph{git}-репозитория, его легко можно использовать как полноценный R-пакет, содержащий в себе в том числе и файлы для веб-приложения.
С помощью функции, представленной на листинге~\ref{serve.R}, можно запускать веб-приложение phantasus непосредственно из R.
\lstinputlisting[float=!h,language=R,label={serve.R},caption={Функция для запуска приложения из R}]{listings/serve.R}

Соответственно, чтобы запустить приложение, необходимо выполнить код, представленный на листинге~\ref{launchR} и перейти на \texttt{http://localhost:8000} в браузере.

\begin{lstlisting}[float=!h,language=R,label={launchR},caption={Запуск приложения в качестве R-пакета}]
library(phantasus)
example(servePhantasus)
\end{lstlisting}

\subsubsection{Docker-образ phantasus}
На \texttt{hub.docker.com} существует автоматический репозиторий, привязанный к \emph{git}-репозиторию \emph{phantasus}. Для каждой перекомпиляции он использует Dockerfile с листинга~\ref{dockerfile}, расположенный в репозитории.

На данный момент существуют две ветки Docker-образа:
\begin{itemize}
\item \emph{master} --- компиляция происходит из \emph{master}-веток составляющих проекта, чаще всего эти скомпилированные образы стабильны и отправляются в открытый доступ для использования;
\item \emph{develop} --- компиляция происходит из \emph{develop}-веток составляющих проекта, эти образы используются для тестирования всего приложения в целом, тестирования нового функционала и не предназначены для использования на серверах.
\end{itemize}

Чтобы загрузить \emph{Docker}-образ, нужно воспользоваться командой с листинга~\ref{dockerpull}
\begin{lstlisting}[float=!h,language=bash,label={dockerpull},caption={Загрузка Docker-образа phantasus}]
  docker pull dzenkova/phantasus
\end{lstlisting}

Для запуска \emph{Docker}-контейнера, необходимо воспользоваться командой с листинге~\ref{launchdocker}.
\begin{lstlisting}[float=!h,language=bash,label={launchdocker},caption={Запуск Docker-контейнера}]
  docker run -t -d -p "800$tag:80" -v /mnt/data/phantasus-cache:/var/phantasus/cache dzenkova/phantasus
\end{lstlisting}

\subsection{Кэш для данных из GEO}
Независимо от способа запуска, данные, загруженные из \emph{GEO}, кэшируются в определенной папке (сейчас это неизменяемое значение: \texttt{/var/phantasus/cache}), чтобы не было необходимости перескачивать их заново.

\section{Настройка с помощью Apache}
В обоих из представленных вариантах запуска, приложение запускается в формате \emph{single-user}, соответственно, настройка деятельности приложения в формате \emph{multi-user} осуществляется возможностями \emph{Apache}. Конфигурационный файл можно увидеть в приложении~\ref{apacheconf}.

\subsection{Переадресация OpenCPU-сервера}
Как было сказано в разделе~\ref{opencpuintro} обзора, \emph{OpenCPU}-сервер работает быстрее в однопользовательском режиме. Сервер запускаетстся на определенном порту (например, 8001), после для корректной работы приложения, происходит переадресация запросов с \texttt{/ocpu} на \texttt{//localhost:8001/ocpu} и наоборот.

\subsection{Балансировщик для multi-user соединения}
Из-за того, что используется однопользовательский \emph{OpenCPU}-сервер, несколько людей, использующих веб-приложение \emph{phantasus} одновременно, вынуждены ждать, пока закончится запрос для одного.

Чтобы избежать такой ситуации, запускается четыре одинаковых экземпляра \emph{OpenCPU}-сервера и с помощью \emph{Apache}-балансировщика можно получать доступ к \emph{R}-серверу параллельно.

\chapterconclusion
В данной главе были рассмотрены реализованные инструменты и методы анализа:\begin{itemize}
\item загрузка данных из \emph{Gene Expression Omnibus};
\item метод главных компонент и его визуализация;
\item кластеризация методом \emph{kmeans};
\item анализ дифференциальной экспрессии. \end{itemize}

Также были описаны технические подробности реализации веб-приложения: инструкции для запуска, варианты использования и детали настройки веб-приложения на сервере.