A Novel Strategy for Large-Scale Metabolomics ... - ACS Publications

Subscriber access provided by La Trobe University Library

Article

A Novel Strategy for Large-Scale Metabolomics Study by Calibrating Gross and Systematic Errors in Gas Chromatography-Mass Spectrometry Yanni Zhao, Zhiqiang Hao, Chunxia Zhao, Jieyu Zhao, Junjie Zhang, Yanli Li, Lili Li, Xin Huang, Xiaohui Lin, Zhongda Zeng, Xin Lu, and Guowang Xu Anal. Chem., Just Accepted Manuscript • DOI: 10.1021/acs.analchem.5b03912 • Publication Date (Web): 12 Jan 2016 Downloaded from http://pubs.acs.org on January 18, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Analytical Chemistry is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Analytical Chemistry

A Novel Strategy for Large-Scale Metabolomics Study by Calibrating Gross and Systematic Errors in Gas Chromatography-Mass Spectrometry

Yanni Zhao1, Zhiqiang Hao2, Chunxia Zhao1, Jieyu Zhao1, Junjie Zhang1, Yanli Li1, Lili Li1, Xin Huang2, Xiaohui Lin2, Zhongda Zeng1，Xin Lu1*, Guowang Xu1*

1

Key Laboratory of Separation Science for Analytical Chemistry, Dalian Institute of

Chemical Physics, Chinese Academy of Sciences, Dalian 116023, China. 2

School of Computer Science & Technology, Dalian University of Technology,

Dalian 116023, China.

* Address correspondence to: Prof. Dr. Guowang Xu, Key Laboratory of Separation Science for Analytical Chemistry, Dalian Institute of Chemical Physics, Chinese Academy of Sciences. Phone / Fax: +86-411-84379530. E-mail: [email protected]. Prof. Dr. Xin Lu, Key Laboratory of Separation Science for Analytical Chemistry, Dalian Institute of Chemical Physics, Chinese Academy of Sciences. Phone / Fax: +86-411-84379559. E-mail: [email protected].

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Page 2 of 29

ABSTRACT:

2

Metabolomics is increasingly applied to discover and validate metabolite

3

biomarkers and illuminate biological variations. Combination of multiple analytical

4

batches in large-scale and long-term metabolomics is commonly utilized to generate

5

robust metabolomics data, but the gross and systematic errors are often observed. The

6

appropriate calibration methods are required before statistical analyses. Here, we

7

develop a novel correction strategy for large-scale and long-term metabolomics study,

8

which could integrate metabolomics data from multiple batches and different

9

instruments by calibrating gross and systematic errors. Gross error calibration method

10

applied various statistical and fitting models of the feature ratios between two

11

adjacent quality control (QC) samples to screen and calibrate outlier variables. Virtual

12

QC of each sample was produced by a linear fitting model of the feature intensities

13

between two neighboring QCs to obtain a correction factor and remove the systematic

14

bias. The suggested method was applied to handle metabolic profiling data of 1197

15

plant samples in nine batches analyzed by two gas chromatography-mass

16

spectrometry instruments. The method was evaluated by the relative standard

17

deviations of all the detected peaks, the average Pearson correlation coefficients and

18

Euclidean distance of QCs and non-QC replicates. The results showed the established

19

approach outperforms the commonly used internal standard correction and total

20

intensity signal correction methods, it could be used to integrate the metabolomics

21

data from multiple analytical batches and instruments, and allows the frequency of

22

QC to one injection of every 20 real samples. The suggested method makes a large

23

amount of metabolomics analysis practicable.

24

Keywords: metabolomics, large-scale, gross error, systematic bias

25 2


Page 3 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


26

INTRODUCTION

27

Metabolomics is an indispensable branch of systems biology research which is

28

focused on the holistically quantitative and qualitative investigations of the changes of

29

the low-molecular-weight metabolites (≤ 1500 Da) related to genetic, environmental

30

and disease perturbations in biological samples.1-3

31

Chromatography coupled with mass spectrometry (MS) and nuclear magnetic

32

resonance (NMR) spectroscopy are the most commonly applied tools in metabolomics

33

study.4-7 Among them, gas chromatography-mass spectrometry (GC-MS) could

34

efficiently separate and detect some of hydrophilic endogenous metabolites (e.g.,

35

sugar, amino acids and organic acids) after derivatization,8 which are essential

36

substances to maintain the life activities.9,10 To identify and validate the physiological

37

perturbations and metabolite biomarkers of different living systems in the primary

38

metabolome, large amount of samples are to be analyzed using a high-throughput

39

approach based on GC-MS platform.11 In our previous work, a pseudotargeted

40

metabolomics approach has been developed in which the characteristic ion

41

information obtained from the nontargeted metabolomics data without the aid of

42

standard references was used to acquire the metabolome data by using the selective

43

ion monitoring (SIM) GC-MS. It displayed wider linear range, higher sensitivity and

44

better repeatability than nontargeted GC-MS in the full-scan mode.10,12 A stability and

45

repeatability platform is required for the analysis of large sample sizes, unfortunately,

46

until now, it is still not easy to analyze all samples in an analytical batch due to MS

47

signal drift in long-term metabolomics studies. It is often necessary to split large-scale

48

samples into multiple batches.13 The occurring trouble is integration of these data

49

because the gross error and systematic bias are often met in metabolomics data of

50

multiple batches and they were caused by instrument cleaning and calibration, 3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

51

Page 4 of 29

accessories changing, etc.

52

In the previous studies,13 both the internal standard (IS) correction and total

53

intensity signal correction (TISC) were considered as sample-based signal calibration

54

procedures, which utilized a correction factor (e.g., IS and TISC) of each sample to

55

normalize all of the analytes in a given sample.10,14 Ideally, the IS or TISC should

56

display similar variation tendency to all of the detected metabolites for the removal of

57

systematic drifts. Actually, for cases in which hundreds to thousands of analytes are

58

detected, it is impossible that TISC or IS shows similar changes in signal patterns to

59

all the metabolic features. To solve this problem, some novel tactics of feature-based

60

signal correction have been suggested to calibrate systemic errors in metabolomics

61

study.13,15,16 Wang et al. utilized a single value regression model with the total

62

abundance information of each sample to develop a Batch Normalizer method, which

63

could remove both batch and injection order effects.16 Kamleh et al.13 employed the

64

local or global mean intensity of QCs (LoMec and GoMec) and the local or global

65

linear regression of QC intensity values (LoReg and GoRec) to formulate a correction

66

factor of each feature and further remove the systematic bias. These suggested

67

protocols showed a remarkable improvement in repeatability compared to TISC and

68

IS methods. Fages et al.15 proposed a grouped-batch profile (GBP) calibration method

69

using the systematic bias of feature in each batch to calibrate the actual metabolic

70

response. This established procedure could adjust systematic variation of NMR

71

metabolomic data sets and improve the prediction of multivariate regression models

72

for robust classification of biological samples. However, all of these methods were

73

concentrated on the systematic bias calibration, few were on gross errors correction.

74

Actually, both of gross and systematic errors are to be calibrated in order to obtain a

75

high quality metabolomics data. 4


Page 5 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


76

In this study, a GC-MS-based pseudotargeted metabolomics approach was

77

applied to acquire the metabolic profiling data of plant leaves. A novel calibration

78

strategy was developed for the integration of large-scale and long-term metabolomics

79

data, in which both systematic and gross errors could be corrected. A metabolomics

80

data set with 1197 plant leaf samples including 130 QCs in nine batches analyzed by

81

using two different instruments were applied to evaluate this suggested calibration

82

method by calculating the percentages of features (POF) with relative standard

83

deviation (RSD) less than 15%, the average Pearson correlation coefficients (r) and

84

the Euclidean distance (ED) of the replicates in principal components analysis (PCA)

85

scores plots. Additionally, the frequency of QC injection was also investigated.

86 87

EXPERIMENTAL SECTION

88

Strategy description

89

In large-scale GC-MS metabolomics analysis, some necessary instrumental

90

maintenances (e.g., changing accessories and columns, cleaning ion source and tuning

91

instrument), are frequently performed to minimize the signal intensity drifts caused by

92

instrument contamination and column age changes. Additionally, several batches are

93

analyzed due to the different sample collection times or repeatability requirement of

94

the instrument. Sometimes data obtained from different GC-MS instruments within a

95

laboratory or inter-laboratory were to be integrated or compared. In considering above

96

these cases, we design an experiment procedure to analyze 1197 plant leaves

97

including 130 QCs in nine batches (Figure 1). The samples in the first four batches

98

were consecutive injection and some changes in instrument accessories (e.g., the

99

injection port liner and silicone septum) between two batches were conducted to avoid

100

the signal decrease. After a long time intervals between batch 4 and batch 5, batch 6 5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 29

101

and batch 7, and batch 8 and batch 9, the MS detector was tuned and the most suitable

102

instrument parameters (e.g., detector voltage, RF Gain, RF offset and lens voltage)

103

were renewed and used to analyze samples. The samples in the first seven batches and

104

last two batches were acquired by GC-QP2010 Plus and GC-MS QP2010,

105

respectively. 188, 209, 25, 25, 156, 97, 354, 84, and 59 samples were analyzed in

106

batches 1, 2, 3, 4, 5, 6, 7, 8 and 9, respectively. The same QC samples and five kinds

107

of non-QC replicates (i.e., S1, S2, S3, S4 and S5) were analyzed in each batch to

108

evaluate the data quality and calibration performance. S1, S2, S3, S4 and S5 consisted

109

of 43, 43, 31, 24 and 36 fresh tobacco leaf samples, respectively.

110

Systematic and gross errors calibration methods were developed using the data of

111

tobacco leaves from first four batches. Systematic error calibration method was

112

established by inserting the virtual QCs to the linear fitting model of feature intensity

113

of two neighboring QCs to correct the each feature in every sample. The gross error

114

correction was conducted by calculating the ratio of each variable between adjacent

115

QCs to screen outliers and establishing a linear fitting model to correct the outliers.

116

Additionally, the developed method was further applied to correct data of all the 1197

117

samples, and the calibration performance was evaluated by investigating the RSD

118

values of peak intensity, the average ED and r of the samples (Figure S-1 in the

119

Supporting Information).

120 121

Experimental

122

Chemicals. Acetonitrile and isopropanol were HPLC-grade from Merck

123

(Darmstadt, Germany). Ultrapure water was obtained by filtration of distilled water

124

using a Milli-Q system (Boston, USA). Pyridine, methoxyamine hydrochloride and

125

MSTFA

(N-methyl-N-(trimethylsilyl)-trifluoroacetamide)

used

as

GC-MS 6


Page 7 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


126

derivatization reagents were purchased from Sigma-Aldrich (St. Louis, USA). Internal

127

standards (ISs), tridecanoic acid and vanilic acid, from Sigma-Aldrich (St. Louis,

128

USA) were used to monitor repeatability of sample pretreatment and instrument.

129

Damascenone, alpha-ionone, methyl palmitate, methyl linoleate and methyl oleate

130

used as external standards (ESs) to monitor stability of the instrument, were obtained

131

from Alfa Aesar (Ward Hill, USA).

132

Sample Preparation. A total of 1197 fresh tobacco leaf samples including 130

133

QC samples were acquired from three main growing districts in China, including

134

Henan, Yunnan and Guizhou provinces. All of fresh leaves were ground to

135

homogeneous powders and stored at -80 oC until metabolite extraction. The equal

136

amount of analytical samples was blended for preparation of the QC to monitor the

137

analytical method and correct the large-scale GC−MS metabolomics data.

138

The preparation workflow of the leaf samples was similar to our previous work

139

with a few modifications.10 About 10 mg leaf powder was soaked in 1.5 mL of the

140

isopropanol/acetonitrile/water (3:3:2, v/v/v) mixed solution (including 20 μL of 0.1

141

mg/mL internal standards and external standards solutions) and adequately

142

vortex-mixed for 4 min to extract leaf metabolites using a Multi-Tube Vortexer

143

(VX-2500). After centrifugation at 14 000 rpm for 10 min at 4 °C, 500 μL of

144

supernatant was collected in a new Eppendorf tube, and then freeze-dried using a

145

Labconco centrifugal concentrator for six hours. The dried residue was immersed in

146

100 μL of 20 mg/mL methoxyamine pyridine solution in 37 oC water bath for 90

147

min to protect the ketones and aldehydes. Subsequently, a 60 min silylation reaction

148

was performed by adding 80 μL of MSTFA into the above solution in 37 oC water

149

bath to increase the volatility of the metabolites. Ultimately, the derivatized sample

150

was transferred to a glass vial with a conical insert for GC-MS analysis. 7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 29

151

GC-MS Pseudotargeted Method. The samples were analyzed by two different

152

GC-MS instruments (QP 2010 and QP 2010 Plus, Shimadzu, Japan) using two

153

different DB-5 MS fused silica capillary columns (30 m × 0.25 mm × 0.25 µm,

154

Agilent, USA) in the constant flow mode with 1.2 mL/min high-purity helium

155

(99.9995%) as the carrier gas. The oven temperature was set to 70 oC, maintained for

156

3 min and then increased at 5 oC/min to 310 oC, held for 5 min. The split ratio was set

157

as 10:1 and the interface temperature was adjusted to 280 oC. The electron impact (EI)

158

at 70 eV was used as the ionization mode. The temperature of inlet and ion source

159

were 300 oC and 240 oC, respectively. The scope of mass signal acquisition under

160

full-scan mode was set to 33-600 m/z with a scan rate of 5 scans/s.

161

To establish a pseudotargeted GC-MS method,10,12 the raw data at full-scan mode

162

of a double concentration QC sample was deconvoluted for peak identification using

163

the AMDIS from the NIST (Gaithersburg, USA) and LECO ChromaTOF (St. Joseph,

164

USA) software programs. Subsequently, the selection of the characteristic ions of the

165

metabolic features was performed by calculating the ions characteristic values, which

166

can reflect the ion uniqueness among the co-eluted components,12 using self-designed

167

software and commercial programs (e.g., AMDIS and ChromaTOF). A metabolic

168

feature is defined as a detected chromatographic peak and a single metabolite can

169

have multiple metabolic features.11 The SIM scan group information was obtained by

170

setting the threshold value of the chromatographic retention time (RT) differences

171

between two adjacent peaks as 0.1 min. According to the group information, RT and

172

characteristic ions of metabolites, the SIM acquired method and quantitative approach

173

of the pseudotargeted GC-SIM-MS were established to analyze the tobacco leaves and

174

obtain the peak integration table for the further study. The instrument parameters of

175

pseudotargeted GC-MS were the same as those of full scan method, except for the MS 8


Page 9 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


176

acquisition mode. Ultimately, 319 metabolites divided into 50 groups were detected

177

by the established pseudotargeted metabolomics method.

178 179

Performance Evaluation

180

The performance of different correction methods could be evaluated using PCA,16

181

which has frequently been applied to extract the main information from

182

high-dimensional data to assess the similarity of analytical samples and the

183

repeatability of analytical method. The parameters R2X value of PCA score plot could

184

reflect the explained fractions of the total variations. To avoid suffering the

185

information loss in the data set, the R2X value of PCA analysis dominated by the

186

number of principal components (PCs), was set more than 85%. The corresponding

187

PCs were applied to calculate the average ED and r between detected replicates.

188

Additionally, the RSD values of the duplicates were another important criterion to

189

evaluate the calibration results.13,16 According to the United States Food and Drug

190

Administration (FDA) guidance,17 the threshold of RSD for the analytical features

191

was set to 15%.

192 193

RESULTS AND DISCUSSION

194

Although pseudotargeted metabolomics method has displayed better repeatability

195

compared to nontargeted method,12 its implementation in the integration of large-scale

196

metabolomics data still faces significant challenges due to the gross or systematic

197

errors derived from MS and/or chromatographic column contaminations and changes

198

of instrument accessories and environmental conditions. In order to solve these

199

problems, we developed different algorithms to calibrate the gross error and

200

systematic bias in different batches and instruments (Figure 1). 9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 29

201 202

Systematic Error Correction

203

Systematic bias is one of the primary limited factors for the integration of large

204

scale metabolomics data. In this study, we established a Virtual-QC method applying

205

a regression model of metabolic abundances of each feature in ω adjacent QC samples

206

(ω: the number of QCs) to remove the systematic variance (Figure S-2A in the

207

Supporting Information). The intensity value Axij of each metabolic feature j in a

208

special sample i was corrected by a correction factor AQCv_ij, which was obtained by

209

the linear a fitting of signal intensities of feature j in ω adjacent QC samples. Each

210

sample i will have a special AQCv_ij, which was closely related with its inserted position

211

(i+1) in a certain fitting line as defined by Eq. 1,

212

= a (i +1) + b

AQC

v _i j

(1)

213

where a and b are the slope and intercept of the regression line, respectively (Figure

214

S-2B in the Supporting Information). The new corrected intensity value Ax'ij was

215

calculated according to Eq. 2.

216

Ax ′ = ij

Ax

ij

AQC

(2)

v _i j

217 218

Virtual-QC method was compared with the commonly used calibration methods,

219

including four kinds of sample-based correction methods (e.g., TISC, TIQSC, IS and

220

ES) and five kinds of feature-based signal correction methods (e.g., LoMec, LoReg,

221

and GoRec etc.)13 (Table S-1). The average ED, r and RSD values of each peak across

222

all QC samples and the non-QC replicates were used to evaluate the calibration

223

performances.

224

It is important to note that the calibration performances of Virtual-QC method 10


Page 11 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


225

may be affected by ω values of QCs, which could be optimized according to the data

226

set characteristics of the analytical samples. For our metabolomics data, the values ω

227

from 2 to 4 displayed the similar calibration performances for QCs and non-QC

228

replicates. Hence, the Virtual-QC method at ω = 2 was adopted in the following

229

comparison.

230

For QC samples, all of systematic error correction (SEC) methods displayed

231

distinct improvements in RSD, ED and r compared to the raw data (Figure 2A). The

232

RSD values of 48.90% of the metabolites were lower than 15% for the IS1 (tridecylic

233

acid) calibration method, which outperformed the other sample-based signal

234

correction methods. Nevertheless, three feature-based signal correction procedures

235

including LoMec, LoReg and Virtual-QC methods, have the shortest ED, greatest r

236

and the strongest repeatability overall with 88.09% POF falling below the 15% RSD

237

threshold (Figure 2A). Encouraged by these results, a further optimized work of these

238

three feature-based signal correction methods was performed by calculating the

239

RSD-value, ED and r of five types of non-QC duplicates (i.e., S1, S2, S3, S4 and S5)

240

in first four batches. One interesting observation is that the Virtual-QC approach

241

displayed shorter ED, higher r and POF with RSD less than 15% compared to LoReg

242

and LoMec methods (Figure 2B). Hence, in the following procedures to study the

243

gross error calibration methods, the Virtual-QC approach was employed as a

244

systematic bias correction method.

245 246

Gross Error Correction

247

Gross error is also an important influencing factor of the correction performance

248

and caused by measurement biases, malfunctioning instruments, and so forth.18,19

249

Because the intensity value of a metabolic feature Axij was corrected by its 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 29

250

neighboring virtual AQCv_ij based on the Virtual-QC method (Figure S-2 in the

251

Supporting Information), the gross error in actual QCs could be transferred to the

252

inserted virtual AQCv_ij, further influenced the other test samples and expanded the

253

scope of the error distribution. Therefore, the gross error correction needs to be

254

conducted before systematic bias calibration. Two vital procedures including

255

screening and calibrating outliers, were optimized by comparing the correction

256

performances across the various algorithms. In the following optimization procedures,

257

the calibration performances of gross error correction method was evaluated by a

258

comparison between the gross error plus systematic bias calibration method (G+S)

259

and only SEC approach based on Virtual-QC.

260

Screening outlier features in QCs. Ideally, the intensity ratio of feature j in two

261

adjacent QCs was equal to 1 with consideration of the same QCs analyzed in our

262

experiment. Nevertheless, some feature ratios between two neighboring QCs could

263

obviously deviate from 1, which probably mainly be induced by the gross errors.

264

Hence, we firstly calculated the feature ratios between two adjacent QCs, then four

265

different methods (i.e., Each 5%, All-5%, Each-box and All-box) were tested to define

266

the potential outlier features (Table S-2 in the Supporting Information). Each-5%

267

method calculated 5% of the feature ratios in two adjacent QCs as outliers, other 95%

268

metabolites were thought to be highly reliable (Figure S-3A in the Supporting

269

Information). All-5% method applied all the ratios to calculate the number of outliers

270

based on the 95% metabolites with statistic confidence (Figure S-3B in the Supporting

271

Information). During screening process, the ratios of metabolic features were firstly

272

sorted from small to large, and then the number of outliers was equally selected from

273

two sides of ratio > 1 and < 1. Boxplot is frequently applied to find outliers in data

274

through their interquartile range (IQR), the first quartile (Q1) and third quartile (Q3).20 12


Page 13 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


275

IQR could be calculated by the Q1 subtracted from Q3. Outliers are observations that

276

fall below Q1 - 1.5(IQR) or above Q3 + 1.5(IQR). Each group ratio and all ratios

277

were used to construct Each-box (Figure S-3C in the Supporting Information) and

278

All-box (Figure S-3D in the Supporting Information) models and detect outlier

279

features, respectively.

280

For QC samples, one interesting observation is that all of the four models of

281

screening potential outlier features showed a remarkable improvement in the

282

repeatability of QCs with more than 99% POF with RSDs falling within 15% (Figure

283

3A). Additionally, the ED-values of QCs in PCA score plots of theses four algorithms

284

combined with the SEC method were also shorter than only SEC method. Further

285

studies found that Each-5% method displayed the best POF with RSD less than 15%,

286

r and the shortest ED than the other methods for five types of non-QC replicates

287

(Figure 3B). Hence, we select Each-5% method to screen potential outlier features.

288

It should be noted that in Each-5% method 5% of the feature ratios in two

289

adjacent QCs were defined as possible outliers, but some of them perhaps have a good

290

repeatability and are not necessary to be corrected. Hence, we need further confirm

291

whether they are the outlier features which are to be calibrated.

292

Firstly, a high quality QC data without outlier variance was needed to further

293

confirm outlier features and calibrate the other QCs. In this experiment, a series of QC

294

samples in a batch are named as QC0, QC1, …, QCn-1, QCn. Based on the practical

295

operation experience, greater gross error possibly is observed in the first QC sample

296

of a batch (i.e. QC0) compared with those in other QCs. Consequently, the second QC

297

(i.e., QC1) in every batch was firstly corrected to obtain a good quality QC by two

298

ratios of AQC /AQC and AQC /AQC (Table S-3 in the Supporting Information). If the ratio

299

of feature j in AQC /AQC was outside the threshold value, but was inside the threshold

2

1

3

2j

2

1j

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

300

Page 14 of 29

range for the AQC /AQC , we speculated variable j was an outlier for QC1. 3j

2j

301

Two different algorithms (i.e., To-QCn and Former-QCn) were developed to

302

further judge the outlier features in the other QCs (Table S-4 in the Supporting

303

Information). To-QCn method applied the ratio of feature j in two adjacent QCs to

304

confirm the outliers. If the ratio of a feature j in AQC /AQC' was outside the threshold

305

value, the variable j was an outlier in QCn. The Former-QCn method applied two

306

group ratios of QCs (AQC /AQC and AQC /AQC ) to further select the outlier points.

307

If the ratio of a feature j in AQC /AQC was outside the threshold value, but was inside

308

the threshold range for the AQC /AQC , we speculated variable j was an outlier for

309

QCn. The calibration performances of these two methods on the QCs and other five

310

kinds of non-QC replicates were comprehensively compared. For QCs, there were

311

100%, 98.75% and 88.09% peaks with RSDs falling within 15% in To-QCn,

312

Former-QCn and Virtual-QC, respectively. A shortest ED and highest r in the PCA

313

scores were noticed using To-QCn method to adjust the QCs (Figure 3C). Figure 3D

314

displayed that To-QCn method also produced the highest r, shortest ED and the

315

strongest repeatability profile for the five kinds of non-QC replicates. Therefore,

316

To-QCn method was considered as the best performance method for further

317

confirming outlier features, which could result in remarkable repeatability

318

improvements of the same samples over other methods.

nj

(n+1) j

(n+2) j

nj

(n+1) j

(n+2) j

(n-1)j

(n+1) j

nj

(n+1) j

319

Correction of outliers. After the outlier features of QCs were defined, we need

320

to establish a model to correct the outlier to a reasonable value. The calibration

321

equations for QC0, QC1 and other QCn (n = 2,3,4……) are shown in Eq. 3, Eq. 4 and

322

Eq. 5, respectively.

323

AQC′

0j

 AQC = AQC′ ×  0j 1j  AQC′  1j

   cor r .

(3)

14


Page 15 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


324

where AQC , AQC' , AQC , and AQC' are the unadjusted QC0, adjusted QC0, unadjusted

325

QC1 and adjusted QC1 responses, respectively.

0j

0j

1j

1j

 AQC AQC′ = AQC ×  1j 1j 2j  AQC  2j

326

   cor r .

(4)

327

where AQC , AQC' and AQC are the unadjusted QC1, adjusted QC1 and unadjusted QC2

328

responses, respectively.

1j

1j

AQC′nj = AQC′

329

( n - 1) j

2j

 A QCnj ×  AQC′  (n - 1) j

   cor r .

(n = 2, 3, 4, …)

330

where AQC , AQC' , AQC

331

uncorrected QC(n-1) and corrected QC(n-1) responses, respectively.

n j

(n-1) j

(n-1) j

are the uncorrected QCn, corrected QCn,

To obtain the calibration parameters (i.e., [AQC /AQC' ]corr., [AQC /AQC ]corr. and [AQC

332 333

n j

and AQC'

(5)

0j

j

/AQC'

(n-1) j

1j

1j

2j

n

]corr.) in Eq. 3, Eq. 4 and Eq. 5, we tried eight methods (Table S-5 in the

334

Supporting Information) based on the ratios of the features between two adjacent QCs

335

without the outlier features. C1, C2 and C3 methods are based on a linear fitting,

336

polynomial fitting of three times and six times, respectively. Additionally, correlation

337

coefficients of polynomial fitting (R2) were used as index to study different methods,

338

among them, the thresholds of R2 for C4, C5 and C6 methods were 0.9, 0.95, and 0.99,

339

respectively. C7 method directly assigned the lowest and largest threshold values to

340

the outliers from the ratios < 1 and > 1, respectively. C8 was established based on

341

principle of Virtual-QC. If the peak j of the QCn was outlier, the calibration value AQC'

342

was corrected by the virtual AQC' obtained based on AQC v_nj

(n-1)j

nj

and AQC . (n+1)j

343

The performance of the above eight methods for correcting outliers was

344

compared by handling metabolomics data of the QCs and non-QC replicates samples

345

(S1~S5). Both of C1 and C8 calibrated the QC samples to shorter average ED in the

346

PCA score plots compared to other correction methods (Figure 3E). In addition, for 15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 29

347

most of replicates, the C1 method displayed greater improvement in repeatability

348

compared to C8 (Figure 3F). Comprehensively considering the correction

349

performance of QCs and the S1~S5 samples, the C1 method was selected as the

350

optimal correction outliers method. Detailed gross and systematic errors calibration

351

algorithms were described in the Supporting Information.

352

PCA score plots were prepared to visualize the calibration performance of the

353

repeatability and similarity of the QC samples. The QC samples of each batch before

354

calibration are separated into four parts along the first two PCs of the PCA score plot

355

depending on batches despite their intrinsic similarity (Figure 4A). On the contrary,

356

the QC samples calibrated by G+S calibration method had a satisfactory performance,

357

which were randomly distributed in the PCA score plot, showing the gross error and

358

batch systematic biases were evidently reduced (Figure 4B). Furthermore, the

359

cumulative RSD curve of QCs peaks indicated that the calibration method based on

360

G+S could significantly decrease RSD value of metabolic features in QC samples and

361

result in a remarkably improved repeatability compared to the crude data (Figure 4C).

362 363

Frequency of QC Injection

364

The frequencies of QC injections could alter the correction factor AQC value,

365

which further influences the performance of each calibration method. In our study, the

366

frequencies of QC injections (F) in analytic sequence were set to one in every 5

367

analyzed samples (1/5), 1/10, and 1/20. The metabolomics data obtained by IS1, TISC,

368

SEC, the optimal G+S method and raw data (RD) without correction, were used to

369

define the optimum frequency of QC injections (Figure 4D).

v_ij

370

Compared to other methods, G+S method showed the strongest improvement for

371

repeatability with RSDs of all features less than 15% under different F (i.e., 1/5, 1/10, 16


Page 17 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


372

and 1/20). The second best performing method, SEC method based on Virtual-QC

373

algorithm without the calibration of gross error, resulted in 98.75%, 97.81% and

374

96.24% of total peak numbers below the 15% RSD threshold at 1/5, 1/10 and 1/20

375

frequencies, respectively. Moreover, both of the IS and TISC methods had modest

376

calibration abilities, and F didn't affect their calibration results. Nevertheless, the

377

lowest frequencies of QC injection (1/20) had a worst performance compared to 1/5

378

and 1/10 for uncorrected raw data. F had less effects on the correction results of IS,

379

TISC, SEC, and G+S methods. Hence, it is reasonable to keep a F value as 1/20, if

380

G+S correction method is used, which had the 100 % POF with RSD less than 15%

381

and could increase the analysis throughput in a given time.

382 383

Method Application in integrating the metabolic profiling data of nine batches of

384

samples

385

The number of analytical samples was enlarged to 1197 in nine batches analyzed

386

by two different instruments and columns (Figure 1) to further validate the established

387

G+S method, which was compared with the two commonly used sample-based

388

correction methods (i.e. IS and TISC) and two feature-based signal calibration

389

methods (i.e. LoReg and LoMec). It was obvious that G+S method calibrated QCs to

390

a shortest ED and highest r compared to IS1, TISC, LoReg and LoMec (Figure S-4A

391

in the Supporting Information). Furthermore, the RSDs distribution showed the

392

significant increase in peaks number of QC samples with RSD less than 10%, 15%

393

and 20% compared to IS1, TISC, LoReg and LoMec (Figure S-4B in the Supporting

394

Information). Table 1 shows that G+S method also improved the repeatability of five

395

kinds of non-QC replicates (S1~S5) due to the highest r, the shortest ED of 85% PCs

396

in the PCA scores plot and the largest POF with RSDs within 10%, 15%, 20% and 17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

397

Page 18 of 29

30%.

398

After calibrated by G+S approach, the RSDs of external standard ES1

399

(Damascenon) and internal standard IS1 were decreased from 58.89% to 8.65%

400

(Figure 5A and 5B) and from 62.90% to 12.41% (Figure 5C and 5D), respectively,

401

which also indicates that G+S method can effectively reduce these gross and

402

systematic error influences.

403

The signal response drifts of instrument induced by changes of instrument

404

conditions, can be reflected by the variation tendency of the raw signals of external

405

standard (ES) across different batches, because the same amount of ES was added to

406

sample vial before instrument analysis. The greatest signal changes of instrument

407

were observed between batch 7 and batch 8, which was probably induced by changing

408

GC-MS instrument or column (Figure 5A). The instrument parameter drifts (e.g.,

409

detector voltage, RF Gain, RF offset and lens voltage) by using different tuning files

410

could generate a moderate signal changes by observing the responses drifts of ES1

411

between batches 4 vs 5 and batches 6 vs 7 (Figure 5A). The signal of ES1 from batch

412

1 to batch 4 showed a continuous decline pattern, which was closely related with

413

injection order effects, and the changing glass insert and silicone septum had smaller

414

effects on the signal drifts (Figure 5A).

415

In order to visualize the effects of calibration method on the QCs and the non-QC

416

samples, PCA analysis of 1197 samples based on UV scaling was conducted (Figure

417

5E-F). All of the QC samples and five kinds of the non-QC replicates in nine batches

418

before calibration displayed obvious linear changes in PCA score plot (Figure 6E),

419

while good aggregation results of QCs and non-QC replicates can be obviously

420

observed by the calibration of G+S method (Figure 5F). Moreover, this method can

421

also be used to correct nontargeted metabolomics data. Taking our previous 18


Page 19 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


422

nontargeted metabolomics data (n=194, one batch) from GC-MS analysis9 as an

423

example. After calibrated with the suggested method, QCs (n=16) have shorter ED,

424

higher r and lower RSD value (Figure S-5). The RSD of IS (tridecanoic acid) in 194

425

analytical samples was decreased from 8.78% to 7.22%. These data indicate that G+S

426

method can effectively improve the repeatability of nontargeted metabolomics

427

analysis. All the results indicate that the developed G+S method is a universal

428

metabolomics data correction method, which is not only suitable for pesudotargeted

429

but also nontargeted metabolomics data.

430 431

CONCLUSIONS

432

The quality of metabolomics data of large-scale samples set in multiple batches

433

analyzed by different instruments and chromatographic columns could be affected by

434

outliers and systematic errors, which are caused by some unstable environment factors,

435

injection order and batch differences. We developed an novel correction strategy for

436

the integration of large-scale and long-term metabolomics data from pseudotargeted

437

GC−MS by calibrating the outliers and systemic errors and minimizing the frequency

438

of QC injection. Our suggested procedures were as follows,

439

1) correct gross error of large-scale metabolomics data by using Each-5% method

440

to select potential outliers, To-QCn method to further confirm the outlier features, and

441

a linear fitting model of metabolic features to calibrate the outlier variables.

442

2) utilize a feature-based correction algorithm, Virtual-QC, to remove the

443

systematic bias of metabolomics data by using a linear regression model of the

444

intensity values of ω neighborhood QCs of each sample to formulate a correction

445

factor of each feature.

446

3) set the frequency of QC injection in the analytical sequence as one in every 20 19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 29

447

real samples (1/20) if the metabolomics data are corrected by the above mentioned

448

G+S calibration method.

449 450

It should be emphasized that although the pseudotargeted method was used in this study, the strategy suggested is also suitable to nontargeted analysis.

451 452

ASSOCIATED CONTENT

453

Supporting Information

454

Additional information as noted in text includes the Algorithm for gross and

455

systematic errors correction, Tables S-1~S5 and Figures S-1~S-5. This material is

456

available free of charge via the Internet at http://pubs.acs.org.

457 458

AUTHOR INFORMATION

459

Corresponding Authors

460

*E-mail: [email protected] (X.L.); [email protected] (G.W.X.).

461 462

Notes

463

The authors declare no competing financial interest.

464 465

Acknowledgements

466

This study has been supported by the National Grand Project (2014ZX08012-002) of

467

Science and Technology of China, the project (No. 21375127), key project (No.

468

21435006) and the creative research group project (No. 21321064) from the National

469

Natural Science Foundation of China.

470 471 20


Page 21 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


472

REFERENCES

473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512

(1) Fiehn, O. Plant Mol. Biol. 2002, 48, 155-171. (2) Goodacre, R.; Vaidyanathan, S.; Dunn, W. B.; Harrigan, G. G.; Kell, D. B. Trends Biotechnol. 2004, 22, 245-252. (3) Huang, Q.; Tan, Y.; Yin, P.; Ye, G.; Gao, P.; Lu, X.; Wang, H.; Xu, G. Cancer Res. 2013, 73, 4992-5002. (4) Lisec, J.; Schauer, N.; Kopka, J.; Willmitzer, L.; Fernie, A. R. Nat. Protoc. 2006, 1, 387-396. (5) Alonso, A. P.; Goffman, F. D.; Ohlrogge, J. B.; Shachar-Hill, Y. Plant J. 2007, 52, 296-308. (6) Marti, G.; Erb, M.; Boccard, J.; Glauser, G.; Doyen, G. R.; Villard, N.; Robert, C. A.; Turlings, T. C.; Rudaz, S.; Wolfender, J. L. Plant Cell Environ. 2013, 36, 621-639. (7) Zhang, J. T.; Zhang, Y.; Du, Y. Y.; Chen, S. Y.; Tang, H. R. J. Proteome Res. 2011, 10, 1904-1914. (8) Ye, G. Z.; Zhu, B.; Yao, Z. Z.; Yin, P. Y.; Lu, X.; Kong, H. W.; Fan, F.; Jiao, B. H.; Xu, G. W. J. Proteome Res. 2012, 11, 4361-4372. (9) Zhou, J.; Zhang, L.; Chang, Y. W.; Lu, X.; Zhu, Z.; Xu, G. W. J. Proteome Res. 2012, 11, 4351-4360. (10) Zhao, Y. N.; Zhao, C. X.; Lu, X.; Zhou, H. N.; Li, Y. L.; Zhou, J.; Chang, Y. W.; Zhang, J. J.; Jin, L. F.; Lin, F. C.; Xu, G. W. J. Proteome Res. 2013, 12, 5072-5083. (11) Dunn, W. B.; Broadhurst, D.; Begley, P.; Zelena, E.; Francis-McIntyre, S.; Anderson, N.; Brown, M.; Knowles, J. D.; Halsall, A.; Haselden, J. N.; Nicholls, A. W.; Wilson, I. D.; Kell, D. B.; Goodacre, R.; Human Serum Metabolome, C. Nat. Protoc. 2011, 6, 1060-1083. (12) Li, Y.; Ruan, Q.; Li, Y. L.; Ye, G. Z.; Lu, X.; Lin, X. H.; Xu, G. W. J. Chromatogr. A 2012, 1255, 228-236. (13) Kamleh, M. A.; Ebbels, T. M.; Spagou, K.; Masson, P.; Want, E. J. Anal. Chem. 2012, 84, 2670-2677. (14) Li, J.; Hoene, M.; Zhao, X.; Chen, S.; Wei, H.; Haring, H. U.; Lin, X.; Zeng, Z.; Weigert, C.; Lehmann, R.; Xu, G. Anal. Chem. 2013, 85, 4651-4657. (15) Fages, A.; Pontoizeau, C.; Jobard, E.; Levy, P.; Bartosch, B.; Elena-Herrmann, B. Anal. Bioanal. Chem. 2013, 405, 8819-8827. (16) Wang, S. Y.; Kuo, C. H.; Tseng, Y. J. Anal. Chem. 2013, 85, 1037-1046. (17) Guidance for Industry Bioanalytical Method Validation; U.S. Department of Health and Human Services, Food and Drug Administration, Center for Drug Evaluation and Research (CDER), Center for

Veterinary

Medicine

(CVM),

2001,

http://www.fda.gov/downloads/drugs/guidancecomplianceregulatoryinformation/guidances/ucm0701 07.pdf. (18) Yuan, Y.; Khatibisepehr, S.; Huang, B.; Li, Z. AICHE J. 2015, 61, 3232-3248. (19) Bretas, N. G.; Bretas, A. S. International Journal of Electrical Power & Energy Systems 2015, 73, 484-490. (20) Calleri, A.; Bono, A.; Bagnardi, V.; Quarna, J.; Mancuso, P.; Rabascio, C.; Dellapasqua, S.; Campagnoli, E.; Shaked, Y.; Goldhirsch, A.; Colleoni, M.; Bertolini, F. Clin. Cancer Res. 2009, 15, 7652-7657.

513

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 29

Table 1. Percentages of features (POF) with RSDs within 10%, 15%, 20% and 30%, Average r and ED using 85% PCs in the PCA analysis of the five types of non-QC replicates (S1, S2, S3, S4 and S5)

Samples Method

S1

S2

S3

S4

S5

Average distance Pearson correlation in PCA coefficient

POF with

G+S

ED 9.53

r 0.84

RSD

A Novel Strategy for Large-Scale Metabolomics ... - ACS Publications

Recommend Documents