Multivariate Data Analysis for Metabolomics Data generated by

Multivariate Data Analysis for Metabolomics Data generated by MS / NMR Spectroscopy Metabolomics Workshop Research Triangle Park / NC July 14 th-15 th H. Thiele, Bruker Daltonik, Bremen

Why do Metabolic Profiling ? Clinical Diagnostics • find metabolic markers for disease progression (e. g. cancer) • diagnose inborn errors or other diseases • study of genetic differences Toxicology • markers for drug toxicity and drug efficacy • analyze time course of toxicological response Food Science • quality control / classification of origin • health/flavor enhancement of agrochemical products MS + NMR together: • Parallel Statistics >> more confidence • Hyphenation >> ultimate characterization tool

Fundamental Issue in MVS: Dimension Reduction Bucketing • several bucketing techniques for optimum design of variables

Dynamic Peak Bucketing Scheme a 1 a 2 a 3 a 4 a 5 b 1 b 2 b 3 b 4 b 5 b 6 a 1 0 a 2 a 3 0 a 4 a 5 b 1 b 2 b 3 0 b 4 b 5 b 6 c 1 c 2 c 3 c 4 c 5 c 6 0 0 a 1 0 a 2 a 3 0 a 4 a 5 0 0 b 1 b 2 b 3 0 b 4 b 5 b 6 c 1 c 2 c 3 0 0 c 4 c 5 c 6 0 • Spectra are bucketed one by one • Bucket table gets a new column whenever a new peak occurs • Spectra not having peaks at new positions get corresponding 0

Kernel Bucketing for LC-MS Data Bucketing Parameter e. g. m/z bucket width = 1, Kernel 0. 3 Da Time bucket width = 60 s, Kernel 10 s LC-MS chromatograms of N samples 3. 2 1 10 26 02 063. 21 00. 8303 47 30 3. 62 1. 1 04. 3 00 5 3 50 8 06 616. 2 7 00 57. 3 67 70. 1. 0 08 00 m /z Intensity m /z 37 Intensity Bucketing 100 Time m/z 1 Da 100 0. 3 Da 10 s 100 30 100 Intensity Rt-m/z . . . 1 min - 1 min - pairs 535 m/z 536 m/z 537 m/z Sample X 500 50 5 Sample Y 0 25 0 Sample Z 500 50 5 . . . 50 100 60 s Time [s] Table of N samples

Which Bucketing Technique to be used ? • Rectangular, equidistant bucketing standard, good compromise if no a priori knowledge • Variable sized bucketing makes shifts ineffective, allows selective usage • Point wise bucketing often used for broad line spectra as a special case of rectangular, equidistant bucketing • Dynamic peak bucketing allows very fine bucketing without getting huge tables, requires stable shifts or masses • Kernelized bucketing variant of rectangle bucketing to reduce effect of shifts

Data Preprocessing : Spectral Background Subtraction • Measured Data are contaminated by solvents and chemical noise • Intensity of contaminants may dominate the relevant data Chemical Noise Solvent at m=75. 2, Baseline and Scaled Noise Estimate Detection of Traces by dynamic grouping

Data Preprocessing : Spectral Background Subtraction • Subtraction of spectral background makes relevant data visible Hidden traces of m=180. 2 and m=208. 2 in BPC Intensity Base Peak Chromatogram (BPC) before and after Background Subtraction time Intensity Visible traces of m=180. 2 and m=208. 2 in BPC Intensity time

Peak Picking Tasks • Find compounds defined by RT, m/z, z and area • Take together isotopic peaks and charge states

Multivariate Statistics in Metabolomics Spectroscopy How to analyze large numbers of complex LC-MS chromatograms or NMR spectra with the target of simple discrimination or grouping? - healthy / non-healthy - high / low quality NMR or MS = sensor LC-MS chromatogram or NMR spectrum= fingerprint Use Pattern Recognition Techniques!

Pattern Recognition (PR) Objectives of PR: • Statistical characterization • Model building • Classification Methods of PR: Exploratory Data Analysis • Statistical Tests • Principle Component Analysis (PCA) >> variance analysis Unsupervised Pattern Recognition • Cluster Analysis Supervised Pattern Recognition • Discriminant Analysis: LDA, PCA-DA • Classification of samples by various means e. g. Genetic Algorithm, SVM, …

Idea of Principal Component Analysis (PCA) y x y x x PC 3 PC 2 PC 1 e. g. Principal Component Analysis

Classification using PCA Input spectra list Bucketing coordinate transformation distance measures comparison to critical values classification model

Coordinate Transformation PC 2 ppm 1 ppm 2 PC 1 Scores Loadings

Baby Urine Samples : PCA - NMR mevalonic aciduria maple syrup disease PC 1, 2 orotic aciduria 200 normal candidates Pattern in PC 1/PC 2 scores plot reveals candidates with inborn errors.

New Born Screening by NMR > 400 baby urines, PCA, disease vectors indicating strength of metabolic disorders PC 1/PC 2 PC 3/PC 4 PC 11/PC 12 Bucket Analysis from 9 to 0. 4 ppm in 0. 04 ppm steps Excluded: 6 to 4. 5 ppm residual water and urea Results of BEST-NMR at 600 MHz 1 D-spectra Noesy presat 64 scans 6 c 6

New Born Screening by NMR Distance from normals distribution is a measure for concentration of the molecule representing an inborn error Hippuric acid vector CH 2 group of hippuric acid 4 c 4

PC 2 -Scores PCA : NMR vs. LC-MS PC 1 -Scores Fig. 1: Scores plot of NMR data from baby urines (born 2003). Fig. 2: Scores plot of LC-MS data from a subset of baby urines (born 2003).

LC-MS data of Samples 114 and 94 Sample 114 BPC: -All MS 8 6 4 -MS, 5. 5 -6. 8 min 263. 1037 2 264. 1062 Intensity * 10 4 0 Sample 94 BPC: -All MS 8 263. 1037 6 C 13 H 15 N 2 O 4 , 263. 10 264. 1068 4 260 261 262 263 Measured Pattern 264 Calculated Pattern 265 266 267 m/z 0 8 6 e. Xpose 94 vs. 114 BPC: -All MS 263. 1 2 4 2 0 2 4 6 8 t [min] Generate Molecular Formula of mass 263. 1037 m/z @ 5. 9 min. The determined formula C 13 H 15 N 2 O 4 corresponds to phenyl-acetylglutamate.

Combining Spectra and Statistical Data Interpretation scores / loadings Loadings in PCA indicate the importance of the original variables (buckets) in the variance space. In ideal cases a set of loadings refers to signals of a compound.

Analysis of Bucket Variables Menu bar / Options Data Viewer Loadings Plot Bucket table

Covariance matrix The covariance matrix looks like a TOCSY, cross - peaks indicate correlated fluctuations. This includes multi molecular fluctuations. Rows at cursor position are shown on top of the 2 D matrix.

Covariance Analysis row from covariance matrix reference spectrum from spectra base Interesting rows can be saved to disk as 1 D NMR spectra and used for spectra base searching as any other 1 D spectrum. Often, a small number of compounds from the spectra base match well while others do not.

PCA analysis of 69 newborn urine LC-MS spectra 1: selecting two LC-MS runs differing in the PC 1 values from scores plot 2: selecting bucket (spectral region) from loadings plot with high PC 1 value Scores plot (PC 1 -PC 2) Loadings plot (PC 1 -PC 2) Generate Sum. Formula no peak LC-MS (run 1) peak LC-MS (run 2)

Sum-Formula Generation Electron Configuration Intensity M+ M+* C/H Ratio, Elemental Limits m/z N-rule; isotope distribution double bond equiv. Isotope Masses, Abundances Experimental Peak Fast, exact Intensities & Masses calculation Molecular Constraints Fast Formula Generator using CHNO Algorithm List of Hits & Mass/Intensity Patterns of isotopic Patterns

Formula Scoring: Isotopic pattern as additional decision criteria for elemental composition List of Hits & Intensity Mass/Intensity Patterns Three independent Scores: • Intensity Ratios • Intensity weighted mean Masses • Intensity weighted Peak Distances m/z Experimental Intensities & Masses { Theoretical Intensities & Masses }

Calculating the elemental composition Simulated mass spectrum of Chlorpyriphos 12 C isotope peak 13 C isotope peak (11% int. ) 3 x Cl isotope peaks

Clinical Proteomics The samples are different The experimental and mathematical techniques are similar But the goal is the same

Workflow Clinical Proteomics Binding Washing Patients Serum Samples Isolation Normal Disease Analysis Elution Detection * * * Clinical Results MALDI-TOF MS W. Pusch et. al. , Pharmacogenomics (2003) 4(4), 463 -476

Data Preparation Aim: Extraction of the same set of features from each individual spectrum. These features will be used for model generation and later for classification of new spectra. As with metabonomics the identification of the features is of large interest Steps: • Quality Checks for spectra • Recalibration, Baseline correction, Noise Reduction • Peak detection and area calculation • Normalization of peak areas

Tasks for Clicinal Proteomics Data Analysis 1. Data preprocessing 2. Peak annotation 3. Statistical characterization 4. Discriminance analysis 5. Most of the tasks are quite similar for both kinds of applications except for the dimensionality of the original data: 6. 1 D MALDI 7. 2 D LCMS-ESI MS

Data Preprocessing : Recalibration Problem: peaks are not aligned to each other as in this example; For LC-MS it is usually the retention time Solution: application of a recalibration algorithm Result: peaks are aligned to each other

Data Preprocessing : Recalibration Intensity Mass shift Above mass tolerance Prominent peak shifted In tol. but not prominent Spec 1 Spec 2 linear mass shift m/z • Selection of prominent peaks (e. g. 30% occurrence) • Use this peak list with average masses as calibrants • Assignment of peaks and calibrants with a mass tolerance • Recalibration of all spectra by solving least square problems for a linear mass shift

Data Transformation Wavelet vs. Fourier Aim: Transformation of spectrum from • time-amplitude (mass spectra) domain • into time-frequency (Fourier) or • time-scale (Wavelet) representation Benefit: - decomposition into distinct frequency/scale bands - significant features (peaks, patterns) occur on specific frequencies/scales Steps for wavelet decomposition: low pass filter approximation coefficients high pass filter detail coefficients

Wavelets for Feature Extraction Aim: Determine features from the spectra which are discriminant for class separation Method: Wavelet-Transformation • gives information about the signal • localized in time (m/z) and frequency • lower freq. corresponds to raw structural information of the spectrum • higher freq. corresponds to fine/detailed information • in contrast to FFT we get knowledge were the feature is located in time (m/z) Feature selection: • we get much features (dep. on time resolution) • Brute force + sophisticated feature selection needed

Peak detection Problem: • Common sets of peaks needed for later model stage • Different peaks vary to a different extent over all spectra • Small peaks , nevertheless giving a good separation between classes, might be overlooked only considering single spectra average spectrum single spectrum

Peak detection Solution - Clin. Pro. Tools: • Determination of peak positions by use of Average-Spectrum • Integration over start and end Masses for detected Peaks Blue areas indicate picked peaks average spectrum Red area for picked peaks in Model (see later)

Average spectra per class Peak at ca. 2022 Da – Idx 24 in Model of most imp. 15 Peaks for GA and SVM Avg. -spec class 1 Avg. -spec class 2 Avg. -spec class 3 Avg. -spec class 4

Univariate Statistics: Getting a basic idea about the data Calculation based on: peak intensities / peak areas • descriptive / robust statistics • Welch`s t-test / Wilcoxen test Statistic peak area sorted according p-value

Algorithms for Discriminate Analysis Some alternatives to classical linear DA: • Feature selection: Genetic Algorithms (GA) + Cluster Analysis • Support Vector Machines (SVM)

GA: Application to MS data • Solution = combinations of peaks • Initial population: randomly generated solutions • Start with multiple initial populations using a migration schema • Each solution is assigned a fitness value according to its ability to separate two or more classes (using centroid or KNNclustering and by determination of between and within class distances) • New generations of population are formed using – Selection: the fitter a solution, the higher the chance for being selected as a parent – Crossover: parents form new solutions by exchanging some of their peaks, new solutions replace parents – Mutation: random changes in solutions • Result: combinations of peaks, which separate classes best

GA: Genetic Evolution Chromosome 1 Start Set Mutation Chromosome 2 1000 1200 1700 2100 2500 1700 1800 2000 2200 2300 1000 1200 1500 2100 2500 1700 1800 2000 2150 2300 1000 1200 1500 2150 2300 1700 1800 2000 2100 2500 50 -500 cycles Cross Over 1000 1200 1500 Selection 2150 Discard disadvantages 2300 1700 1800 2000 Fitness Test using k-NN 2100 2500 Keep advantages

KNN: k-nearest neighbor clustering • Spectrum = point in Rn (e. g. areas of selected peaks) • Determination of k nearest neighbors for each spectrum • Classification of all points using classes of neighboring points • Example: point A is classified as class 2, point B as class 1 • Fitness value: – percentage of correctly classified points – calculation of between/within distances A B Legend: class 1 class 2

Centroid clustering • Spectrum = point in Rn (e. g. areas of selected peaks) • Spectrum by spectrum is analyzed (iterative process): – if it is the first spectrum or too far away from all existing clusters, a new cluster with just this spectrum is created – otherwise it is assigned to the nearest cluster, the centroid is recalculated • Fitness value: – pureness & #clusters (optimal: k clusters, with all spectra of one class) – calculation of between/within distances 3 4 1 2

GA: Results Prediction capability of GA (plot for best 2 peaks) ~76% pred. * ~91% pred. * ~66% pred. * * Prediction acc. for a model with 25 peaks

SVM: Support Vector Machine • SVM: Calculation of direction in Rn, which separates best between two classes (supervised method) • PCA (principal component analysis): calculation of direction in Rn, which best explains variability (unsupervised method, i. e. without looking at class memberships) PCA = PCA SVM support vectors

What is SVM – Basic problem • Assume we have two classes of data points with two peaks • Now I look for that line which – optimal - separates these 2 classes Class 2 Many decision boundaries can separate these two classes. Which should be chosen? Class 1 The green boundaries are valid but bad ones

What is SVM – Basic idea The decision boundary should be as far away from the data of both classes as possible. We should maximize the margin, m: Class 2 Class 1 m This problem can be solved by mathematical optimization theory

SVM: Application to MS data • SVM: quadratic optimization problem, solved by an iterative process using Sequential Minimal Optimization • In simplest case a hyperplane separating classes is calculated • Therefrom contribution of individual peaks is calculated From Spectra in 3 classes : We get: • Separating hyperplanes • Recognition & Prediction accuracy • Peak ranking highlighting potential biomarker-patterns

SVM: Application to MS data - results Note its plotted in 2 D but in fact it is high dimenensional Rec. Pred. Class 1 90 69 Class 2 86 83 Class 3 90 76 Rank 1 2 3 … Index 7 17 18 … Mass 1212 1450 1469 …

SVM: Results Prediction capability of SVM (plot for best 2 peaks) ~93% pred. * ~70% pred. * ~75% pred. * * Prediction acc. for a model with 25 peaks