Tutorial¶
1: Chemical calculator¶
import imp
import itertools
import pprint
from lipyd import pprint_namedtuple
pprint.PrettyPrinter = pprint_namedtuple.PrettyPrinter
from lipyd import mass
from lipyd import formula
The mass
module knows exact masses of isotopes, isotopic abundances,
weights, etc. The MassBase
class is able to process chemical
formula, calculate masses and do arithmetics.
mass.get_mass('Na')
22.989769282
mass.MassBase('C2H6') - mass.MassBase('H') + mass.MassBase('OH')
46.04186481376
Expression evaluation with mass.expr
:
mass.expr('C2H6 - H + OH')
46.04186481376
mass.expr('C6H12O6 - water')
162.0528234236
Make a deuterium:
mass.expr('H + n')
2.0164899481400003
Hydrogensulphate ion:
mass.expr('HSO4 + e')
96.96010326786924
Additional attributes can be provided in keyword arguments to carry metadata.
lactic_acid = formula.Formula('CH3CHOHCOOH', name = 'lactic acid')
lactic_acid.attrs.name
'lactic acid'
A galactose:
((2 * formula.Formula('C6H12O6')) - 'H2O').formula
'C12H22O11'
2: Calculations with adducts¶
from lipyd import mz
This is an oleic acid. What is the mass of the [M-H]- adduct?
formula.Formula('C18H34O2').remove_h()
281.24860387047
We have seen a mass and wondering what is the exact mass if it is an [M+NH4]+ adduct:
mz.Mz(854.576482597791).remove_nh4()
836.5426570444603
Calculate the [M+Li]+ adduct for the same molecule:
mz.Mz(836.5426570444603).adduct(mass.MassBase('Li', charge = 1))
843.558111907551
3: Metabolite model¶
Metabolites consist of a core and optionally substituents. Substituents might be formulas or moieties with aliphatic chains.
from lipyd import metabolite
from lipyd import substituent
Make all combinations of halogenated methanes:
halo_methanes = metabolite.AbstractMetabolite(
core = 'C',
subs = [('H', 'F', 'Cl', 'Br', 'I')] * 4
)
Check the first 3:
[(m.formula, m.mass) for m in halo_methanes][:3]
[('C1H4', 16.031300129039998),
('C1F1H3', 34.02187826038),
('C1Cl1H3', 49.99232782678)]
Do the same with all alcohols up to 1-8 carbon count with 0-2 unsaturated bonds:
chain = metabolite.AbstractSubstituent(c = (1, 8), u = (0, 2))
alcohols = metabolite.AbstractMetabolite(subs = (chain, ('OH',)))
[(m.formula, m.mass) for m in alcohols][:3]
[('C1H4O1', 32.02621474924),
('C2H6O1', 46.04186481376),
('C3H8O1', 60.057514878279996)]
Make some ceramides:
# fatty acyls of length 16, 18 or 20 and one or no unsaturation:
fattyacyl = substituent.FattyAcyl(c = (16, 18, 20), u = (0, 1))
lcb = substituent.Sphingosine(c = 18, u = 1)
ceramides = metabolite.AbstractMetabolite(core = 'H', subs = (lcb, fattyacyl), name = 'Cer')
# name, formula, mass, [M+H]+, [M+NH4]+
[(cer.name, cer.formula, cer.mass, cer.add_nh4()) for cer in ceramides]
[('Cer(18:1/16:0)', 'C34H67N1O3', 537.51209502622, 555.5459205795507),
('Cer(18:1/16:1)', 'C34H65N1O3', 535.4964449617, 553.5302705150308),
('Cer(18:1/18:0)', 'C36H71N1O3', 565.54339515526, 583.5772207085907),
('Cer(18:1/18:1)', 'C36H69N1O3', 563.52774509074, 581.5615706440708),
('Cer(18:1/20:0)', 'C38H75N1O3', 593.5746952843, 611.6085208376307),
('Cer(18:1/20:1)', 'C38H73N1O3', 591.55904521978, 609.5928707731108)]
4: Lipid definitions¶
In the lipyd.lipid
module more than 150 lipid varieties are
predefined.
from lipyd import lipid
d_ceramides = lipid.CeramideD(
sph_args = {'c': 18, 'u': (0, 1)},
fa_args = {'c': (14, 24), 'u': (0, 4), 'even': True},
)
[(cer.name, cer.mass) for cer in d_ceramides][:3]
[('Cer(DH18:0/14:0)', 511.49644496170004),
('Cer(DH18:0/14:1)', 509.48079489718003),
('Cer(DH18:0/14:2)', 507.46514483266003)]
5: External databases¶
The lipyd.moldb
module provides access to SwissLipis and
LipidMaps. Automatically downloads and processes the databases which
you can search by various names and identifiers, you can also access
structures as OpenBabel objects, InChI and SMILE strings.
from lipyd import moldb
from lipyd import name
from lipyd.lipproc import *
5.1: SwissLipids¶
swl = moldb.SwissLipids(levels = {'species'})
Indexing SwissLipids -- finished: 100%|██████████| 623M/623M [00:13<00:00, 47.4Mit/s]
Get phosphatidylethanolamines as OpenBabel objects:
swl.reload()
pe = swl.get_hg_obmol('PE')
pe0 = next(pe)
pe0.draw()
[m.title for m in itertools.islice(pe, 0, 3)]
['Phosphatidylethanolamine (28:1)',
'Phosphatidylethanolamine (39:3)',
'Phosphatidylethanolamine (O-29:1)']
5.2: LipidMaps¶
lm = moldb.LipidMaps()
:: Indexed 42556 records from cache/LMSDFDownload12Dec17FinalAll.sdf.
gibberellin = list(lm.get_record('LMPR0104170034', typ = 'mainkey'))[0]
gibberellin['name']
{'COMMON_NAME': 'gibberellin A17',
'SYNONYMS': 'GA17; gibberellin 17',
'PUBCHEM_CID': '5460657',
'CHEBI_ID': '24236',
'INCHI': 'InChI=1S/C20H26O7/c1-10-8-18-9-19(10,27)7-4-11(18)20(16(25)26)6-3-5-17(2,15(23)24)13(20)12(18)14(21)22/h11-13,27H,1,3-9H2,2H3,(H,21,22)(H,23,24)(H,25,26)/t11-,12-,13-,17-,18+,19+,20-/m1/s1'}
LipidMaps too is able to yield OpenBabel objects:
tag = list(lm.get_obmol('TAG(15:0_20:4_20:5)', 'synonym'))[0]
tag.draw()
6: Lipid name parser¶
In order to make the databases computationally useful and use them as a
combined database, we need to process their nomenclature. The
lipyd.name
module is able to recognize dozens of lipid names used in
SwissLipids and LipidMaps.
from lipyd import name
nameproc = name.LipidNameProcessor(database = 'swisslipids', iso = True)
processed_name = nameproc.process(
['Phosphatidylethanolamine (16:0/20:4(5Z,8Z,11Z,14Z))']
)
pprint.pprint(processed_name)
(
Headgroup(main='PE', sub=()),
ChainSummary(
c = 36,
u = 4,
typ = ('FA', 'FA'),
attr = (
ChainAttr(sph='', ether=False, oh=()),
ChainAttr(sph='', ether=False, oh=())
),
iso = None
),
(
Chain(
c = 16,
u = 0,
typ = 'FA',
attr = ChainAttr(sph='', ether=False, oh=()),
iso = ()
),
Chain(
c = 20,
u = 4,
typ = 'FA',
attr = ChainAttr(sph='', ether=False, oh=()),
iso = ('5Z', '8Z', '11Z', '14Z')
)
)
)
It understands even greek names:
nameproc.process(['eicosapentaenoate'])
(Headgroup(main='FA', sub=()),
ChainSummary(c=20, u=5, typ=('FA',), attr=(ChainAttr(sph='', ether=False, oh=()),), iso=None),
[Chain(c=20, u=5, typ='FA', attr=ChainAttr(sph='', ether=False, oh=()), iso=())])
7: Combined molecule database¶
db = moldb.MoleculeDatabaseAggregator()
Indexing SwissLipids -- finished: 100%|██████████| 623M/623M [00:13<00:00, 16.0Mit/s]
:: Indexed 42556 records from cache/LMSDFDownload12Dec17FinalAll.sdf.
Generating metabolites -- finished: 100%|██████████| 44.0/44.0 [00:09<00:00, 4.50it/s]
Generating metabolites -- finished: 100%|██████████| 18.0/18.0 [00:05<00:00, 4.10it/s]
Generating metabolites -- finished: 100%|██████████| 106/106 [00:09<00:00, 11.4it/s]
Generating metabolites -- finished: 100%|██████████| 1.00/1.00 [00:00<00:00, 11.9it/s]
Generating metabolites -- finished: 100%|██████████| 1.00/1.00 [00:00<00:00, 48.9it/s]
Either exact masses or adducts can be searched in the database by
lookup
and adduct_lookup
methods, respectively.
result = db.adduct_lookup(757.549011, ionmode = 'pos')
Take a closer look at one of the resulted records:
pprint.pprint(result['[M+NH4]+'][1][2])
LipidRecord(
lab = LipidLabel(db_id=None, db='lipyd.lipid', names=('PC(33:4)',)),
hg = Headgroup(main='PC', sub=()),
chainsum = ChainSummary(
c = 33,
u = 4,
typ = ('FA', 'FA'),
attr = (
ChainAttr(sph='', ether=False, oh=()),
ChainAttr(sph='', ether=False, oh=())
),
iso = None
),
chains = ()
)
The exact masses and errors for all hits are also provided. Errors in ppm:
result['[M+NH4]+'][2]
array([-1.69750819e-02, -1.69750819e-02, -2.69161082e-02, -2.69161082e-02,
-2.69161082e-02, -2.69161082e-02, -3.04974547e-02, -3.04974547e-02,
-3.04974547e-02, -3.04974547e-02, -3.04974547e-02, -3.04974547e-02,
-3.04974547e-02, -3.04974547e-02, -3.04974547e-02, -3.04974547e-02,
-3.04974547e-02, -3.04974547e-02, -3.04974547e-02, -3.04974547e-02,
-3.04974547e-02, -3.04974547e-02, -3.04974547e-02, -3.04974547e-02,
-3.04974547e-02, -3.04974547e-02, -3.04974547e-02, -3.04974547e-02,
-3.04974547e-02, -3.04974547e-02, -3.04974547e-02, -3.23033863e+00,
-3.23033863e+00, -3.23033863e+00, -1.11723392e+01, 1.73995818e+01,
1.73995818e+01, 1.73995818e+01])
Repeat the lookup with lower tolerance, and the items with high ppm disappear:
result = db.adduct_lookup(757.549011, ionmode = 'pos', tolerance = 5)
result['[M+NH4]+'][2]
array([-0.01697508, -0.01697508, -0.02691611, -0.02691611, -0.02691611,
-0.02691611, -0.03049745, -0.03049745, -0.03049745, -0.03049745,
-0.03049745, -0.03049745, -0.03049745, -0.03049745, -0.03049745,
-0.03049745, -0.03049745, -0.03049745, -0.03049745, -0.03049745,
-0.03049745, -0.03049745, -0.03049745, -0.03049745, -0.03049745,
-0.03049745, -0.03049745, -0.03049745, -0.03049745, -0.03049745,
-0.03049745, -3.23033863, -3.23033863, -3.23033863])
8: MS2 fragment definitions¶
The fragment database provided by lipyd.fragment
and
lipyd.fragdb
modules works similar way as lipyd.lipid
and
lipyd.moldb
. lipyd.fragment
contains near 100 predefined
aliphatic chain derived fragments. In addition 140 headgroup derived
fragments are included like for example 184 for choline.
from lipyd import fragment
from lipyd import fragdb
As an example take a look at a [Sph-NH2-OH]- fragment:
sphfrag = fragment.Sph_mNH2_mOH(c = 18, u = 1)
At this fragment type the constraints tell us which lipid varieties this fragment can be observed. In this case dCer and DHCer.
sphfrag.constraints
(FragConstraint(hg='Cer', family=None, sub=None, sph='d', oh=0, chaintype='Sph'),
FragConstraint(hg='Cer', family=None, sub=None, sph='DH', oh=0, chaintype='Sph'))
list(sphfrag)[0].charge, list(sphfrag)[0].mass
(-1, 263.2380392001693)
9: MS2 fragment database¶
Look up a negative mode fragment m/z in the database. It results an array with mass, fragment name, fragment type, aliphatic chain type, carbon count, unsaturation and charge in each row. At neutral losses the charge is 0.
fragdb.lookup_neg(283.26)
array([[283.2642539494093, '[FA(18:0)-H]-', 'FA-H', 'FA', 18, 0, -1],
[283.2642539494093, '[Sph(20:0)-C2H4-NH2-2H]-', 'Sph-C2H4-NH2-2H',
'Sph', 20, 0, -1]], dtype=object)
Now let’s annotate an MS2 scan with possible fragment identifications.
To do this we open an example MGF file included in the module. The
lipyd.mgf
module serves MS2 scans from MGF files on demand. Btw the
lipyd.settings
module gives easy access for and control over near
100 customizable parameters.
from lipyd import mgf
from lipyd import settings
mgffile = settings.get('mgf_example')
mgfreader = mgf.MgfReader(mgffile)
precursor = 590.45536 # this is a Cer-1P
idx, rtdiff = mgfreader.lookup_scan_ids(precursor)
We found the following scans for precursor 590.455:
idx
array([1941, 1929, 427, 423, 589, 645, 642, 308, 481, 478, 586,
368, 696, 535, 532, 755, 752, 700, 721])
Select a scan from the ones above and annotate its fragments:
scan = mgfreader.scan_by_id(1941)
annot = fragdb.FragmentAnnotator(
mzs = scan[:,0],
ionmode = 'pos',
precursor = precursor
)
One example of the annotations. This fragment ranks 25 by intensity.
pprint.pprint(list(annot)[24])
(
FragmentAnnotation(
mz = 228.23219101229077,
name = '[Sph(14:0)-H2O+H]+',
fragtype = 'Sph-H2O+H',
chaintype = 'Sph',
c = 14,
u = 0,
charge = 1
),
FragmentAnnotation(
mz = 228.23219101229077,
name = '[FA(14:0)+NH2-O]+',
fragtype = 'FA+NH2-O',
chaintype = 'FA',
c = 14,
u = 0,
charge = 1
),
FragmentAnnotation(
mz = 228.23219101229077,
name = '[Sph(14:1)-O+H]+',
fragtype = 'Sph-O+H',
chaintype = 'Sph',
c = 14,
u = 1,
charge = 1
)
)
10: MS2 spectrum analysis¶
The lipyd.ms2.Scan
class is able to perform the entire
identification workflow. By an alternative constructor method it can be
initialized by providing and MGF file and scan ID.
from lipyd import ms2
mgfname = settings.get('mgf_pos_examples')
scan_id = 3626
scan = ms2.Scan.from_mgf(mgfname, scan_id, 'pos')
If not provided the Scan
instance performs the database lookup of
the precursor ion. Here are the results:
pprint.pprint(scan.ms1_records['[M+H]+'][1][0])
LipidRecord(
lab = LipidLabel(db_id=None, db='lipyd.lipid', names=('Hex2-Cer(t42:2)',)),
hg = Headgroup(main='Cer', sub=('Hex2',)),
chainsum = ChainSummary(
c = 42,
u = 2,
typ = ('Sph', 'FA'),
attr = (
ChainAttr(sph='t', ether=False, oh=()),
ChainAttr(sph='', ether=False, oh=())
),
iso = None
),
chains = ()
)
The identify
method attempts to confirm each of the records by
analysing the MS2 spectrum.
identity = scan.identify()
The results are grouped by lipid species and come with a score. Hex2-Cer(t42:2) got a score of 45, which is the highest at this scan:
pprint.pprint(identity['Hex2-Cer(t42:2)'][0])
MS2Identity(
score = 45,
hg = Headgroup(main='Cer', sub=('Hex2',)),
chainsum = ChainSummary(
c = 42,
u = 2,
typ = ('Sph', 'FA'),
attr = (
ChainAttr(sph='t', ether=False, oh=()),
ChainAttr(sph='', ether=False, oh=())
),
iso = None
),
chains = (
Chain(
c = 18,
u = 1,
typ = 'Sph',
attr = ChainAttr(sph='t', ether=False, oh=()),
iso = ()
),
Chain(
c = 24,
u = 1,
typ = 'FA',
attr = ChainAttr(sph='', ether=False, oh=()),
iso = ()
)
),
details = ChainIdentificationDetails(rank = (0, None), i = (1.0, None), fragtype = ('Sph-2xH2O-H', None))
)
At the same time there were attempts to confirm for example Hex-Cer(d53:9-2OH) but it resulted a score of 0.
pprint.pprint(identity['Hex-Cer(d53:9-2OH)'][0])
MS2Identity(
score = 0,
hg = Headgroup(main='Cer', sub=('Hex',)),
chainsum = ChainSummary(
c = 53,
u = 9,
typ = ('Sph', 'FAOH'),
attr = (
ChainAttr(sph='d', ether=False, oh=()),
ChainAttr(sph='', ether=False, oh=('2OH',))
),
iso = None
),
chains = (
Chain(
c = 17,
u = 1,
typ = 'Sph',
attr = ChainAttr(sph='d', ether=False, oh=()),
iso = ()
),
Chain(
c = 36,
u = 8,
typ = 'FAOH',
attr = ChainAttr(sph='', ether=False, oh=('2OH',)),
iso = ()
)
),
details = ChainIdentificationDetails(
rank = (1, None),
i = (0.18172325900094663, None),
fragtype = ('Sph-2xH2O+H', None)
)
)
Let’s see one more example.
mgfname = settings.get('mgf_neg_examples')
scan_id = 2516
scan = ms2.Scan.from_mgf(mgfname, scan_id, 'neg')
identity = scan.identify()
We see that this is a PI(34:1) with score 11 and both fatty acyl
fragments are confirmed by [FA-H]- ions (see the
ChainIdentificationDetails
object). These fragments are the 1st and
2nd most abundant with relative intensities of 100% and 99%.
pprint.pprint(identity['PI(34:1)'][0])
MS2Identity(
score = 11.0,
hg = Headgroup(main='PI', sub=()),
chainsum = ChainSummary(
c = 34,
u = 1,
typ = ('FA', 'FA'),
attr = (
ChainAttr(sph='', ether=False, oh=()),
ChainAttr(sph='', ether=False, oh=())
),
iso = None
),
chains = (
Chain(
c = 18,
u = 1,
typ = 'FA',
attr = ChainAttr(sph='', ether=False, oh=()),
iso = ()
),
Chain(
c = 16,
u = 0,
typ = 'FA',
attr = ChainAttr(sph='', ether=False, oh=()),
iso = ()
)
),
details = ChainIdentificationDetails(rank = (0, 1), i = (1.0, 0.9897746748655405), fragtype = ('FA-H', 'FA-H'))
)