Retrospective Urine Metabolomics of Clinical Toxicology Samples Reveals Features Associated With Cocaine Exposure

1. Introduction

Cocaine is one of the most widely used illicit stimulants globally, associated with diverse toxic effects, including cardiovascular complications, neuropsychiatric syndromes, and long-term addiction [1,2,3,4,5]. Despite its clinical and forensic relevance, the broader metabolic alterations associated with cocaine use remain incompletely understood, especially in real-world clinical settings. Routine urine drug screening relies on immunoassays and subsequent mass spectrometry-based assays to detect known drugs and their primary metabolites [6]. These approaches are practical for confirming exposure but do not capture the full extent of biochemical changes induced by drug use.

Metabolomics, the global profiling of small molecules in biological systems, enables more profound insight into drug-induced metabolic perturbations. Untargeted metabolomics has been increasingly applied in toxicology research to identify novel biomarkers, uncover altered pathways, and better understand drug mechanisms [7,8,9,10,11,12,13]. With advances in liquid chromatography-high-resolution mass spectrometry (LC-HRMS), especially liquid chromatography–quadrupole time-of-flight mass spectrometry (LC-qToF-MS), untargeted data acquisition is now routinely incorporated into clinical toxicology workflows, creating opportunities for retrospective analysis [6].

As part of routine clinical practice, we conduct urine comprehensive drug screening (UCDS) (> 500 cases per month) to identify causative drugs and chemicals for intoxication and to monitor compliance of prescription drugs at the University of Pittsburgh Medical Center Clinical Toxicology Laboratory. In urine comprehensive drug screening, we acquire mass spectral data of analytes in an untargeted manner using LC-qToF-MS with all ion fragmentation scans [6,14]. LC-qToF-MS reveals 5000 - 10000 features per specimen, but our clinical practice focuses only on 50 features of xenobiotics (drugs and their metabolites) per sample at best, leaving the remaining features unannotated and untouched. These unannotated features represent an underutilized resource that may harbor clinically actionable biomarkers and provide novel insights into the metabolic effects of drug exposure.

In this study, we repurposed data from our UCDS, which employs LC-qToF-MS for qualitative toxicological assessment. Although the primary goal of this platform is drug screening, the untargeted data generated by LC-qToF-MS are suitable for downstream metabolomics. Using an archived mass spectrometry dataset, we conducted a retrospective metabolomics study, classifying specimens using qualitative cocaine-enzyme immunoassay (EIA) results as metadata. By comparing the metabolic profiles of EIA-positive and EIA-negative samples, we aimed to identify discriminating features associated with cocaine use.

While cocaine itself is not a novel analyte, our approach demonstrates the value of mining routine clinical toxicology data for metabolomic insights. Specifically, we show how retrospective feature-level analysis can reveal potential biomarkers of cocaine exposure, identify co-eluting compounds and metabolites, and shed light on broader biochemical alterations. We also address analytical challenges, including annotating isomeric features and possible in-source fragmentation, which are common in LC-MS-based metabolomics [15,16,17,18,19].

This study highlights the translational potential of integrating metabolomics into clinical toxicology workflows. By leveraging existing diagnostic data, we provide a model for low-barrier biomarker discovery and a deeper understanding of drug-induced metabolic profiles using real-world clinical samples.

2. Materials and Methods

2.1. Specimens

The UCDS datasets were originally derived from the urine specimens from 363 patients (242 patients with age ≥ 20 years and 121 patients with age ≤ 19 year-old) between Mar and May 2021 and in June 2022 (Table 1). The specimens were sent from emergency departments of UPMC hospitals and clinics in Pittsburgh, PA. These specimens underwent the initial EMIT-II-based qualitative drug screening panel including the one for cocaine (cocaine-EIA) with the cutoff at 300 ng/mL for benzoylecgonine (Siemens Healthineers USA, Malvern, PA) before mass spectrometry-based analysis. One volume of urine specimen was mixed with 4 volumes of distilled water spiked with nitrazepam and tenoxicam as internal standards (40 ng/mL final). Diluted urine specimen was directly injected to LC-qToF-MS (3 μL).

2.2. LC-qToF-MS conditions/settings

Xevo® G2-S QToF mass spectrometer and an ACQUITY UPLC® I-Class system with UPLC HSS C18 column (Waters) was used for UCDS. The acquired mass spectrometry datasets were analyzed using UNIFI® Scientific Information System (Waters) to screen compounds through retention time, monoisotopic mass accuracy, and fragment ions match in the UCDS.

Gradient-elution chromatography in reverse-phase systems was achieved using the mobile phase A (5 mM ammonium formate pH 3.0) and mobile phase B (Acetonitrile containing 0.1% formic acid). The gradient elution was performed as follows in the 15 min run: the column was equilibrated with 87% A for 0.5 minutes, then the concentration of A was decreased to 50% over 9.5 minutes, then decreased to 5% over 0.75 minutes and held for 1.5 minutes, then increased back to 87% over 2 minutes. The column temperature was set at 50.0 ⁰C.

Electron spray ionization (ESI) in the positive ionization mode was used. The ESI settings are as follows: capillary 0.80 kV, sample cone voltage 25 V, source temperature 150 ⁰C, desolvent gas 800 L/h at 400 ⁰C, and cone gas 20 L/h. MS data of both precursor and product ions of 50 – 1000 Da mass range were collected in a non-targeted manner using MSE technology with the following settings (Collision cell setting: low energy 6 eV, high energy ramp 10 to 40 eV; Scan time: 0.1 sec; Data format: Continuum).

2.3. Data export and preprocessing

An overview of the overall data processing and analysis workflow, from raw mass spectrometry data export through statistical analysis, is presented in Figure 1. MS datasets in (.uep) files were first exported from UNIFI as (.raw) files. After removing any personal identifiers and lock mass data, these files were converted to (mzML) files in the centroid mode using ProteoWizard msConvert [20].

From the MS datasets (.mzML), peak detection, deconvolution, alignment, and putative annotation of the aligned features were performed using MS-DIAL version 4.90 [21], resulting in a peak aligned table comprising 14,883 features across 363 urine specimens. As for MS libraries, 48189 MS records (both experimental and predicted) of human urine metabolites were downloaded from Human Metabolome Database Ver 5 [22], converted into the NIST msp format library, and added to the existing 56694 MS records in MS libraries (MassBank, MassBank EU, ReSpect, GNPS, Fiehn HILIC, CASMI2016 MetaboBASE, RIKEN PlaSMA authentic standards, RIKEN PlaSMA bio-MS/MS (MSI level 1,2,3, or 4) from plant tissues, Karolinska institute and Gunma (GIAR) zic-HILIC deconvoluted MS2 spectra in data independent acquisition, Fiehn Pathogen Box, and Fiehn/Vaniya natural product library) available for MS-DIAL. Mass accuracy parameters for feature detection and spectral deconvolution were set to 0.01 Da for MS1 and 0.05 Da for MS2. For MS/MS spectral annotation, the accurate mass tolerances were set to 0.025 Da for MS1 and 0.25 Da for MS2, reflecting the modest mass accuracy due to the absence of lock mass correction. Peak detection was performed with a minimum peak height of 5,000 amplitude and a mass slice width of 0.1 Da. For peak alignment, the retention time tolerance was set to 0.25 minutes, and the MS1 tolerance was 0.025 Da. From the peak-aligned table generated by MS-DIAL, the peak intensity data matrix was constructed by applying row-wise normalization using probabilistic quotient normalization (PQN). The median of the entire dataset was used as the reference to minimize the effects of varying urine concentrations and improve comparability across specimens [23].

2.4. Data analysis

The peak intensity data matrix was further processed for data exploration and analysis using MetaboAnalyst (both the web version and R-version or MetaboAnalystR), a unified platform developed for comprehensive metabolomics data analysis[24,25,26,27,28]. Statistical analyses were selected to balance sensitivity for feature discovery (e.g. volcano plot) with robustness to complex, high-dimensional data and class imbal-ance (e.g., Random Forest), ensuring comprehensive identification and valida-tion of cocaine-associated metabolic features. Both column-wise normalization and data filtering were first applied to the peak intensity data matrix to generate the normalized data matrix, before statistical analyses were performed using MetaboAnalyst. Column-wise normalization with log transformation and autoscaling was performed to make each feature comparable. Data filtering was performed to remove non-informative features with very small or near-constant values to reduce the number of features to 5000, the maximum number of the features that MetaboAnalyst can analyze. The cocaine-EIA results were used as metadata to categorize the specimens for further statistical analysis.

Point biserial correlation analysis was performed with MetaboAnalystR to explore the correlations between metabolomics features and EIA-results (binary metadata). The significant features with FDR < 0.05 were selected.

Volcano plot was performed with MetaboAnalystR to combine the results from fold change analysis and t-tests. Fold changes are calculated as the ratios between two group means > 2 using the original data before normalization. T-tests were also conducted for each feature with p-value threshold < 0.05.

Significance analysis of microarrays and metabolites (SAM) [29] and Empirical Bayesian analysis of microarrays and metabolites (EBAM) [30] were performed with MetaboAnalystR to select the significant features based on false discovery rate (FDR). The SAM plot is a scatter plot of the observed relative difference versus the expected relative difference estimated by data permutation. The EBAM was performed to select the significant features at posterior delta of 0.9 and FDR < 0.05.

Partial least squares discriminant analysis (PLS-DA) was performed with MetaboAnalystR to select significant features to predict EIA results (metadata-based classification) by calculating the weighted sum of absolute regression coefficient for the overall components. The selected features were evaluated further only when statistically significant results were obtained by permutation testing 2000 times.

Random Forest was performed using R packages (caret and randomForest) after class imbalance of samples was corrected by combination of random oversampling of the non-dominant class samples and undersampling of the dominant class samples using R package ROSE.

Cluster analysis was performed by calculating the feature-to-feature Pearson correlation coefficients from the transposed peak intensity matrix of cocaine-associated features for the cocaine-EIA positive specimens. The resulting correlation matrix was visualized as a heatmap using the pheatmap package in R and hierarchical clustering was performed using Ward’s method (ward.D2). Features were grouped into seven clusters (Clusters A–G; Table 2A–G) by high-level dendrogram slicing (k = 7). Middle-level (k = 14) and low-level (k = 50) dendrogram slicing were additionally applied to identify subclusters of co-varying features that may reflect structurally related compounds or analytical associations.

Cluster C (Table 2C)			Cluster D (Table 2D)
91.05293_0.83	135.05078_2.818		91.05488_2.605	94.03997_0.693	116.22058_0.802
168.11021_0.847	230.1496_1.057		149.06212_3.642	105.03697_0.875	200.09056_1.02
200.12935_0.795	230.29869_1.051		186.01846_4.664	290.10568_6.327	200.09056_1.294
441.17444_0.915	168.10132_1.248		121.04194_0.825	522.40491_12.074	168.06374_1.28
117.05989_0.864	168.10931_1.568		122.07076_0.782	107.04818_1.115	186.08173_1.203
139.12474_1.347	184.10083_1.021		137.19272_0.733	546.40143_11.75	200.09158_0.796
163.12866_0.849	186.11212_1.025		150.05968_0.839	131.11659_0.693
177.34515_0.976	186.11491_0.834		152.07004_2.343	402.37381_11.955
198.11752_1.083	192.10089_1.04		191.02417_0.942	534.41858_11.965
214.10306_1.408	208.09537_1.018		140.07437_0.729	172.10063_1.902
216.12822_1.072	177.10405_0.988		185.07774_1.451	173.08467_2.303
150.09793_0.812	177.14621_0.68		201.08185_1.004	187.06279_0.952
162.108_2.008	177.23868_0.976		265.11493_0.706	238.84851_0.721
186.11717_1.818	179.12364_0.927		155.07788_0.851	256.82074_0.732
258.08966_1.448	193.10092_0.876		182.08345_0.819	384.85068_0.716
339.16333_0.701	193.13478_0.697		268.09927_0.854	400.82483_0.724
	193.24315_0.872
	357.22891_0.88
	214.15097_0.935
	240.10555_1.076
	264.12234_0.734
	475.2739_9.493
	504.31165_11.765
	218.08513_1.244
	221.08717_1.055
	282.15509_4.618
	309.13354_1.511
	331.22858_4.093
	388.22858_3.758

Cluster E (Table 2E)			Cluster F (Table 2F)			Cluster G (Table 2G)
104.10567_0.733	104.10567_11.995		106.06516_1.994	153.13197_4.525		116.05332_11.649
104.2101_0.734	104.10775_11.868		256.18127_1.385	311.15509_4.525		135.04501_11.694
152.065_2.776	104.10775_11.447		258.08966_1.149	378.18292_2.252		311.34421_11.694
207.15585_2.426	149.05659_1.326		300.20895_1.719	165.07965_0.797		312.16464_12.233
233.17049_3.226	213.12489_0.836		330.18475_1.651	224.12656_0.967		312.16464_11.693
125.06116_2.653	395.29956_11.649		338.26901_11.447	277.12286_2.952		527.40106_12.155
223.11052_1.41	104.10567_0.978		498.27695_6.524	292.15628_2.97		593.33722_10.015
323.14807_2.849	142.08922_1.821		554.4505_12.18	519.40143_12.072		659.48175_11.901
382.83783_0.703	149.06207_4.254		680.44757_11.809	544.40845_11.628		142.12054_0.995
449.16858_4.928	169.10295_0.889		124.07899_1.053	639.43219_12.067		203.06929_2.133
470.35092_11.745	165.07159_4.713		253.05179_3.2	180.08086_2.536		237.07332_1.753
462.8587_0.722	165.07463_5.316		278.14926_2.887	235.13065_0.829		303.16507_3.432
467.10083_11.731	218.0378_1.446		301.11752_4.474	195.12511_1.716		146.05814_1.87
	250.86957_0.719		329.21564_3.3	246.10471_0.754		188.06906_1.861
	256.17471_1.742		332.25024_4.576	256.18103_6.388		247.14194_1.873
	350.02979_0.703		183.09117_0.77	239.15465_1.732		247.29958_1.876
	468.81296_0.725		201.16733_0.774	283.18033_2.464		182.123_0.87
	498.89853_0.711		237.131_1.815	300.21008_1.962		200.13109_4.253
			290.20178_3.419	344.23004_1.752		200.13182_1.007
			362.12997_4.903			200.13376_1.911
			220.12494_1.151			292.1264_4.051
			220.27379_1.155			304.15652_4.399
			242.0983_1.15
			247.09369_1.997
			254.17285_0.906
			261.15073_2.827
			330.17572_5.651
			354.23468_5.676
			493.37396_11.814
			544.42969_11.838
			553.40015_11.93

2.5. Feature annotations

MS-FINDER version 3.61[31,32] was utilized to further elucidate the putative chemical structures of selected features based on MS/MS fragmentation patterns and elemental composition. Minor metabolites and impurities of recreational drugs, including cocaine, tobacco, opioids, and benzodiazepines, and major novel psychoactive substance (NPS), were also fed into MS-FINDER through the user-defined database prepared manually. SMILES and InChIKey of these metabolites were prepared manually using ChemDraw (Ver 23.1.2) (PerkinElmer, Waltham, MA) if SMILES and InChIKey are not available in the public database (Human Metabolome Database and Pubchem).

MS1 mass spectra and extracted ion chromatograms (EICs) were visualized using MZmine 4.78 [33,34]. The retention times were also predicted using R package Retip [35], which was trained using our retention time data of the compounds.

3. Results

Among the 14,883 features in the dataset, 262 were initially identified as significantly associated with cocaine EIA results by at least one statistical analysis. However, initial putative annotation of these features using MS-DIAL and MS-FINDER met with only limited success, presumably because of the limited inclusion of minor cocaine metabolites, cocaine-related impurities and their metabolites including phase II metabolites in the public databases [36]. Thus, we created a user-defined database covering these analytes for MS-FINDER. In-source fragments might also be mistaken as precursor ions [15,37]. To better interpret these features, we performed cluster analyses on the 262 features in the normalized data matrix, based on the assumption that structurally related compounds—such as parent drugs, their immediate metabolites, and associated ion signatures—would cluster together. The features were grouped into seven high-level clusters (Clusters A–G), corresponding to Table 2A through 2G. This hierarchical clustering approach was designed to reflect co-variation patterns and facilitate the interpretation of unannotated or partially annotated signals (Table 2 and Figure 2).

Based on these analyses, we attempted to annotate 37 more significant features that were selected by more than one analysis and, thus, more strongly associated with the cocaine-EIA results. These 37 significant features and their annotations are summarized in Table 3. These 37 features are classified as cocaine and its metabolites, cocaine impurities, nicotine and its metabolites, opioid (norfentanyl) and tryptophan/serotonine metabolite (3-hydroxy-tryptophan) (Table 3).

3.1. Features 200.12935_0.795, 200.13182_1.007 and 200.13109_4.253

These features likely correspond to the molecular formula of C10H17NO3 ([M+H]+ 200.12812). Possible chemical structures for C10H17NO3 elucidated by MS-FINDER include ethyl norecgonine or methyl ecgonine among the cocaine metabolites, and methyl ecgonine is the top candidate for these three features. The RT of ethyl norecgonine is predicted to be slightly longer than that of methyl ecgonine (RT: 0.80 min, determined experimentally) (Table 4) and thus, the feature 200.13182_1.007 is putatively annotated to be ethyl norecgonine. Even though the annotation of the feature 200.13109_4.253 is unclear, this feature is clustered in the low slicing together with cocaine 304.15652_4.399 and its metabolites including the feature 200.13182_1.007 (putative ethyl norecgonine) in Cluster G (Table 2G), and this feature 200.13109_4.253 is likely another cocaine metabolite or possibly an in-source fragment of another cocaine metabolite.

3.2. Features 304.15652_4.399, 304.16803_4.254

These features likely correspond to the molecular formula of C17H21NO4 ([M+H]+ 304.15434). Among the structures proposed by MS-FINDER, one of the top candidates includes cocaine, which has the RT of 4.46 min, consistent with these features. The feature 304.15652_4.399 is clustered in the low slicing together with other cocaine metabolites such as methyl ecgonine 200.12935_0.795 in Cluster G (Table 2G), as discussed previously, whereas the feature 304.16803_4.254 is clustered together in the low slicing with an illicit cocaine impurity (Cinnamoylcocaine, 330.17252_5.513) and in the middle slicing with a cocaine metabolite (Benzoylecgonine, 290.14951_3.314) in Cluster A (Table 2A), as discussed below.

These two closely eluting features with almost the same m/z values are likely representation of cocaine. Manual inspection of the extracted ion chromatograms (EICs) confirmed that these features share a single chromatographic peak shape with slightly offset apexes. Given the absence of lock mass correction, minute m/z drift and centroid rounding likely caused the algorithm to assign separate features to the same molecular ion.

3.3. Feature 182.123_0.87

This feature likely corresponds to the molecular formula of C10H15NO2 ([M+H]+ 182.11756). Possible chemical structures for C10H15NO2 elucidated by MS-FINDER include ethyl norecgonidine or methyl ecgonidine (anhydroecgonine methylester) among the cocaine metabolites. While methyl ecgonidine is the top match based on structure scoring, its experimentally determined RT is 1.10 min, which does not match the observed RT of 0.87 min. Since ethyl norecgonidine is predicted to have an even longer RT than methyl ecgonidine (Table 4), it is also an unlikely match for this feature.

An alternative and more plausible explanation is that this feature represents an in-source loss of water ([M+H–H2O]+) from a parent ion at m/z 200.1291, corresponding to methyl ecgonine 200.12935_0.795. Although methyl ecgonine and its in-source fragment are expected to co-elute, they can be assigned slightly different retention times (0.80 min vs. 0.87 min) by MS-DIAL. This difference likely reflects the way MS-DIAL independently detects and deconvolutes each ion feature based on apex intensity within its respective m/z window, rather than a true chromatographic separation, in the peak alignment process [21]. Furthermore, this feature was observed in the same peak-slicing cluster as cocaine 304.15652_4.399 and other known cocaine-related metabolites in Cluster G (Table 2G), supporting its biological relevance to cocaine metabolism.

3.4. Features 186.11491_0.834, 186.11212_1.025 and 186.11717_1.818

These features likely correspond to the molecular formula of C9H15NO3 ([M+H]+ 186.11247). Possible chemical structures for C9H15NO3 elucidated by MS-FINDER include ecgonine or methyl norecgonine among the cocaine metabolites. Ecgonine is the top candidate for these three features, whereas methyl norecgonine is one of the top hits for the feature (186.11212_1.025). The RT of methyl norecgonine is predicted to be longer than that of ecgonine (RT: 0.80 min determined experimentally) (Table 4). Thus, ecgonine is likely the molecule of the feature 186.11491_0.834, and methyl norecgonine would be the molecule of the feature 186.11212_1.025.

To seek alternative explanation for the feature 186.11717_1.818, we manually reviewed the aligned peak table generated by MS-DIAL for other cocaine metabolites [8,38] eluting around 1.8 min and identified another feature 306.14551_1.781 for which hydroxy-benzoylecgonine (C16H19NO5, [M+H]+ 306.13360) is elucidated to be the top candidate by MS-FINDER. This feature was filtered out by MetaboAnalyst; therefore, it is not included in the normalized data matrix or cluster analysis. The MS2 spectrum of this feature includes an ion 186.1163, which can be generated after the neutral loss of C7H4O2, which corresponds to the hydroxybenzaldehyde moiety, leaving the ecgonine moiety (Figure 3A). Peak patterns of these features look similar to each other in the EIC of MS1 around 1.8 min (Figure 3B,C).

In-source fragmentation is a rather common phenomenon in the ESI-based analysis of natural molecules [39]. Ion fragments generated through in-source fragmentation and a collision cell are comparable [40,41]. Thus, we speculate that the feature 186.11717_1.818 might be an in-source fragment of hydroxybenzoylecgonine.

3.5. Features 330.17252_5.513 and 330.17572_5.651

These features likely correspond to the molecular formula of C19H23NO4 ([M+H]+ 330.16998). Possible chemical structures for C19H23NO4 elucidated by MS-FINDER include cis- and trans-cinnamoylcocaine, tropane alkaloids in coca leaves and impurities of illicit cocaine products [42,43].

3.6. Features 168.11021_0.847, 168.10132_1.248, and 168.10931_1.568

These features likely correspond to the molecular formula of C9H13NO2 ([M+H]+ 168.10191), which is either ecgonidine or methyl norecgonidine among the cocaine metabolites. The RT of methyl norecgonine is predicted to be slightly longer than that of ecgonidine (RT: 0.84 min, determined experimentally)(Table 4); thus, the feature 168.11021_0.847 would be annotated to be ecgonidine whereas the feature 168.10132_1.248 or 168.10931_1.568 would be methyl norecgonine.

Both features (168.10132_1.248 and 168.10931_1.568) are grouped together in the low slicing with a feature 186.11212_1.025, which is suspected to be methyl norecgonine in Cluster C (Table 2C) (see above). Notably, the mass difference between this feature and the 168 m/z features is approximately 18 Da, suggesting a possible in-source loss of water ([M+H–H2O]+) for methyl norecgonine. This raises the possibility that 186.11212_1.025 may be the intact parent ion, while 168.11021_0.847 and 168.10132_1.248 could include contributions from its in-source fragment of 186.11212_1.025.

Furthermore, this m/z 168 signal also matches several isomeric compounds unrelated to cocaine metabolism, including p-synephrine, phenylephrine (m-synephrine), and 3-methoxytyramine. Synephrine is occasionally reported as an adulterant in street cocaine and is also found in weight-loss supplements [44], while phenylephrine is a common over-the-counter decongestant. 3-Methoxytyramine, an immediate dopamine metabolite, is biologically relevant in the context of cocaine use, which alters dopaminergic signaling. Increased urinary levels of 3-methoxytyramine have been reported during early abstinence [45].

Given the physicochemical similarity with the shared formula, overlapping MS/MS spectra, and similar RTs (p- and m-synephrine for 0.83 min and 3-methoxytyramine for 0.98 min experimentally) and possible occurrence of in-source fragmentation, these features, especially 168.11021_0.847 and 168.10132_1.248, might possibly represent composite signals from multiple co-eluting compounds and an in-source fragment of another ion.

3.7. Feature 214.15097_0.935

These features likely correspond to the molecular formula of C11H19NO3 ([M+H]+ 214.14377), for which MS-FINDER suggests ethyl ecgonine as a top candidate molecule among the cocaine metabolites. The RT of ethyl ecgonine is 0.93 min (experimentally determined), matching that of this feature. Thus, the feature 214.15097_0.935 is likely the representation of ethyl ecgonine.

3.8. Feature 233.17049_3.226

These features likely correspond to the molecular formula of C14H20N2O ([M+H]+ 233.16484). Possible chemical structures for C14H20N2O elucidated by MS-DIAL and MS-FINDER include norfentanyl, a dealkylated metabolite of fentanyl, but not any cocaine metabolites. The RT of norfentanyl is 3.25 min (experimentally determined), matching that of this feature. Thus, the feature 233.17049_3.226 is likely the representation of norfentanyl. Opioids such as fentanyl are often co-ingested with cocaine in the fourth wave of the opioid crisis [46,47], supporting this finding.

3.9. Feature 292.1264_4.051

This feature likely corresponds to the molecular formula of C15H17NO5 ([M+H]+ 292.11795), for which hydroxy-benzoylnorecgonine is suggested among the cocaine metabolites by MS-FINDER. There are four isomers possible, based on the hydroxylation site (o-, m-, p-, or N-), which was not discernible by mass spectra in this case. There are four isomers possible, based on the hydroxylation site (o-, m-, p-, or N-); MS-FINDER lists N-hydroxy-benzoylnorecgonine as the top molecule among cocaine metabolites, but it also lists other isomers among the top metabolites as well. N-hydroxy-benzoylnorecgonine elutes at 3.74 min and (o-, m-, or p-)hydroxy-benzoylnorecgonine elutes at 2.48 min in a similar but slightly shorter gradient RPLC program (10 min instead of 15 min in our program) [38]. Thus, N-hydroxy-benzoylnorecgonine would be the putative annotation for the feature 292.1264_4.051, but (o-, m-, or p-)hydroxy-benzoylnorecgonine cannot be ruled out as well.

3.10. Feature 290.14951_3.314

This feature likely corresponds to the molecular formula C16H19NO4 ([M+H]+ 290.13869), which matches both benzoylecgonine (experimentally observed RT: 2.90 min) and norcocaine (RT: 4.54 min) among cocaine metabolites. Based on the MS2 pattern and RT, this feature aligns more closely with benzoylecgonine. However, the RT is approximately 0.4 min later than expected for benzoylecgonine.

The closer examination suggests that this feature may represent the tailing portion of a large benzoylecgonine peak (the feature 290.14487_2.854). Notably, this main peak was removed during data filtering by MetaboAnalyst, likely due to its detection across both cocaine-EIA-positive and -negative specimens, possibly due to its rather ubiquitous presence of this feature across both cocaine-EIA-positive and -negative specimens because of its relatively high cut-off level (300 ng/mL). In the meantime, MS-DIAL seems to detect the tail of the large benzoylecgonine peak as a distinct feature, a known limitation of peak detection algorithms in the presence of broad or asymmetric chromatographic peaks [48,49].

3.11. Feature 221.08717_1.055

This feature likely corresponds to the molecular formula of C11H12N2O3 ([M+H]+ 221.092069). Possible chemical structures for C11H12N2O3 elucidated by MS-FINDER include 5-hydroxy-tryptophan (RT: 0.90 experimentally), an immediate precursor of the neurotransmitter serotonin generated from tryptophan [50], but no cocaine metabolites are suggested by MS-FINDER. Cocaine use is associated with altered tryptophan and serotonin metabolism [51], with elevated plasma serotonin levels observed in abstinent cocaine users [52]. Thus, the feature 221.08717_1.055 is annotated as 5-hydroxy-tryptophan.

3.12. Features 193.10092_0.876, and 193.24315_0.872

This feature likely corresponds to the molecular formula of C10H12N2O2 ([M+H]+ 193.09715). Possible chemical structures for C10H15NO2 elucidated by MS-FINDER include hydroxycotinines and cotinine N-oxide. Among these metabolites, 3-hydroxycotinine is the dominant nicotine metabolite in the urine [53]. The RT is experimentally determined to be 0.90 min, whereas that of cotinine N-oxide is 0.96 min. Thus, the feature 193.10092_0.876 is likely 3-hydroxycotinine, but minor contributions from other isomeric nicotine metabolites (e.g., 5-hydroxycotinine and cotinine-N-oxide) cannot be ruled out.

Regarding the feature 193.24315_0.872, both MS-DIAL and MS-FINDER do not suggest any candidate molecules. This feature has an almost identical RT to that of the feature 193.24315_0.872 (3-hydroxycotinine) but exhibits a slightly higher m/z (Figure 4A,C). It is also grouped together with the features of nicotine metabolites, including 3-hydroxycotinine 193.10092_0.876, in the low slicing in Cluster C (Table 2C). Furthermore, the injection of authentic 3-hydroxycotinine yields both m/z 193.0861 and m/z 193.2396 at RT 0.90 min (Figure 5A). Based on these facts, we speculate that this higher-mass feature 193.24315_0.872 would be derived from 3-hydroxycotinine, rather than being an unrelated artefact, even though there is no plausible explanation for this possible mass shift of ~0.14 Da.

3.13. Feature 179.12364_0.927

This feature likely corresponds to the molecular formula of C10H14N2O ([M+H]+ 179.117889). Possible chemical structures for C10H14N2O elucidated by MS-FINDER include 2-hydroxynicotine and nicotine-N-oxide. Among these metabolites, nicotine-N-oxide has been known as a minor primary metabolite of nicotine excreted in urine, whereas 2-hydroxynicotine has not been identified in human urine. The RT of nicotine-N-oxide is experimentally determined to be 0.96 min. Thus, this feature is identified as nicotine-N-oxide.

3.14. Features 177.10405_0.988, and 177.23868_0.976

This feature likely corresponds to the molecular formula of C10H12N2O ([M+H]+ 177.102239). Possible chemical structures for C10H12N2O elucidated by MS-FINDER include cotinine, the dominant nicotine metabolite in the urine [53]. The RT is experimentally determined to be 1.0 min; thus, the feature 177.10405_0.988 is cotinine.

Regarding the feature 177.23868_0.976, both MS-DIAL and MS-FINDER do not suggest any candidate molecules. This feature has an almost identical RT to that of the feature 177.10405_0.988 (cotinine) but exhibits a slightly higher m/z (Figure 4B,D). It is also grouped together with the features of nicotine metabolites, including 3-hydroxycotinine 193.10092_0.876 and cotinine 177.10405_0.988, in the low slicing in Cluster C (Table 2C). Furthermore, the injection of authentic cotinine yields both m/z 177.1052 and m/z 177.2387 at RT 1.05 min (Figure 5). Based on these facts, we speculate that this higher-mass feature 177.23868_0.976 would be derived from cotinine, similar to the feature 193.24315_0.872 for 3-hydroxycotinine (see above). This phenomenon appears to be characteristic of cotinine and 3-hydroxycotinine, even though there is no plausible explanation for this possible mass shift of ~0.14 Da.

4. Discussion

This retrospective untargeted metabolomics analysis of routine UCDS data revealed the features significantly associated with cocaine exposure. Table 3 summarizes annotated features, including known cocaine metabolites and related impurities, nicotine metabolites, opioid metabolite (norfentnayl), and tryptophan/serotonin metabolite (3-hydroxy-tryptophan).

Various features were annotated as known or putative cocaine metabolites, consistent with expected urinary excretion patterns in cocaine users. Notably, one significant feature is likely ethyl ecgonine. Cocaine and alcohol are frequently used together by cocaine users [54,55,56]. While ethyl ecgonine has been previously reported in the context of ethanol and cocaine co-exposure [8], it has not been widely utilized as a biomarker of co-ingestion compared to cocaethylene, a toxic cocaine metabolite formed in the presence of ethanol [57]. However, the corresponding feature 318.17764_5.403 for cocaethylene (RT 5.55 experimentally) was filtered out by MetaboAnalyst during preprocessing as a non-informative feature due to low prevalence and signal intensity across most specimens. Our data indicates ethyl ecgonine as a useful biomarker of ethanol and cocaine co-exposure.

In addition to cocaine-related compounds, other classes of exogenous metabolites were also annotated among the features significantly associated with cocaine-EIA. For example, co-use of tobacco products was evident in the dataset, with features corresponding to hydroxycotinine, cotinine, nicotine-N-oxide, and related metabolites. This aligns with known behavioral co-use patterns [58]. Additionally, norfentanyl, a major fentanyl metabolite, was also included among the significant features. Recent studies indicate both cross-contamination of fentanyl within illicit cocaine products and intentional co-use of cocaine and fentanyl, consistent with this finding [46,47,59,60,61,62,63,64].

Besides exogeneous metabolites, the detection of 5-hydroxy-L-tryptophan as a discriminating feature between cocaine-EIA-positive and cocaine-EIA-negative groups is noteworthy. Cocaine alters tryptophan and serotonin metabolism [51]. Plasma serotonin levels are elevated among abstinent cocaine users, and psychiatric comorbidity among them is related to higher plasma serotonin levels than the abstinent cocaine users without psychiatric comorbidity[52]. As 5-hydroxy-L-tryptophan is an immediate precursor of serotonin, it would be worth investigating if urinary 5-hydroxy-L-tryptophan levels could be used for stratification of cocaine-addicted patients for risk of psychiatric symptoms during abstinence from cocaine.

This study has several limitations inherent to untargeted metabolomics workflows by LC-HRMS. The first issue is annotation uncertainty, an intrinsic limitation to metabolomics research [65,66], caused by the limited availability of authentic standards for most metabolites and the presence of isomeric molecules with shared molecular formulas. These isomeric molecules often yield similar MS/MS fragmentation due to conserved substructures, complicating confident differentiation. For example, the molecule formula of C9H13NO2 ([M+H]+ 168.1019) can be exogeneous metabolites (e.g., cocaine metabolites (ecgonidine or methyl norecgonidine) and p/m-synephrine) or endogenous metabolites (e.g., 3-methoxytyramine), and these molecules have similar MS/MS spectra and RTs. Indeed, synephrine was identified along with coca-alkaloids in the Mariani wine, a popular 19th century tonic wine with coca leaf extracts, by LC-high resolution-MS [67], but an alternative annotation of that feature might be either ecgonidine or ISF of ecgonine. Even targeted data acquisition with the multiple reaction monitoring mode might lead to misidentification of isomeric molecules in LC-MS/MS analysis[68,69]. Authentic chemicals, orthogonal separation techniques (e.g., ion mobility), or confirmatory NMR analysis are required to overcome this challenge in LC-HRMS-based metabolomics.

Potential analytical artifacts are other issues. In-source fragmentation can produce artifactual features whereas peak tailing can slightly shift RT values. Both phenomena can lead to false annotation of the features [15,70]. Furthermore, our MS dataset were centroided without lock mass correction, contributing to reduced mass accuracy and potential feature duplication. This was why we adapted rather generous mass tolerance parameters in MS-DIAL and MS-FINDER.

Other limitations include potential bias in the collected data and data preprocessing. Our MS dataset was generated using positive ESI mode, which is suitable for providing comprehensive coverage of xenobiotics, most of which are basic and neutral drugs. However, positive ESI mode might not be ideal for analysis of acidic metabolites, as they typically ionize better in negative ESI mode [71,72]. Thus, our dataset might not provide comprehensive coverage of acidic metabolites.

Other limitations of our study relate to potential bias secondary to sample preparation and normalization strategy. The mass spectrometry dataset was acquired using a dilute-and-shoot approach (1:4 dilution) of urine, which avoids analyte loss and facilitates broad metabolome coverage. However, this approach is also susceptible to matrix effects and ion suppression, especially in complex biological matrices like urine [14]. Although positive ESI mode is generally less prone to ion suppression than negative ESI mode when analyzing urine specimens for drug testing [73], matrix-related variability may still affect features. Additionally, we applied PQN to mitigate variability in urine concentration across specimens [23], but no urinary creatinine-based normalization was performed. As a result, residual bias due to differences in hydration status may persist despite normalization efforts, and thus, the peak heights might not accurately reflect the amounts of the analytes in some specimens.

Despite these limitations, our study underscores the translational potential of repurposing clinical toxicology data generated by LC-qTof-MS for metabolomics research, offering a feasible approach to identifying clinically relevant but often overlooked biomarkers of drug exposure and drug-related metabolic disturbances in real-world settings.

5. Conclusions

This study demonstrates the utility of repurposing routine LC-qToF-MS data from clinical toxicology workflows for untargeted metabolomics analysis. By retrospectively analyzing UCDS data, we identified features associated with cocaine exposure, including known metabolites, potential indicators of co-exposures (e.g., ethyl ecgonine), and alterations in endogenous metabolic pathways such as tryptophan metabolism. While many of these markers are not entirely novel, they remain underutilized in clinical and forensic toxicology. Our findings emphasize the potential of routine mass spectrometry data as a resource for exposure biomarker discovery and for expanding our understanding of drug-related metabolic perturbations in real-world populations. Further validation studies using authentic standards and complementary analytical platforms are warranted to refine the biological interpretation of these features and assess their clinical relevance.