\documentclass[]{article}
\usepackage{lmodern}
\usepackage{amssymb,amsmath}
\usepackage{ifxetex,ifluatex}
\usepackage{fixltx2e} % provides \textsubscript
\ifnum 0\ifxetex 1\fi\ifluatex 1\fi=0 % if pdftex
  \usepackage[T1]{fontenc}
  \usepackage[utf8]{inputenc}
\else % if luatex or xelatex
  \ifxetex
    \usepackage{mathspec}
  \else
    \usepackage{fontspec}
  \fi
  \defaultfontfeatures{Ligatures=TeX,Scale=MatchLowercase}
\fi
% use upquote if available, for straight quotes in verbatim environments
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
% use microtype if available
\IfFileExists{microtype.sty}{%
\usepackage{microtype}
\UseMicrotypeSet[protrusion]{basicmath} % disable protrusion for tt fonts
}{}
\usepackage[margin=1in]{geometry}
\usepackage{hyperref}
\hypersetup{unicode=true,
            pdftitle={Recount Methylation manuscript supplemental information},
            pdfauthor={Sean Maden},
            pdfborder={0 0 0},
            breaklinks=true}
\urlstyle{same}  % don't use monospace font for urls
\usepackage{graphicx}
% grffile has become a legacy package: https://ctan.org/pkg/grffile
\IfFileExists{grffile.sty}{%
\usepackage{grffile}
}{}
\makeatletter
\def\maxwidth{\ifdim\Gin@nat@width>\linewidth\linewidth\else\Gin@nat@width\fi}
\def\maxheight{\ifdim\Gin@nat@height>\textheight\textheight\else\Gin@nat@height\fi}
\makeatother
% Scale images if necessary, so that they will not overflow the page
% margins by default, and it is still possible to overwrite the defaults
% using explicit options in \includegraphics[width, height, ...]{}
\setkeys{Gin}{width=\maxwidth,height=\maxheight,keepaspectratio}
\IfFileExists{parskip.sty}{%
\usepackage{parskip}
}{% else
\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt plus 2pt minus 1pt}
}
\setlength{\emergencystretch}{3em}  % prevent overfull lines
\providecommand{\tightlist}{%
  \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}
\setcounter{secnumdepth}{0}
% Redefines (sub)paragraphs to behave more like sections
\ifx\paragraph\undefined\else
\let\oldparagraph\paragraph
\renewcommand{\paragraph}[1]{\oldparagraph{#1}\mbox{}}
\fi
\ifx\subparagraph\undefined\else
\let\oldsubparagraph\subparagraph
\renewcommand{\subparagraph}[1]{\oldsubparagraph{#1}\mbox{}}
\fi

%%% Use protect on footnotes to avoid problems with footnotes in titles
\let\rmarkdownfootnote\footnote%
\def\footnote{\protect\rmarkdownfootnote}

%%% Change title format to be more compact
\usepackage{titling}

% Create subtitle command for use in maketitle
\providecommand{\subtitle}[1]{
  \posttitle{
    \begin{center}\large#1\end{center}
    }
}

\setlength{\droptitle}{-2em}

  \title{Recount Methylation manuscript supplemental information}
    \pretitle{\vspace{\droptitle}\centering\huge}
  \posttitle{\par}
    \author{Sean Maden}
    \preauthor{\centering\large\emph}
  \postauthor{\par}
      \predate{\centering\large\emph}
  \postdate{\par}
    \date{2/2/2020}


\begin{document}
\maketitle

\hypertarget{obtaining-and-compiling-geo-data}{%
\section{Obtaining and compiling GEO
data}\label{obtaining-and-compiling-geo-data}}

We identified DNAm array samples published to GEO and available on the
GEO Data Sets (GDS) database as of March 2019. We queried GDS, which
provides FTP access to GEO data, using the \href{}{Entrez Utilities}
software. This allowed us to obtain available samples by DNAm array
platform and graph sample and study availability by year (Figures 1 and
S1).

Upon recognition of samples with IDATs available for the HM450K
platform, we programmatically managed data downloads with versioning
using \texttt{recount-methylation-server}, a software we developed
utilizing the celery job que manager and UTP timestamps for file
versioning. Downloads were organized as tasks around GSE IDs containing
valid samples, where each task comprised of downloading available GSM
IDATs and the GSE SOFT file containing experiment and sample metadata.

We used additional scripts to scrape GSM-specific metadata from the
GSE-specific SOFT files and stored scraped data as GSM-specific JSON
files. Certain GSE-specific data was also removed from the resulting
JSON files, such as methods and experiment title, to avoid learning
study-specific metadata and increase sample-specific learning in
metadata processing downstream.

\hypertarget{processing-sample-metadata}{%
\section{Processing sample metadata}\label{processing-sample-metadata}}

We performed extensive processing of sample metadata in order to clean
and harmonize data originating from the SOFT files.

\hypertarget{learning-metadata-annotation-labels}{%
\subsection{Learning metadata annotation
labels}\label{learning-metadata-annotation-labels}}

We initially reviewed metadata encoding patterns manually on a
GSE-specific basis. Key sample-specific terms useful for research,
including terms pertaining to tissue of origin, disease context, sex,
and age, were then coerced under corresponding variables or under a
``miscellaneous'' variable. Considering this as a preprocessing step, we
then conducted postprocessing of the coerced metadata that used regular
expression string pattern matching and logic.

Postprocessing resulted in uniformly coded and cleaned metadata labels.
For instance, for the variables \texttt{tissue} and \texttt{disease},
labels are mainly lower-case and separated by semicolons. Term labels
are then discoverable and queryable with regular expression pattern
mapping. We conducted additional annotation of sample storage practices
(e.g.~frozen or formalin fixed paraffin-embedded) by mining of the
``miscellaneous'' term from preprocessing and manual review of the
largest available studies.

\hypertarget{metasra-pipeline}{%
\subsection{MetaSRA-pipeline}\label{metasra-pipeline}}

We used the MetaSRA-pipeline to learn sample type predictions from
GSM-specific JSON files mined from SOFT files. This pipeline uses
machine learning to map metadata to ENCODE ontology terms and form a
probability for each of several possible sample types. We retained the
highest-probability sample types under a \texttt{sampletype} variable in
the postprocessed metadata, along with the returned probability.

\hypertarget{model-based-metadata-predictions}{%
\subsection{Model-based metadata
predictions}\label{model-based-metadata-predictions}}

We performed model-based predictions of sex, age, and cell type
fractions using noob-normalized DNAm (Beta-values) compiled from
available GSM IDATs. Sex and cell fractions were predicted using methods
from the \texttt{minfi} package, and age predictions were performed
using the clock model implemented in the \texttt{wateRmelon} package and
first presented in Horvath 2013. Prior to calculating predictions for
sex and cell fractions, light preprocessing is performed that may result
in prediction variation depending on the number and type of samples
included. When we tested prediction concordance from randomizing samples
to calculations, we found \textgreater95\% concordance across
randomizations.

\hypertarget{deriving-dnam-array-data}{%
\section{Deriving DNAm array data}\label{deriving-dnam-array-data}}

The HM450K Illumina BeadArray DNAm array platform utilizes two probes
(Type I and Type II) to probe over 480,000 CpG loci genome-wide. It also
includes specially designed control probes as well as probes targeting
high-frequency SNPs for certain types of genetic inference. Owing to the
amount of samples and size of data files here, we opted to compile
fairly large (\textasciitilde120 Gb) tables of various vital DNAm and
related data, and much of our analyses were implemented in parallel
across rows of these large data files.

\hypertarget{dnam-signal-and-beta-values}{%
\subsection{DNAm signal and
Beta-values}\label{dnam-signal-and-beta-values}}

From available IDATs, we used \texttt{minfi} to extract and store signal
tables for each channel color (red and green), and we used these to
derive raw Beta-values. The Beta-value reflects the fraction of
methylated DNA strands in a sample, and is thus constrained from 0 - 1
(see Methods). We further performed within-sample normalization using
the noob method in \texttt{minfi} to produce normalized Beta-values and
extact methylated and unmethylated signal tables.

\hypertarget{deriving-approximate-genome-wide-and-autosome-wide-distance-matrices}{%
\subsection{Deriving approximate genome-wide and autosome-wide distance
matrices}\label{deriving-approximate-genome-wide-and-autosome-wide-distance-matrices}}

Aging owing to the amount of data involved, we needed to perform an
intermediate dimensionality-reduction step prior to principal component
analyses (PCAs). We used feature hashing, or the hashing trick, to
preserve between-probe distances while sufficiently reducing data
dimensionality so that PCA could be run in active memory (see below). We
performed feature hashing on raw and normalized Beta-values for
genome-wide and autosomal probes, respectively.

\hypertarget{deriving-control-metric-signals}{%
\subsection{Deriving control metric
signals}\label{deriving-control-metric-signals}}

Consulting published literature and recommendations, we wrote code to
calculate 17 BeadArray control metrics from the HM450K control probe
signal from sample IDATs (see values in Table S2, and BeadArray metric
descriptions in Table S3). We especially based this work off of
functions in the \texttt{ewastools} package that also derives the
BeadArray metrics, we few minor discrepancies where we felt our code was
an improvement. This included inclusion of a denominator offset of 1
where certain metrics would otherwise divide by 0, and inclusion of
certain control probe type that were omitted prior. We note these
changes impacted a few of the metrics and may have slightly impacted
sample signal distribution relations to published metric thresholds.

We further derived methylated and unmethylated log2 medians from noob
Beta-values. Examination of samples in the 2d plane formed by the
metrics has been used to identify likely low-quality or unreliable
samples prior.

\hypertarget{preprocessing-dnam-across-7-distinct-tissues}{%
\subsection{Preprocessing DNAm across 7 distinct
tissues}\label{preprocessing-dnam-across-7-distinct-tissues}}

We implemented additional filters and preprocessing for a subset of
samples from one of 7 types of tissues. We initially removed all samples
containing likely from cancer tissue or cancer patients. We then
reviewed available studies linked to remaining GSE IDs, noting the
disease studied and included sample types in Table S6. We also removed
samples showing study-specific low signal for both methylated and
unmethylated log2 medians (\textless{} 5th quantile) for studies with
\textgreater4 samples. Next, we filtered likely replicates using the
probabilitic model implemented in \texttt{ewastools}, which utilizes the
array probes mapping to high-frequency SNPs. Finally, we only retained
tissues with \textgreater99 remaing samples from at least 2 studies.

For our cross-study analyses, we reasoned that cross-study batch effects
would have greater confounding impact than within-study batch effects.
After preliminarily performing linear model adjustment of DNAm across
samples from all 7 tissues, we observed most probes showed drastic
``polarity swapping'' after adjustment, which indicates failed
adjustment due far greater study-specific than tissue-specific impact on
variances. For this reason, we decided to perform adjustments within
tissues before further within-tissue preprocessing.

After linear adjustment on GSE ID, we removed sex chromosome probes.
Within each tissue, we then performed array-wide probe-specific ANOVAs
from multiple regressions adjusting on GSE ID, predicted sex, predicted
age, and predicted cell type fractions. Probes we significant and
high-magnitude variance contributions from any of these variables were
the removed.

\hypertarget{analyses-of-dnam-data}{%
\section{Analyses of DNAm data}\label{analyses-of-dnam-data}}

\hypertarget{geo-year-wise-sample-and-study-availability}{%
\subsection{GEO year-wise sample and study
availability}\label{geo-year-wise-sample-and-study-availability}}

Using the data mined with FTP calls to GDS using E-utilities, we plotted
year-wise cumulative abundances of GSMs and GSEs using the
\texttt{ggplot2} package.

\hypertarget{metadata-concordance-tests-and-analyses}{%
\subsection{Metadata concordance tests and
analyses}\label{metadata-concordance-tests-and-analyses}}

We evaluated concordance by using base R functions to calculate
fractional concordances for sex, and linear modeling, R-squared, and
Spearman's Rho for age. Scatterplots showing age concordance (Figure S3)
were generated using base R.

\hypertarget{analyzing-control-metric-signal}{%
\subsection{Analyzing control metric
signal}\label{analyzing-control-metric-signal}}

We assessed control metric signal (BeadArray metrics and methylated or
unmethylated log2 medians) using base R and certain packages. We
calculated metadata term enrichment for the \texttt{tissue} and
\texttt{disease} annotations using the Bimodal test with multiple
testing adjustment using the Benjamini-Hotchberg method (p-adjusted
\textless{} 1e-3 significance).

Stacked barplots, confidence intervals, and violin plots in Figure 2
were generated using \texttt{ggplot}. Study-wise sub-threshold frequency
heatmaps were generated using \texttt{ComplexHeatmap}.

\hypertarget{analyzing-pca-across-samples-and-probes}{%
\subsection{Analyzing PCA across samples and
probes}\label{analyzing-pca-across-samples-and-probes}}

We performed PCA using the feature-hashed probe datasets and the
\texttt{prcomp} function with default arguments (see above). Plots were
generated using the \texttt{factoextra} package and base R.

\hypertarget{analyzing-hypo--and-hypervariable-probes-across-7-tissues}{%
\subsection{Analyzing hypo- and hypervariable probes across 7
tissues}\label{analyzing-hypo--and-hypervariable-probes-across-7-tissues}}

For analysis, we used several approaches to identify hypo- and
hypervariant CpGs within and across tissues. Due to preprocessing within
individual tissues, our analyses were somewhat limited to set-based
characterizations rather than direct differential variance tests.


\end{document}
