PMD based Reactomics tried to evaluate untargeted HRMS profiles at reaction level. Reactomics is the extension of structure/reaction directed analysis for HRMS data and PMDs are treated as relationship unit for further discussion as compounds for metabolomics. Reaction level evaluation could be found in PMD network analysis, Source appointment and biomarker reaction discovery.
PMD network analysis
In untargeted metabolites profiles from HRMS, two ions or peaks could be treated together as long as they could be connected by relationship. Regular untargeted workflow prefers intensity correlation between compounds. However, PMD could also be the chemical bridge between two compounds or ions. For example, oxidation would add oxygen atom to the parent compound and introduce a PMD 15.995 Da. Meanwhile, One peak or compound could be involved in multiple reactions. In this case, we could build PMD network for certain ion or compound.
However, we need to the PMDs list to build such network for either one compound or one sample. One way is using the high frequency PMDs list from previous reported reactions and another way is using the high frequency PMDs within certain data set such as KEGG or HMDB. The former way try to focus on known reactions such as Phase I reactions for exogenous compounds while the latter way is useful to explorer new reactions or unknown reaction patterns within the data set. The latter way is actually the structure/reaction directed analysis.
PMD network analysis is the analysis to check or explorer the PMD relationship among co-existed ions from one sample or multiple samples. The edge between ions in the network means certain PMD relationship are valid in this data set. Meanwhile, the intensity correlation between paired ions could also be considered to connect the vertices. Such network is build based on local recursive search for all possible connections.
PMD network analysis is useful to screen metabolites. Regular metabolites discovery method try to predict metabolites’ MS2 spectra and then match the data. In PMD network analysis, metabolites are predicted by high frequency PMD or preferred PMD list within the MS1 data and such prediction could be extended to the metabolites of metabolites. Such PMD search will stop when no new metabolites could be connected to the network. Such searching method is much easier and quick to get the overview of metabolites networks. Identification could follow the discovery from MS1 data.
If you have a specific compound and want to check the metabolites of
certain PMD, you could use getchain to extract the network
of that compounds
library(igraph)
#>
#> Attaching package: 'igraph'
#> The following objects are masked from 'package:stats':
#>
#> decompose, spectrum
#> The following object is masked from 'package:base':
#>
#> union
# check metabolites of C18H39NO
# Use common PMDs for biological reactions
chain <- getchain(spmeinvivo,diff = c(2.02,14.02,15.99,58.04,13.98),mass = 286.3101,digits = 2,corcutoff = 0)
# show as network
net <- graph_from_data_frame(chain$sdac,directed = F)
pal <- grDevices::rainbow(5)
plot(net,vertex.label=round(as.numeric(V(net)$name),2),vertex.size =5,edge.width = 3,edge.color = pal[as.numeric(as.factor(E(net)$diff2))],vertex.label.dist=1,vertex.color=ifelse(round(as.numeric(V(net)$name),4) %in% 286.3101,'red','black'), main = 'PMD network')
legend("topright",bty = "n",
legend=unique(E(net)$diff2),
fill=unique(pal[as.numeric(as.factor(E(net)$diff2))]), border=NA,horiz = F)
# Consider the correlation coefficient cutoff
chain <- getchain(spmeinvivo,diff = c(2.02,14.02,15.99,58.04,13.98),mass = 286.3101,digits = 2,corcutoff = 0.6)
# show as network
net <- graph_from_data_frame(chain$sdac,directed = F)
pal <- grDevices::rainbow(5)
plot(net,vertex.label=round(as.numeric(V(net)$name),2),vertex.size =5,edge.width = 3,edge.color = pal[as.numeric(as.factor(E(net)$diff2))],vertex.label.dist=1,vertex.color=ifelse(round(as.numeric(V(net)$name),4) %in% 286.3101,'red','black'), main = 'PMD network')
legend("topright",bty = "n",
legend=unique(E(net)$diff2),
fill=unique(pal[as.numeric(as.factor(E(net)$diff2))]), border=NA,horiz = F)
Here only three PMD relationship could be found for C18H39NO. The duplicate edges between two vertices or self-loop edges mean isomer related PMD reactions. If we consider the correlation, the network would be trimmed. Since reaction might not always involve correlation, PMD network analysis could found more potential metabolites.
If you want to see all the independent peaks’ high frequency PMDs as networks for certain sample, the following code will help. This part will use the high frequency PMDs cutoff 12 from the data to build the networks for all the independent peaks.
std <- globalstd(spmeinvivo,sda = F)
#> 75 retention time cluster found.
#> 369 paired masses found
#> 5 unique within RT clusters high frequency PMD(s) used for further investigation.
#> The unique within RT clusters high frequency PMD(s) is(are) 28.03 21.98 44.03 17.03 18.01.
#> 719 isotopologue(s) related paired mass found.
#> 492 multi-charger(s) related paired mass found.
#> 8 retention group(s) have single peaks. 14 23 32 33 54 55 56 75
#> 11 group(s) with multiple peaks while no isotope/paired relationship 4 5 7 8 11 41 42 49 68 72 73
#> 9 group(s) with multiple peaks with isotope without paired relationship 2 9 22 26 52 62 64 66 70
#> 4 group(s) with paired relationship without isotope 1 10 15 18
#> 43 group(s) with paired relationship and isotope 3 6 12 13 16 17 19 20 21 24 25 27 28 29 30 31 34 35 36 37 38 39 40 43 44 45 46 47 48 50 51 53 57 58 59 60 61 63 65 67 69 71 74
#> 291 std mass found.
sda <- getsda(std,freqcutoff = 12)
#> 8 groups were found as high frequency PMD group.
#> 0 was found as high frequency PMD.
#> 2.02 was found as high frequency PMD.
#> 13.98 was found as high frequency PMD.
#> 14.05 was found as high frequency PMD.
#> 15.99 was found as high frequency PMD.
#> 28.03 was found as high frequency PMD.
#> 30.05 was found as high frequency PMD.
#> 126.14 was found as high frequency PMD.
df <- sda$sda
net <- graph_from_data_frame(df,directed = F)
pal <- grDevices::rainbow(length(unique(E(net)$diff2)))
plot(net,vertex.label=NA,vertex.size = 5,edge.width = 3,edge.color = pal[as.numeric(as.factor(E(net)$diff2))],main = 'PMD network')
legend("topright",bty = "n",
legend=unique(E(net)$diff2),
fill=unique(pal[as.numeric(as.factor(E(net)$diff2))]), border=NA,horiz = F)
Here we could find clusters of metabolites. We could detect such network community structure.
# network community structure detection
ceb <- cluster_edge_betweenness(net,weights = abs(E(net)$cor), directed = F)
#> Warning in cluster_edge_betweenness(net, weights = abs(E(net)$cor), directed =
#> F): At vendor/cigraph/src/community/edge_betweenness.c:498 : Membership vector
#> will be selected based on the highest modularity score.
plot(ceb, net,vertex.label=NA,vertex.size = 5,edge.width = 3,) 
# output membership
head(cbind(ceb$membership,ceb$names))
#> [,1] [,2]
#> [1,] "1" "146.118268664153"
#> [2,] "1" "162.112813715075"
#> [3,] "2" "284.295501708984"
#> [4,] "2" "284.295597233849"
#> [5,] "1" "192.160430908203"
#> [6,] "3" "326.342517752117"Such network could also be build on correlation directed analysis which use correlation between paired peaks to build network.
cbp <- enviGCMS::getfilter(std,rowindex = std$stdmassindex)
cda <- getcda(cbp)
df <- cda$cda
# filter based on retention time differences larger than 2 mins
df <- df[df$diffrt>120,]
netc <- graph_from_data_frame(df,directed = F)
plot(netc,vertex.label=NA,vertex.size = 5,edge.width = 3,main = 'Correlation network')
As shown above, correlation network without PMD might merge into one big network, which lose the details of chemical reactions.
Shiny application
The PMD network for certain compound could be generated by run
runPMDnet().
Source appointment
Peaks from samples could be from endogenous compounds or exogenous compounds. However, it’s hard to tell for untargeted analysis. In terms of PMD, if one peak belongs to a high frequency PMD network, it means a relatively high activity. If such sample belongs to a biological specimen, it might be endogenous compound. If a peak show no PMD network with other peaks, the biological system might not have enzyme to make reaction happen. Exogenous compounds will show a lower degree since they are xenobiotics. Since most of the peaks will show a low degree, the median of the degree could be used as cutoff. Then we could make source appointment if the assumption is hold.
In this case, we will have 31 endogenous compounds while 51 exogenous compounds. When you find a peak show differences between groups, you could check the degree to infer its sources.
Another parameter would be the average network distances. Endogenous compounds could form a larger network with long average network distances while exogenous compounds will connected to network with small average network distances. Such parameter could be used to determine the source of unknown compound by checking the average network distances of the compounds PMD network.
Be careful, one compound could be endogenous for one sample while exogenous for another sample. In this case, PMD network would give hints on the sources based on the context of the samples.
Biomarker reaction
Biomarker always means biomarker compounds. However, if we could
quantify the reaction relationship, we could use biomarker reaction to
trace certain biological process. You could use getreact to
extract the ion pairs shared the same PMD and intensity ratio RSD% lower
than certain cutoff. Then the sum of the intensity of all PMD pairs’
ions could be used to compare the reaction level changes among
samples.
pmd <- getreact(spmeinvivo,pmd=15.99)
# show the ions with the same PMD
head(pmd$pmd)
#> ms1 ms2 diff rt1 rt2 diffrt rtg1 rtg2 rtgdiff
#> 109378 174.1489 158.1546 15.99434 337.418 470.364 132.946 44 31 13
#> 114898 175.1481 159.1575 15.99055 614.413 470.364 144.049 17 31 14
#> 212776 209.1552 193.1632 15.99202 611.412 337.633 273.779 17 44 27
#> 219286 211.1695 195.1745 15.99491 614.625 594.485 20.140 17 46 29
#> 255283 228.1973 212.2025 15.99473 453.157 639.100 185.943 45 13 32
#> 290503 244.1921 228.1973 15.99486 614.842 453.157 161.685 17 45 28
#> cor diff2 r rh rl
#> 109378 0.9247231 15.99 25.01946 48.53964 28.03499
#> 114898 0.9164059 15.99 18.12926 35.13135 25.48822
#> 212776 0.7567381 15.99 29.09434 38.95531 41.61674
#> 219286 0.7378409 15.99 26.62548 33.46165 61.56227
#> 255283 0.9749689 15.99 6.76019 30.01754 26.76093
#> 290503 0.9670674 15.99 17.96722 42.20945 30.01754
# show the corresponding quantitative PMD data across samples, each row show the sum of intensity of paired masses qualified for stable mass pairs
head(pmd$pmddata)
#> 1405_Fish1_F1 1405_Fish1_F2 1405_Fish1_F3 1405_Fish2_F1 1405_Fish2_F2
#> 174.1/337 18327.863 11388.643 15368.711 9068.557 19759.098
#> 175.1/614 6746.288 4237.569 6581.992 3085.251 6108.837
#> 209.2/611 30644.459 19496.458 33519.253 11647.966 31169.817
#> 211.2/615 5958.334 4270.901 6285.724 2757.028 6065.346
#> 228.2/453 29245.928 21199.245 38668.407 16822.114 36020.058
#> 244.2/615 41582.379 30025.253 57302.819 20846.835 51716.437
#> 1405_Fish2_F3 1405_Fish3_F1 1405_Fish3_F2 1405_Fish3_F3
#> 174.1/337 12995.786 16878.397 18647.770 25557.292
#> 175.1/614 3420.075 6435.200 7535.579 8806.608
#> 209.2/611 14064.021 29421.515 38893.482 43313.165
#> 211.2/615 3041.366 5809.946 8795.206 8089.021
#> 228.2/453 19484.297 32281.370 36797.052 39304.713
#> 244.2/615 24315.314 45535.922 59156.064 61720.089If your data don’t have retention time, reaction level change can also be checked.
spmeinvivo$rt <- NULL
pmd <- getreact(spmeinvivo,pmd=15.99)
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
#> Warning in as.numeric(ratio1)/as.numeric(ratio2): longer object length is not a
#> multiple of shorter object length
# show the ions with the same PMD
head(pmd$pmd)
#> ms1 ms2 diff cor diff2 r rh rl
#> 109378 174.1489 158.1546 15.99434 0.9247231 15.99 25.01946 48.53964 28.03499
#> 114898 175.1481 159.1575 15.99055 0.9164059 15.99 18.12926 35.13135 25.48822
#> 212776 209.1552 193.1632 15.99202 0.7567381 15.99 29.09434 38.95531 41.61674
#> 219286 211.1695 195.1745 15.99491 0.7378409 15.99 26.62548 33.46165 61.56227
#> 255283 228.1973 212.2025 15.99473 0.9749689 15.99 6.76019 30.01754 26.76093
#> 290503 244.1921 228.1973 15.99486 0.9670674 15.99 17.96722 42.20945 30.01754
# show the corresponding quantitative PMD data across samples, each row show the sum of intensity of paired masses qualified for stable mass pairs
head(pmd$pmddata)
#> 1405_Fish1_F1 1405_Fish1_F2 1405_Fish1_F3 1405_Fish2_F1 1405_Fish2_F2
#> 174.1/337 18327.863 11388.643 15368.711 9068.557 19759.098
#> 175.1/614 6746.288 4237.569 6581.992 3085.251 6108.837
#> 209.2/611 30644.459 19496.458 33519.253 11647.966 31169.817
#> 211.2/615 5958.334 4270.901 6285.724 2757.028 6065.346
#> 228.2/453 29245.928 21199.245 38668.407 16822.114 36020.058
#> 244.2/615 41582.379 30025.253 57302.819 20846.835 51716.437
#> 1405_Fish2_F3 1405_Fish3_F1 1405_Fish3_F2 1405_Fish3_F3
#> 174.1/337 12995.786 16878.397 18647.770 25557.292
#> 175.1/614 3420.075 6435.200 7535.579 8806.608
#> 209.2/611 14064.021 29421.515 38893.482 43313.165
#> 211.2/615 3041.366 5809.946 8795.206 8089.021
#> 228.2/453 19484.297 32281.370 36797.052 39304.713
#> 244.2/615 24315.314 45535.922 59156.064 61720.089Now we have two methods to compute the quantitative PMD responses and user should select method depending on research purposes. ‘static’ will only consider the stable mass pairs across samples and such reactions will be limited by the enzyme or other factors than substrates. ‘dynamic’ will consider the unstable paired masses by normalization the relatively unstable peak with stable peak between paired masses and such reactions will be limited by one or both peaks in the paired masses.
data("spmeinvivo")
pmd <- getreact(spmeinvivo,pmd=15.99,method = 'dynamic')
# show the ions with the same PMD
head(pmd$pmd)
#> ms1 ms2 diff rt1 rt2 diffrt rtg1 rtg2 rtgdiff
#> 43332 147.9900 132.0050 15.98503 145.0880 49.4910 95.5970 20 4 16
#> 77336 162.1128 146.1183 15.99455 256.1585 288.1300 31.9715 5 23 18
#> 103842 172.1705 156.1777 15.99279 478.9360 405.3890 73.5470 43 30 13
#> 121772 177.0914 161.0967 15.99477 574.3390 490.7225 83.6165 39 37 2
#> 140886 183.0809 167.0901 15.99080 533.7950 511.2940 22.5010 48 12 36
#> 143601 184.9858 168.9988 15.98695 85.4930 135.4660 49.9730 2 41 39
#> cor diff2 r rh rl
#> 43332 0.41231393 15.99 34.86972 27.176089 39.84553
#> 77336 0.95185232 15.99 30.42748 89.272157 84.24422
#> 103842 0.64609844 15.99 56.74081 50.385313 42.86262
#> 121772 -0.07743575 15.99 100.61906 8.080905 43.90465
#> 140886 0.54183684 15.99 38.29089 16.749724 42.53298
#> 143601 -0.36694067 15.99 46.87352 18.159277 40.61084
# show the corresponding quantitative PMD data across samples, each row show the sum of intensity of paired masses qualified for stable mass pairs
head(pmd$pmddata)
#> 1405_Fish1_F1 1405_Fish1_F2 1405_Fish1_F3 1405_Fish2_F1 1405_Fish2_F2
#> 1 0.01999437 0.016242675 0.01880472 0.006371798 0.024304585
#> 2 1.53528194 0.962295635 1.10219074 0.982190726 1.604363528
#> 3 0.46404167 1.194387570 1.06934597 0.250054581 0.509676599
#> 4 1.20386141 1.669736528 1.32468366 1.213218611 1.083057203
#> 5 0.01461837 0.006770954 0.01055026 0.009490853 0.008958319
#> 6 9.26463385 12.763610562 3.97502609 14.657795223 4.041131973
#> 1405_Fish2_F3 1405_Fish3_F1 1405_Fish3_F2 1405_Fish3_F3
#> 1 0.01309609 0.011523649 0.011596937 0.017668097
#> 2 1.43917818 1.064822325 1.846236771 2.609787311
#> 3 1.01508469 0.406659176 0.874861257 0.621074214
#> 4 2.05169832 7.765315730 0.832736564 2.117740805
#> 5 0.01549754 0.007511291 0.005490916 0.005883117
#> 6 9.12757512 7.779539564 7.970752923 4.004298286You can also output the quantitative results of all high frequency PMDs existing in the data.
data("spmeinvivo")
# remove redundant peaks
list <- globalstd(spmeinvivo,sda = T)
#> 75 retention time cluster found.
#> 369 paired masses found
#> 5 unique within RT clusters high frequency PMD(s) used for further investigation.
#> The unique within RT clusters high frequency PMD(s) is(are) 28.03 21.98 44.03 17.03 18.01.
#> 719 isotopologue(s) related paired mass found.
#> 492 multi-charger(s) related paired mass found.
#> 8 retention group(s) have single peaks. 14 23 32 33 54 55 56 75
#> 11 group(s) with multiple peaks while no isotope/paired relationship 4 5 7 8 11 41 42 49 68 72 73
#> 9 group(s) with multiple peaks with isotope without paired relationship 2 9 22 26 52 62 64 66 70
#> 4 group(s) with paired relationship without isotope 1 10 15 18
#> 43 group(s) with paired relationship and isotope 3 6 12 13 16 17 19 20 21 24 25 27 28 29 30 31 34 35 36 37 38 39 40 43 44 45 46 47 48 50 51 53 57 58 59 60 61 63 65 67 69 71 74
#> 291 std mass found.
#> PMD frequency cutoff is 6 by PMD network analysis with largest network average distance 6.67 .
#> 53 groups were found as high frequency PMD group.
#> 0 was found as high frequency PMD.
#> 1.98 was found as high frequency PMD.
#> 2.01 was found as high frequency PMD.
#> 2.02 was found as high frequency PMD.
#> 6.97 was found as high frequency PMD.
#> 11.96 was found as high frequency PMD.
#> 12 was found as high frequency PMD.
#> 13.98 was found as high frequency PMD.
#> 14.02 was found as high frequency PMD.
#> 14.05 was found as high frequency PMD.
#> 15.99 was found as high frequency PMD.
#> 16.03 was found as high frequency PMD.
#> 19.04 was found as high frequency PMD.
#> 28.03 was found as high frequency PMD.
#> 30.05 was found as high frequency PMD.
#> 31.99 was found as high frequency PMD.
#> 33.02 was found as high frequency PMD.
#> 37.02 was found as high frequency PMD.
#> 42.05 was found as high frequency PMD.
#> 48.04 was found as high frequency PMD.
#> 48.98 was found as high frequency PMD.
#> 49.02 was found as high frequency PMD.
#> 54.05 was found as high frequency PMD.
#> 56.06 was found as high frequency PMD.
#> 56.1 was found as high frequency PMD.
#> 58.04 was found as high frequency PMD.
#> 58.08 was found as high frequency PMD.
#> 58.11 was found as high frequency PMD.
#> 63.96 was found as high frequency PMD.
#> 66.05 was found as high frequency PMD.
#> 68.06 was found as high frequency PMD.
#> 70.04 was found as high frequency PMD.
#> 70.08 was found as high frequency PMD.
#> 74.02 was found as high frequency PMD.
#> 80.03 was found as high frequency PMD.
#> 82.08 was found as high frequency PMD.
#> 88.05 was found as high frequency PMD.
#> 91.1 was found as high frequency PMD.
#> 93.12 was found as high frequency PMD.
#> 94.1 was found as high frequency PMD.
#> 96.09 was found as high frequency PMD.
#> 101.05 was found as high frequency PMD.
#> 108.13 was found as high frequency PMD.
#> 110.11 was found as high frequency PMD.
#> 112.16 was found as high frequency PMD.
#> 116.08 was found as high frequency PMD.
#> 122.15 was found as high frequency PMD.
#> 124.16 was found as high frequency PMD.
#> 126.14 was found as high frequency PMD.
#> 144.18 was found as high frequency PMD.
#> 148.04 was found as high frequency PMD.
#> 150.2 was found as high frequency PMD.
#> 173.18 was found as high frequency PMD.
newlist <- enviGCMS::getfilter(list,rowindex = list$stdmassindex)
# get high frequency pmd
hfpmd <- unique(newlist$sda$diff2)
# generate quantitative results
pmd <- getreact(newlist,pmd=hfpmd)
# output the kegg pmd in the data
table(pmd$pmd$diff2)
#>
#> 0 1.98 2.01 2.02 6.97 11.96 12 13.98 14.02 14.05 15.99
#> 1 3 2 6 1 1 3 4 3 6 3
#> 16.03 19.04 28.03 30.05 31.99 33.02 42.05 48.04 48.98 49.02 54.05
#> 3 2 5 1 2 1 2 1 1 3 1
#> 56.06 58.04 58.08 63.96 66.05 68.06 70.04 70.08 74.02 80.03 82.08
#> 1 6 1 2 2 1 1 2 7 3 2
#> 88.05 91.1 93.12 94.1 96.09 101.05 108.13 110.11 112.16 116.08 122.15
#> 3 2 3 2 1 1 2 3 4 6 2
#> 124.16 126.14 148.04
#> 3 6 2
# output quantitative result for each PMD
head(pmd$pmddata)
#> 1405_Fish1_F1 1405_Fish1_F2 1405_Fish1_F3 1405_Fish2_F1 1405_Fish2_F2
#> 200.2/522 76781.54 40625.84 77922.53 28510.920 79221.24
#> 228.2/453 38231.65 26164.21 44139.68 21243.456 45048.10
#> 240.2/639 31881.67 19485.44 31170.56 13105.650 33620.20
#> 268.3/639 39620.62 22678.74 40315.71 19270.755 39626.76
#> 254.2/639 27475.32 16563.06 26416.10 12568.091 26500.43
#> 261.2/557 24604.42 14415.48 21970.93 9424.969 25368.75
#> 1405_Fish2_F3 1405_Fish3_F1 1405_Fish3_F2 1405_Fish3_F3
#> 200.2/522 40736.98 68844.35 83690.73 91978.11
#> 228.2/453 26938.08 40441.85 45445.05 52223.88
#> 240.2/639 20422.46 29060.96 32886.11 39062.02
#> 268.3/639 24927.91 36500.59 41108.70 48308.70
#> 254.2/639 15106.66 27064.28 30026.87 32384.46
#> 261.2/557 11698.38 21809.08 37162.69 37774.73
# output quantitative result for unique PMD
upmd <- aggregate(pmd$pmddata, by=list(pmd$pmd$diff2),sum)
# column for samples and row for unique PMD
head(upmd)
#> Group.1 1405_Fish1_F1 1405_Fish1_F2 1405_Fish1_F3 1405_Fish2_F1 1405_Fish2_F2
#> 1 0.00 72539.75 46921.79 83489.08 28793.02 80299.40
#> 2 1.98 112810.76 66245.18 146623.69 47083.72 139250.82
#> 3 2.01 221243.45 93950.41 271935.00 52129.93 225979.34
#> 4 2.02 8315494.48 6518031.80 8671833.02 5135730.18 8679495.21
#> 5 6.97 26432.13 20929.65 35015.28 16707.43 33324.33
#> 6 11.96 214066.24 112212.13 270152.44 79751.55 240539.00
#> 1405_Fish2_F3 1405_Fish3_F1 1405_Fish3_F2 1405_Fish3_F3
#> 1 35190.01 73796.58 101046.21 105464.49
#> 2 63606.69 127173.02 148943.11 149585.77
#> 3 89406.65 141388.51 155992.98 166769.47
#> 4 6196265.89 8057717.46 8798917.44 9524655.11
#> 5 19249.17 30033.02 33716.19 34905.22
#> 6 104862.63 211044.91 265591.16 314792.74You can also output the quantitative results of all PMDs existing in current KEGG database.
# output all existing PMD in KEGG
keggpmd <- unique(round(keggrall$pmd,2))
data("spmeinvivo")
# remove redundant peaks
list <- globalstd(spmeinvivo)
#> 75 retention time cluster found.
#> 369 paired masses found
#> 5 unique within RT clusters high frequency PMD(s) used for further investigation.
#> The unique within RT clusters high frequency PMD(s) is(are) 28.03 21.98 44.03 17.03 18.01.
#> 719 isotopologue(s) related paired mass found.
#> 492 multi-charger(s) related paired mass found.
#> 8 retention group(s) have single peaks. 14 23 32 33 54 55 56 75
#> 11 group(s) with multiple peaks while no isotope/paired relationship 4 5 7 8 11 41 42 49 68 72 73
#> 9 group(s) with multiple peaks with isotope without paired relationship 2 9 22 26 52 62 64 66 70
#> 4 group(s) with paired relationship without isotope 1 10 15 18
#> 43 group(s) with paired relationship and isotope 3 6 12 13 16 17 19 20 21 24 25 27 28 29 30 31 34 35 36 37 38 39 40 43 44 45 46 47 48 50 51 53 57 58 59 60 61 63 65 67 69 71 74
#> 291 std mass found.
newlist <- enviGCMS::getfilter(list,rowindex = list$stdmassindex)
# generate quantitative results
pmd <- getreact(newlist,pmd=keggpmd)
# output the kegg pmd in the data
table(pmd$pmd$diff2)
#>
#> 0 0.04 1.98 2.02 2.04 2.98 3.05 3.98 4.03 6 6.05
#> 1 3 3 6 1 1 2 1 1 1 2
#> 7.03 7.06 7.2 9.95 9.98 10.02 11.96 12 12.04 12.07 12.95
#> 1 1 1 1 2 1 1 3 2 1 1
#> 12.99 13.03 13.98 14.02 14.05 15 15.01 15.02 15.05 15.79 15.94
#> 2 1 4 3 6 1 1 1 1 1 1
#> 15.99 16 16.01 16.03 16.09 16.93 17 17.03 17.99 18.01 18.06
#> 3 2 1 3 1 1 3 1 2 1 1
#> 19 19.04 19.87 20.06 21 21.98 22.05 22.76 22.8 24 25.02
#> 1 2 1 1 1 1 1 1 2 3 1
#> 25.98 26.02 26.06 28 28.03 28.07 28.1 28.99 29.03 29.96 31
#> 3 5 1 2 5 2 1 1 1 1 3
#> 31.01 31.04 31.95 31.99 33.93 34.01 34.02 34.31 35.02 35.03 38.06
#> 1 2 1 2 1 2 2 1 1 1 1
#> 39.99 40.03 40.07 40.96 42.01 42.05 42.08 43.01 43.17 45.03 45.06
#> 1 3 1 1 1 2 2 1 1 1 2
#> 45.96 46.04 47.94 48 50.02 51.05 52.01 52.07 53.05 55.06 56.03
#> 1 1 2 2 1 1 1 2 3 1 1
#> 56.04 56.06 56.27 56.94 57.02 57.06 57.97 58.04 58.08 58.75 59.04
#> 1 1 1 1 1 2 1 6 1 1 2
#> 60.08 61.02 61.22 63.03 63.96 64.02 64.03 66.05 67.03 67.07 68.03
#> 1 3 1 1 2 1 1 2 4 1 2
#> 68.06 70.01 70.04 71.04 71.98 74 75.03 75.05 75.34 77.03 77.05
#> 1 1 1 1 1 2 1 1 2 1 1
#> 78.09 78.96 79.04 80.03 80.06 80.83 81.95 81.97 82.04 82.09 83.07
#> 2 2 1 3 1 1 1 1 1 1 1
#> 84.02 84.06 86.11 88.02 88.05 90.05 90.14 90.55 93.02 93.04 94.02
#> 1 1 2 1 3 3 1 1 1 1 2
#> 94.12 94.31 95.04 95.97 98 98.07 98.08 98.14 100.08 100.95 101.05
#> 1 2 1 1 2 2 1 2 1 1 1
#> 101.96 102.01 102.14 102.25 102.97 103.03 103.27 103.83 103.91 103.99 104.06
#> 1 1 1 1 1 1 1 1 1 1 2
#> 106.03 108.02 108.06 108.09 108.13 109.03 109.05 110.04 110.83 112.12 113.07
#> 1 1 1 1 2 1 1 1 1 1 1
#> 115.83 116.95 117.06 118.03 119.04 119.06 120.06 121.87 122.01 122.02 123.01
#> 1 1 3 1 1 1 2 1 1 1 1
#> 123.99 125.02 125.88 126.03 126.1 126.18 126.96 127.12 128.01 128.12 129.08
#> 1 1 2 1 2 2 1 3 1 1 1
#> 130.85 131.09 131.86 132.94 135 135.96 137.06 139.96 142.04 142.75 142.97
#> 1 1 1 1 1 2 1 1 1 1 1
#> 143.1 143.94 144.04 144.06 145.04 146.06 147 149.94 150.07 152.12 154
#> 2 1 1 1 1 1 1 1 1 1 1
#> 157.01 157.05 157.11 158.02 158.13 158.17 158.95 159.09 159.14 160.16 161.15
#> 1 1 1 1 1 1 2 1 1 1 1
#> 161.95 162.03 162.05 163.09 165.06 168.02 168.07 170 172.05 177.06 178.18
#> 1 1 1 1 2 1 1 1 1 1 1
#> 182.17 183.13 186.03 191.99 193.01 194.04 194.24 195.02 195.03 197.05 198.01
#> 1 2 1 1 1 1 1 1 1 1 1
#> 199.96 200.18 201.04 204.19 206.06 206.08 207.02 210.2 215.05 215.97 217.07
#> 1 1 1 2 1 1 2 3 1 1 1
#> 218.17 222.2 229.04 230.09 236.93 241.05 243.02 249.09 250.09 252.22 252.92
#> 1 1 1 1 1 1 1 1 1 1 1
#> 265.03 268.22 269.1 275.12 278.02 278.16 289.17 291.1 297 300.08 300.9
#> 1 1 1 1 2 1 1 1 1 1 1
#> 301.91 302.95 304.08 306.17 311.87 315.18 316.9 339.1 343.71 344.95 345.86
#> 1 1 1 1 1 2 1 1 1 1 1
#> 346.74 348.01 349.89 369.8 375.91 381.02 390.28 405.02 409.03 409.42 409.98
#> 1 1 1 1 1 1 1 1 1 1 1
#> 431.06 437.11 439.02 440.89 447.07 454.17 454.3 455.05 456.23 463.06 466.1
#> 1 1 1 1 1 1 1 1 1 1 1
#> 472.97 473.05 475.05 476.07 478.04 488.07 501.07 501.08 505.13 516.12 517.06
#> 1 1 1 1 1 1 1 1 1 1 1
#> 518.06 526.05 537.08 546.13 546.98 552.06 554.09 586.06 599.18 612.56
#> 1 1 1 1 1 1 1 1 1 1
# output quantitative result for each PMD
head(pmd$pmddata)
#> 1405_Fish1_F1 1405_Fish1_F2 1405_Fish1_F3 1405_Fish2_F1 1405_Fish2_F2
#> 159/155 35181.202 38958.595 32009.376 51346.078 28559.260
#> 636.2/819 45321.672 33160.944 70592.139 26804.485 70922.188
#> 173.1/611 21060.045 12412.283 21212.641 8046.332 19278.585
#> 577.1/820 2860.706 1645.137 3285.761 1937.062 3347.514
#> 593.2/819 6353.084 4155.004 8805.660 3774.374 8029.871
#> 159/155 25512.422 38824.191 33660.574 48794.630 20651.311
#> 1405_Fish2_F3 1405_Fish3_F1 1405_Fish3_F2 1405_Fish3_F3
#> 159/155 41888.784 29068.087 25006.374 33043.005
#> 636.2/819 33615.276 29465.045 65983.634 82993.419
#> 173.1/611 8715.165 19432.080 24901.629 26559.379
#> 577.1/820 2016.036 2521.557 3219.349 4223.686
#> 593.2/819 4331.519 5995.990 7702.984 11416.088
#> 159/155 40127.066 35015.350 29515.079 30543.740
# output quantitative result for unique PMD
upmd <- aggregate(pmd$pmddata, by=list(pmd$pmd$diff2),sum)
# column for samples and row for unique PMD
head(upmd)
#> Group.1 1405_Fish1_F1 1405_Fish1_F2 1405_Fish1_F3 1405_Fish2_F1 1405_Fish2_F2
#> 1 0.00 72539.754 46921.790 83489.076 28793.020 80299.395
#> 2 0.04 172890.767 102824.406 209381.165 67517.693 200256.694
#> 3 1.98 112810.762 66245.180 146623.690 47083.719 139250.821
#> 4 2.02 8315494.479 6518031.801 8671833.018 5135730.184 8679495.213
#> 5 2.04 5588.112 2835.595 6882.802 1988.072 4584.455
#> 6 2.98 70741.800 36761.679 71959.565 29974.327 78883.528
#> 1405_Fish2_F3 1405_Fish3_F1 1405_Fish3_F2 1405_Fish3_F3
#> 1 35190.011 73796.583 101046.212 105464.491
#> 2 87815.305 184631.030 232005.956 239070.075
#> 3 63606.694 127173.025 148943.105 149585.769
#> 4 6196265.889 8057717.463 8798917.442 9524655.115
#> 5 2931.107 4411.049 5688.054 5438.388
#> 6 40546.116 62671.685 79751.686 82474.088Reactomics analysis for MS only data
When retention time is not provided, m/z vector can still be used to
check reaction level changes. You might use getrda to find
the high frequency PMDs.
data(spmeinvivo)
# get the m/z
mz <- spmeinvivo$mz
# get the m/z intensity for all m/z, the row order is the same with mz
insms <- spmeinvivo$data
# check high frequency pmd
sda <- getrda(mz)
#> 164462 pmd found.
#> 20 pmd used.
colnames(sda)
#> [1] "0" "1.001" "1.002" "1.003" "1.004" "2.015" "2.016"
#> [8] "14.015" "17.026" "18.011" "21.982" "28.031" "28.032" "44.026"
#> [15] "67.987" "67.988" "88.052" "116.192" "135.974" "135.975"
# save them as numeric vector
hfpmd <- as.numeric(colnames(sda))Then getpmddf function can be used to extract all the
paired ions for certain PMD.
# get details for certain pmd
pmddf <- getpmddf(mz,pmd=18.011,digits = 3)
# add intensity for all the paired ions
mz1ins <- insms[match(pmddf$ms1,mz),]
mz2ins <- insms[match(pmddf$ms2,mz),]
# get the pmd pair intensity
pmdins <- mz1ins+mz2ins
# get the pmd total intensity across samples
pmdinsall <- apply(pmdins,2,sum)
# show the PMD intensity
pmdinsall
#> 1405_Fish1_F1 1405_Fish1_F2 1405_Fish1_F3 1405_Fish2_F1 1405_Fish2_F2
#> 9898514 7801273 10363201 5847334 10479551
#> 1405_Fish2_F3 1405_Fish3_F1 1405_Fish3_F2 1405_Fish3_F3
#> 7021375 10584976 12989961 12559649You can also calculate the static or dynamic PMD intensity for m/z only data.
# get the ratio of larger m/z over smaller m/z
ratio <- mz2ins/mz1ins
# filter PMD based on RSD% across samples
# cutoff 30%
cutoff <- 0.3
# get index for static PMD
rsdidx <- apply(ratio,1,function(x) sd(x)/mean(x)<cutoff)
# get static PMD
pmddfstatic <- pmddf[rsdidx,]
# get static intensity
pmdinsstatic <- pmdins[rsdidx,]
# normalize the ions pair intensity to avoid influences from large response factors
pmdinsstaticscale <- t(scale(t(pmdinsstatic)))
# get the pmd static intensity across samples
pmdinsstaticall <- apply(pmdinsstaticscale,2,sum)
# show the PMD static intensity for each sample
pmdinsstaticall
#> 1405_Fish1_F1 1405_Fish1_F2 1405_Fish1_F3 1405_Fish2_F1 1405_Fish2_F2
#> 1.027100 -16.703628 2.373645 -27.241358 12.433566
#> 1405_Fish2_F3 1405_Fish3_F1 1405_Fish3_F2 1405_Fish3_F3
#> -17.757951 7.924290 19.802820 18.141514
# get index for dynamic PMD
rsdidx <- apply(ratio,1,function(x) sd(x)/mean(x)>=cutoff)
# get dynamic PMD
pmddfdynamic <- pmddf[rsdidx,]
# get dynamic intensity for ms1 and ms2
pmdinsdynamicms1 <- apply(mz1ins[rsdidx,],1,function(x) sd(x)/mean(x))
pmdinsdynamicms2 <- apply(mz2ins[rsdidx,],1,function(x) sd(x)/mean(x))
# find the stable ms and use ratio as intensity
idx <- pmdinsdynamicms1>pmdinsdynamicms2
pmdinsdynamic <- ratio[rsdidx,]
pmdinsdynamic[idx,] <- 1/ratio[rsdidx,][idx,]
# get the pmd dynamic intensity across samples
pmdinsdynamicall <- apply(pmdinsdynamic,2,sum)
# show the PMD dynamic intensity for each sample
pmdinsdynamicall
#> 1405_Fish1_F1 1405_Fish1_F2 1405_Fish1_F3 1405_Fish2_F1 1405_Fish2_F2
#> 374.1957 315.6090 388.0063 207.8236 233.4380
#> 1405_Fish2_F3 1405_Fish3_F1 1405_Fish3_F2 1405_Fish3_F3
#> 199.9064 283.4682 327.9651 256.2044You can also use getpmddf function extract all the
paired ions for multiple PMDs. Then you could generate the network based
on the output.
# get details for certain pmd
pmddf <- getpmddf(mz,pmd=hfpmd,digits = 3)
# viz by igraph package
library(igraph)
net <- graph_from_data_frame(pmddf,directed = F)
pal <- grDevices::rainbow(length(unique(E(net)$diff2)))
plot(net,vertex.label=NA,vertex.size = 5,edge.width = 3,edge.color = pal[as.numeric(as.factor(E(net)$diff2))],main = 'PMD network')
legend("topright",bty = "n",
legend=unique(E(net)$diff2),
fill=unique(pal[as.numeric(as.factor(E(net)$diff2))]), border=NA,horiz = F)
If you prefer to get a pmd network for a specific mass. You can still
use getchain function.
data(spmeinvivo)
spmeinvivo$rt <- NULL
chain <- getchain(spmeinvivo,diff = c(2.02,14.02,15.99,58.04,13.98),mass = 286.3101,digits = 2,corcutoff = 0)
# show as network
net <- graph_from_data_frame(chain$sdac,directed = F)
pal <- grDevices::rainbow(5)
plot(net,vertex.label=round(as.numeric(V(net)$name),2),vertex.size =5,edge.width = 3,edge.color = pal[as.numeric(as.factor(E(net)$diff2))],vertex.label.dist=1,vertex.color=ifelse(round(as.numeric(V(net)$name),4) %in% 286.3101,'red','black'), main = 'PMD network')
legend("topright",bty = "n",
legend=unique(E(net)$diff2),
fill=unique(pal[as.numeric(as.factor(E(net)$diff2))]), border=NA,horiz = F)
PMD Reaction Database
To check the pmd reaction database:
# all reaction
data("omics")
head(omics)
#> X KEGG RHEA_ID DIRECTION MASTER_ID ec ecocyc macie
#> 1 1 R00001 22455 BI 22452 <NA> <NA> <NA>
#> 2 2 R00004 24579 BI 24576 <NA> <NA> <NA>
#> 3 3 R00005 19032 BI 19029 <NA> <NA> <NA>
#> 4 4 R00008 22751 BI 22748 <NA> <NA> <NA>
#> 5 5 R00009 20312 BI 20309 <NA> <NA> <NA>
#> 6 6 R00010 20871 BI 20868 <NA> <NA> <NA>
#> metacyc reactome
#> 1 ENDOPOLYPHOSPHATASE-RXN <NA>
#> 2 INORGPYROPHOSPHAT-RXN <NA>
#> 3 <NA> <NA>
#> 4 <NA> <NA>
#> 5 <NA> <NA>
#> 6 <NA> <NA>
#> compounds pmd
#> 1 [phosphate](n+1) + n H2O <=> n H(+) + (n+1) phosphate NA
#> 2 diphosphate + H2O <=> H(+) + 2 phosphate NA
#> 3 3 H(+) + H2O + urea-1-carboxylate <=> 2 CO2 + 2 NH4(+) NA
#> 4 4-hydroxy-4-methyl-2-oxoglutarate <=> 2 pyruvate NA
#> 5 2 H2O2 <=> 2 H2O + O2 NA
#> 6 alpha,alpha-trehalose + H2O <=> 2 D-glucose 162.053
# kegg reaction
data("keggrall")
head(keggrall)
#> rid rc ez ms1 formula1
#> 1 R00002 RC00002 1.18.6.1 506.9957 C10H16N5O13P3
#> 2 R00010 RC00049 3.2.1.28 342.1162 C12H22O11
#> 3 R00012 RC00002 2.7.7.45 522.9907 C10H16N5O14P3
#> 4 R00014 RC00027;RC02744 1.2.4.1;2.2.1.6;4.1.1.1 425.0450 C12H19N4O7P2S
#> 5 R00015 RC00077;RC00247 2.4.1.99 342.1162 C12H22O11
#> 6 R00015 RC00077;RC00247 2.4.1.99 342.1162 C12H22O11
#> ms2 formula2 pmd C H O N P S
#> 1 427.0294 C10H15N5O10P2 79.966 0 1 3 0 1 0
#> 2 180.0634 C6H12O6 162.053 6 10 5 0 0 0
#> 3 868.0381 C20H28N10O21P4 345.047 10 12 7 5 1 0
#> 4 469.0712 C14H23N4O8P2S 44.026 2 4 1 0 0 0
#> 5 180.0634 C6H12O6 162.053 6 10 5 0 0 0
#> 6 504.1690 C18H32O16 162.053 6 10 5 0 0 0
# literature reaction for mass spectrometry
data("sda")
head(sda)
#> PMD origin
#> 1 0.984016 OH ↔ NH2, e.g. de-amidiation, CHNO compounds
#> 2 1.995663 F ↔ OH, halogen exchange with hydroxy group (typically -F + OH)
#> 3 2.015650 ± 2H, opening or forming of double bond
#> 4 7.004671 F ↔ CN, halogen exchange with cyano group
#> 5 8.965779 Cl ↔ CN, halogen exchange with cyano group
#> 6 13.979265 O ↔ 2H, e.g. Oxidation follwed by H2O elimination
#> Ref. mode
#> 1 https://doi.org/10.1016/S1044-0305(99)00090-2 both
#> 2 https://doi.org/10.1016/S1044-0305(99)00090-2 both
#> 3 https://doi.org/10.1016/S1044-0305(99)00090-2 both
#> 4 https://doi.org/10.1016/S1044-0305(99)00090-2 both
#> 5 https://doi.org/10.1016/S1044-0305(99)00090-2 both
#> 6 https://doi.org/10.1016/S1044-0305(99)00090-2 bothTo check the HMDB pmd database:
data("hmdb")
head(hmdb)
#> C H O N P S pmd pmd2 percentage
#> 1 -3 4 2 0 0 0 0.021 0.02 0.9623060
#> 2 1 4 -1 0 0 0 0.036 0.04 0.9503645
#> 3 5 4 -4 0 0 0 0.052 0.05 0.9412861
#> 4 2 8 -2 0 0 0 0.073 0.07 0.9617834
#> 5 6 8 -5 0 0 0 0.088 0.09 0.9247730
#> 6 -1 12 0 0 0 0 0.094 0.09 0.9648936To extract any compound KEGG compound’s pmd network with known PMD:
plotcn('C6H12O6','Glucose',c(2.016,14.016,15.995))
#> Average distance of PMD network is 9.42710759198882
#> Average degree 2.8421052631579