<jupyter_start><jupyter_text>Roest Lab Jupyter Notebook for Exploratory Data AnalysisThis notebook is configured to the `roest-lab` virtual environment. Please ensure you have the `roest-lab` kernel activated when using this notebook. Required PyPI packages- `numpy`- `pyopenms` Auxiliary PyPI packages for data analysis- `matplotlib`<jupyter_code>print("Welcome to the Roest Lab Exploratory MS Jupyter Notebook!")

# Set common directories
import os

BASE_DIR = os.path.abspath(os.getcwd())
DATA_DIR = os.path.join(BASE_DIR, "data", "mzML")

# From pyOpenMS tutorial page - https://pyopenms.readthedocs.io/en/latest/first_steps.html
from pyopenms import MSExperiment, MzMLFile

exp = MSExperiment()
MzMLFile().store(os.path.join(DATA_DIR, "testfile.mzML"), exp)

# Get help with Python objects in pyOpenMS
# help(MzMLFile)
# help(MSExperiment)<jupyter_output><empty_output><jupyter_text>File reading<jupyter_code>import os
from urllib.request import urlretrieve

# Download small mzML file
github_content = "https://raw.githubusercontent.com/OpenMS/pyopenms-extra/master"
tiny_mzml_filepath = os.path.join(DATA_DIR, "tiny.mzML")
urlretrieve(f"{github_content}/src/data/tiny.mzML", tiny_mzml_filepath)

from pyopenms import MSExperiment, MzMLFile

exp = MSExperiment()
# Load file
MzMLFile().load(tiny_mzml_filepath, exp)

# tiny_mzml_filepath is loaded onto exp by reference

# View help of exp object
# help(exp)

# We can get information about spectra and chromatograms using the following methods
print("exp.getNrSpectra() -->", exp.getNrSpectra())
print("exp.getNrChromatograms() -->", exp.getNrChromatograms())

# Iterate through spectra
for spectra in exp:
    print("MS Level:", spectra.getMSLevel())

# Use index to get spectrum
print("MS Level:", exp[1].getMSLevel())

# Obtain raw peaks
spec = exp[1]
mz, intensity = spec.get_peaks()
print(list(zip(mz, intensity)))
print(sum(intensity))<jupyter_output>[(100.08453527754324, 8568.404), (101.07125201017237, 19943.791), (101.08801011452718, 16450.7), (102.05533609439806, 20074.756), (110.07160430250622, 22257.0), (111.1936654237387, 4178.48), (112.09515353098757, 5585.1104), (114.11091451437406, 643005.56), (115.10800882752966, 458708.97), (116.11143034474966, 182238.38), (117.11465545916923, 206543.3), (141.06555481567577, 9189.294), (143.09742966176154, 4464.8716), (145.10718248218157, 25135.135), (145.113110026963, 28470.342), (159.07599691296252, 35307.57), (160.49979928401467, 4270.829), (163.1178193193383, 4036.727), (165.06681811261896, 4227.5884), (167.08098942386886, 4911.2964), (171.77136438028302, 4491.086), (174.13371022164173, 43575.543), (174.13955413286106, 54625.406), (181.06051384686805, 7942.771), (183.07650026110298, 7059.7827), (183.11235221579705, 20127.709), (187.07118767410628, 23547.564), (197.12820284477095, 7912.1187), (199.0711043070118, 25906.55), (200.1028951149372, 21087.969), (200.1391767484103, 7945.1797)[...]<jupyter_text>Total ion current calculationHere, we will apply what we have learned to calculate the total ion current (TIC). The TIC represents the summed intensity across the entire range of masses being detected at every point in the analysis. Basically, we calculate the total ion current of the whole experiment.<jupyter_code># Define a function that calculates TIC for a given ms level


def calcTIC(exp, mslevel):
    tic = 0
    # Iter through all spectra of the experiment
    for spec in exp:
        # Only calc TIC for matching MS1 spectra
        if spec.getMSLevel() == mslevel:
            mz, i = spec.get_peaks()
            tic += sum(i)
    return tic

# Now let's use our function to calculate the TIC
print(calcTIC(exp, 1))
# A list comprehension version of the calcTIC function above
print(sum([sum(s.get_peaks()[1]) for s in exp if s.getMSLevel() == 1]))
print(calcTIC(exp, 2))<jupyter_output>240.0
240.0
110.0<jupyter_text>TIC chromatogramThe total ion current is visualized over the retention time, to allow for the inspection of areas with general high intensity (usually multiple analytes were measured there). This can help the experimentalist to optimize the chromatography for a better seperation in a specific area.<jupyter_code># Let's use matplotlib to visualize our TIC chromatogram
import matplotlib.pyplot as plt
from urllib.request import urlretrieve

# retrieve MS data
gh = "https://raw.githubusercontent.com/OpenMS/pyopenms-extra/master"
output_ms_data = os.path.join(DATA_DIR, "ms_data.mzML")
urlretrieve(gh + "/src/data/FeatureFinderMetaboIdent_1_input.mzML", output_ms_data)

# load MS data into MSExperiment()
exp = MSExperiment()
MzMLFile().load(output_ms_data, exp)

# choose one of the following three methods to access the TIC data
# 1) recalculate TIC data with the calculateTIC() function
tic = exp.calculateTIC()
retention_times, intensities = tic.get_peaks()

# 2) get TIC data using list comprehensions
retention_times = [spec.getRT() for spec in exp]
intensities = [sum(spec.get_peaks()[1]) for spec in exp if spec.getMSLevel() == 1]

# 3) get TIC data looping over spectra in MSExperiment()
retention_times = []
intensities = []
for spec in exp:
    if spec.getMSLevel() == 1:
        retention_times.append(spec.getRT())
        intensities.append(sum(spec.get_peaks()[1]))

# plot retention times and intensities and add labels
plt.plot(retention_times, intensities)

plt.title("TIC")
plt.xlabel("time (s)")
plt.ylabel("intensity (cps)")

plt.show()<jupyter_output><empty_output><jupyter_text>Spectrum`MSSpectrum` is the most important container for raw data and peaks. `MSSpectrum` is a container for 1D peak data (container for `Peak1D`). You can directly access these objects, however it is faster to use `get_peaks()` and `set_peaks` functions. These functions use numpy arrays for raw data access. Access spectrum metadata through `SpectrumSettings`.In the example below, let's instantiate the `MSSpectrum` with peaks, sort the peaks according to m/z ratio, and finally let's display peak positions.<jupyter_code>from pyopenms import MSSpectrum

spectrum = MSSpectrum()
mz = range(1500, 500, -100)
# i --> intensity
i = [0 for mass in mz]
spectrum.set_peaks([mz, i])

# Sort peak according to ascending m/z ratio
spectrum.sortByPosition()

# Iterate over spectrum of those peaks
# This process is expensive and should NOT be used in production code
# Direct access to the Peak1D object is expensive and should NOT be done
for p in spectrum:
    # p --> Peak1D object
    print(p.getMZ(), p.getIntensity())

# More efficient peak access can be achieved with get_peaks()
# Unpack tuple to allow 1-to-1 zipping
for mz, i in zip(*spectrum.get_peaks()):
    print(mz, i)

# Again, you can access a peak by its index
print(spectrum[2].getMZ(), spectrum[2].getIntensity())

from pyopenms import (
    MSExperiment,
    MSSpectrum,
    Precursor,
    InstrumentSettings,
    FloatDataArray,
    StringDataArray,
    IonSource,
)

# create spectrum and set properties
spectrum = MSSpectrum()
spectrum.setDriftTime(25)  # 25 ms
spectrum.setRT(205.2)  # 205.2 s
spectrum.setMSLevel(3)  # MS3

# Add peak(s) to spectrum
spectrum.set_peaks(([401.5], [900]))

# create precursor information
p = Precursor()
p.setMZ(600)  # isolation at 600 +/- 1.5 Th
p.setIsolationWindowLowerOffset(1.5)
p.setIsolationWindowUpperOffset(1.5)
p.setActivationEnergy(40)  # 40 eV
p.setCharge(4)  # 4+ ion

# and store precursor in spectrum
spectrum.setPrecursors([p])

# set additional instrument settings (e.g. scan polarity)
IS = InstrumentSettings()
IS.setPolarity(IonSource.Polarity.POSITIVE)
spectrum.setInstrumentSettings(IS)

# get and check scan polarity
polarity = spectrum.getInstrumentSettings().getPolarity()
if polarity == IonSource.Polarity.POSITIVE:
    print("scan polarity: positive")
elif polarity == IonSource.Polarity.NEGATIVE:
    print("scan polarity: negative")

# Optional: additional data arrays / peak annotations
fda = FloatDataArray()
fda.setName("Signal to Noise Array")
fda.push_back(15)
sda = StringDataArray()
sda.setName("Peak annotation")
sda.push_back("y15++")
spectrum.setFloatDataArrays([fda])
spectrum.setStringDataArrays([sda])

# Add spectrum to MSExperiment
exp = MSExperiment()
exp.addSpectrum(spectrum)

# Add second spectrum to the MSExperiment
spectrum2 = MSSpectrum()
spectrum2.set_peaks(([1, 2], [1, 2]))
exp.addSpectrum(spectrum2)

output_test_file = os.path.join(DATA_DIR, "testfile.mzML")
# store spectra in mzML file
MzMLFile().store(output_test_file, exp)<jupyter_output>scan polarity: positive<jupyter_text>What just happened above?- We created a single spectrum and set basic spectrum properties    - drift time    - RT    - MS level    - precursor charge    - isolation window    - activation energy- We also set instrument settings (i.e. polarity of ion source)- We then added actual peaks to the spectrum    - 401.5 m/z with 900 intensity- We then instantiated an MSExperiment object and saved our results to an .mzML fileWe can view this .mzML file in TOPPView, an application that is installed with OpenMS.We can also view our spectrum with matplotlib, as follows:<jupyter_code>import matplotlib.pyplot as plt


def plot_spectrum(spectrum):
    # plot every peak in spectrum and annotate with it's m/z
    for mz, i in zip(*spectrum.get_peaks()):
        plt.plot([mz, mz], [0, i], color="black")
        plt.text(mz, i, str(mz))

    # for the title add RT and Precursor m/z if available
    title = ""
    if spectrum.getRT() >= 0:
        title += "RT: " + str(spectrum.getRT())
    if len(spectrum.getPrecursors()) >= 1:
        title += "   Precursor m/z: " + str(spectrum.getPrecursors()[0].getMZ())

    plt.title(title)
    plt.ylabel("intensity")
    plt.xlabel("m/z")
    plt.ylim(bottom=0)

    plt.show()


# plotting out spectrum that was defined earlier
plot_spectrum(spectrum)<jupyter_output><empty_output><jupyter_text>ChromatogramWe can use the MSChromatogram class, which is similar to the MSSpectrum class, but contains ChromatogramPeak and is derived from ChromatogramSettings<jupyter_code>import numpy as np
from pyopenms import MSChromatogram, Product


def gaussian(x, mu, sig):
    return np.exp(-np.power(x - mu, 2.0) / (2 * np.power(sig, 2.0)))


# Create new chromatogram
chromatogram = MSChromatogram()

# Set raw data (RT and intensity)
rt = range(1500, 500, -100)
i = [gaussian(rtime, 100, 150) for rtime in rt]
chromatogram.set_peaks([rt, i])

# Sort peaks according to ascending retention time
chromatogram.sortByPosition()

# Iterate over chromatogram of those peaks
for p in chromatogram:
    print(p.getRT(), p.getIntensity())

# Here's a more efficient peak access using the get_peaks() method
for rt, i in zip(*chromatogram.get_peaks()):
    print(rt, i)

# Again, let's access a peak by its index
print(chromatogram[2].getRT(), chromatogram[2].getIntensity())

# Let's add metadata to our chromatogram
chromatogram.setNativeID("Trace XIC@405.2")

# Let's add a precursor ion for the chromatogram
p = Precursor()
p.setIsolationWindowLowerOffset(1.5)
p.setIsolationWindowUpperOffset(1.5)
p.setMZ(405.2)  # isolation at 405.2 +/- 1.5 Th
p.setActivationEnergy(40)  # 40 eV
p.setCharge(2)
p.setMetaValue("description", chromatogram.getNativeID())
p.setMetaValue("peptide_sequence", chromatogram.getNativeID())
chromatogram.setPrecursor(p)

# Also store a product ion for the chromatogram (e.g., for SRM)
p = Product()
p.setMZ(603.4)  # Transition from 405.2 --> 603.4
chromatogram.setProduct(p)

# Finally, instantiate a new experiment and store chromatogram as an mzML file
exp = MSExperiment()
exp.addChromatogram(chromatogram)
output_file_test3 = os.path.join(DATA_DIR, "testfile3.mzML")
MzMLFile().store(output_file_test3, exp)

# Now let's visualize our exported file in matplotlib

import matplotlib.pyplot as plt

for chrom in exp.getChromatograms():
    retention_times, intensities = chrom.get_peaks()
    plt.plot(retention_times, intensities, label=chrom.getNativeID())

plt.xlabel("time (s)")
plt.ylabel("intensity (cps)")
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>What we see above shows that the `MSExperiment` class can hold:        - spectra    - chromatograms  Now let's take a look at working with LC-MS/MS ExperimentIn OpenMS, LC-MS/MS injections are represented as so-called peak maps (we can use the `MSExperiment` class here).`MSExperiment` can hold peak maps as well as metadata about the injections.Access to individual spectra is performed through `MSExperiment.getSpectrum` and `MSExperiment.getChromatogram`Let's see an example of this below<jupyter_code># Create an MSExperiment instance with six MSSpectrum instances

from pyopenms import MSExperiment, MSSpectrum, Peak1D

exp = MSExperiment()
for i in range(6):
    spectrum = MSSpectrum()
    spectrum.setRT(i)
    spectrum.setMSLevel(1)
    for mz in range(500, 900, 100):
        peak = Peak1D()
        peak.setMZ(mz + i)
        peak.setIntensity(100 - 25 * abs(i - 2.5))
        spectrum.push_back(peak)
    exp.addSpectrum(spectrum)

# Now iterate over spectra
for spect in exp:
    for peak in spect:
        print(spectrum.getRT(), peak.getMZ(), peak.getIntensity())<jupyter_output>5.0 500.0 37.5
5.0 600.0 37.5
5.0 700.0 37.5
5.0 800.0 37.5
5.0 501.0 62.5
5.0 601.0 62.5
5.0 701.0 62.5
5.0 801.0 62.5
5.0 502.0 87.5
5.0 602.0 87.5
5.0 702.0 87.5
5.0 802.0 87.5
5.0 503.0 87.5
5.0 603.0 87.5
5.0 703.0 87.5
5.0 803.0 87.5
5.0 504.0 62.5
5.0 604.0 62.5
5.0 704.0 62.5
5.0 804.0 62.5
5.0 505.0 37.5
5.0 605.0 37.5
5.0 705.0 37.5
5.0 805.0 37.5<jupyter_text>What we did above is create 6 instances of `MSSpectrum`, populate it with peaks @ 500, 900, and 100 m/z, and append the group to `MSExperiment`.<jupyter_code># We can instead iterate over the spectra using a list comprehension and sum the intensities
# Sum intensity of all spectra between RT 2.0 and 3.0
print(sum([p.getIntensity() for s in exp if s.getRT() >= 2.0 and s.getRT() <= 3.0 for p in s]))

# Again, we can store our experiment as an .mzML file
output_test_data2 = os.path.join(DATA_DIR, "testfile2.mzML")
MzMLFile().store(output_test_data2, exp)<jupyter_output><empty_output><jupyter_text>You can view the above file (`testfile2.mzML`) in TOPPView (3D mode)OR you can view the data directly with Python using `matplotlib`.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors


def plot_spectra_2D(exp, ms_level=1, marker_size=5):
    exp.updateRanges()
    print("collecting peak data...")
    for spec in exp:
        if spec.getMSLevel() == ms_level:
            mz, intensity = spec.get_peaks()
            p = intensity.argsort()  # sort by intensity to plot highest on top
            rt = np.full([mz.shape[0]], spec.getRT(), float)
            plt.scatter(
                rt,
                mz[p],
                c=intensity[p],
                cmap="afmhot_r",
                s=marker_size,
                norm=colors.LogNorm(exp.getMinInt() + 1, exp.getMaxInt()),
            )
    plt.clim(exp.getMinInt() + 1, exp.getMaxInt())
    plt.xlabel("time (s)")
    plt.ylabel("m/z")
    plt.colorbar()
    print("showing plot...")
    plt.show()  # slow for larger data sets


from urllib.request import urlretrieve

gh = "https://raw.githubusercontent.com/OpenMS/pyopenms-extra/master"
output_filepath = os.path.join(DATA_DIR, "test.mzML")
urlretrieve(gh + "/src/data/FeatureFinderMetaboIdent_1_input.mzML", output_filepath)

exp = MSExperiment()
MzMLFile().load(output_filepath, exp)

plot_spectra_2D(exp)<jupyter_output>collecting peak data...
showing plot...<jupyter_text>For larger datasets this will be too slow, as every individual peak is rendered.We can instead use `BilinearInterpolation` to produce an overview image of our spectra.This is useful for brief inspection of your sample for QC<jupyter_code>import numpy as np
import matplotlib.pyplot as plt

from pyopenms import BilinearInterpolation


def plot_spectra_2D_overview(experiment):
    rows = 200.0
    cols = 200.0
    exp.updateRanges()

    bilip = BilinearInterpolation()
    tmp = bilip.getData()
    tmp.resize(int(rows), int(cols), float())
    bilip.setData(tmp)
    bilip.setMapping_0(0.0, exp.getMinRT(), rows - 1, exp.getMaxRT())
    bilip.setMapping_1(0.0, exp.getMinMZ(), cols - 1, exp.getMaxMZ())
    print("collecting peak data...")
    for spec in exp:
        if spec.getMSLevel() == 1:
            mzs, ints = spec.get_peaks()
            rt = spec.getRT()
            for i in range(0, len(mzs)):
                bilip.addValue(rt, mzs[i], ints[i])

    data = np.ndarray(shape=(int(cols), int(rows)), dtype=np.float64)
    for i in range(int(rows)):
        for j in range(int(cols)):
            data[i][j] = bilip.getData().getValue(i, j)

    plt.imshow(np.rot90(data), cmap="gist_heat_r")
    plt.xlabel("retention time (s)")
    plt.ylabel("m/z")
    plt.xticks(
        np.linspace(0, int(rows), 20, dtype=int),
        np.linspace(exp.getMinRT(), exp.getMaxRT(), 20, dtype=int),
    )
    plt.yticks(
        np.linspace(0, int(cols), 20, dtype=int),
        np.linspace(exp.getMinMZ(), exp.getMaxMZ(), 20, dtype=int)[::-1],
    )
    print("showing plot...")
    plt.show()


plot_spectra_2D_overview(exp)<jupyter_output>collecting peak data...
showing plot...<jupyter_text>Inspecting Precursor PurityWhen an MS/MS spectrum is generated, the precursor from the MS1 spectrum is gathered, fragmented, and measured.In reality, the m/z isolation is actually a window - everything within this window is fragmented.It goes to say that there will be contaminants in this precursor isolation window. This can lead to issues in quantification of samples.Let's assess precursor purity and filter spectra that fall below our defined score threshold.<jupyter_code>from pyopenms import MSExperiment, Precursor, PrecursorPurity

gh = "https://raw.githubusercontent.com/OpenMS/pyopenms-extra/master"
ouput_filepath = os.path.join(DATA_DIR, "PrecursorPurity_input.mzML")
urlretrieve(gh + "/src/data/PrecursorPurity_input.mzML", output_filepath)

exp = MSExperiment()
MzMLFile().load(output_filepath, exp)

# for this example, we check which are MS2 spectra and choose one of them
for element in exp:
    print(element.getMSLevel())

# get the precursor information from the MS2 spectrum at index 3
ms2_precursor = Precursor()
ms2_precursor = exp[3].getPrecursors()[0]

# get the previous recorded MS1 spectrum
isMS1 = False
i = 3  # start at the index of the MS2 spectrum
while isMS1 == False:
    if exp[i].getMSLevel() == 1:
        isMS1 = True
    else:
        i -= 1

ms1_spectrum = exp[i]

# calculate the precursor purity in a 10 ppm precursor isolation window
purity_score = PrecursorPurity().computePrecursorPurity(ms1_spectrum, ms2_precursor, 10, True)

print(purity_score.total_intensity)  # 9098343.890625
print(purity_score.target_intensity)  # 7057944.0
print(purity_score.signal_proportion)  # 0.7757394186070014
print(purity_score.target_peak_count)  # 1
print(purity_score.residual_peak_count)  # 4<jupyter_output>1
2
2
2
2
2
1
9098343.890625
7057944.0
0.7757394186070014
1
4<jupyter_text>In the above example, we can say that the signal of the isotopic peaks correspond to roughly 78% of all intensities in the precursor isolation window.  Taking a look at filtering spectraHere we will go through the various filtering methods we can use on our spectra data.<jupyter_code># Load data
from pyopenms import MSExperiment, MzMLFile

gh = "https://raw.githubusercontent.com/OpenMS/pyopenms-extra/master"
output_filepath = os.path.join(DATA_DIR, "test.mzML")
urlretrieve(gh + "/src/data/tiny.mzML", output_filepath)

inp = MSExperiment()
MzMLFile().load(output_filepath, inp)

# Filter spectra by MS level
# Let's filter all MS spectra that IS NOT MS1
filtered = MSExperiment()
for s in inp:
    if s.getMSLevel() > 1:
        filtered.addSpectrum(s)

# We can also use scan numbers to filter spectra we are interested in.
scan_nrs = [0, 2, 5, 7]

filtered = MSExperiment()
for k, s in enumerate(inp):
    if k in scan_nrs:
        filtered.addSpectrum(s)

# If we are only interested in a small m/z window of our fragment
# ion spectra, we can filter our data like so:

mz_start = 6.0
mz_end = 12.0
filtered = MSExperiment()
for s in inp:
    if s.getMSLevel() > 1:
        filtered_mz = []
        filtered_int = []
        for mz, i in zip(*s.get_peaks()):
            if mz > mz_start and mz < mz_end:
                filtered_mz.append(mz)
                filtered_int.append(i)
        s.set_peaks((filtered_mz, filtered_int))
        filtered.addSpectrum(s)

# `filtered` will now only contain spectra with peaks in range [`mz_start`, and `mz_end`]
# In real-world applications, we'd set `mz_start` and `mz_end` to a relevant range<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1106504025842534, 'fraction_numerical': 0.06504025842534111, 'mean_word_length': 4.287280142747897, 'pattern_counts': {'":': 0, '<': 48, '<?xml version=': 0, '>': 57, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '217425', 'license': 'MIT', 'n_tokens_mistral': 8071, 'n_tokens_neox': 7045, 'n_words': 2091, 'path': 'pyopenms-demo.ipynb', 'repo_name': 'irahorecka/roest-ms', 'size': '20741'}