diff --git a/assays/canal-1_Phenotyping_2019/dataset/canal-1_yield_data.xlsx b/assays/canal-1_Phenotyping_2019/dataset/canal-1_yield_data.xlsx
new file mode 100644
index 0000000000000000000000000000000000000000..04fab0fd8d9f30827bd95147a8a1eede84431b98
--- /dev/null
+++ b/assays/canal-1_Phenotyping_2019/dataset/canal-1_yield_data.xlsx
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3bc80d6593d3bd97a204ded64ecb2ec569149d38ed14592fe529491647f74221
+size 9944
diff --git a/runs/Met_whole/Bubble_comb.png b/runs/Met_whole/Bubble_comb.png
new file mode 100644
index 0000000000000000000000000000000000000000..6c9d248db18ba12059fc519c99e1624341bb4498
Binary files /dev/null and b/runs/Met_whole/Bubble_comb.png differ
diff --git a/runs/Met_whole/Combined_bubble_metabolites.RDS b/runs/Met_whole/Combined_bubble_metabolites.RDS
new file mode 100644
index 0000000000000000000000000000000000000000..01ead4804f45085c2b90f73716b28adabf5009a0
Binary files /dev/null and b/runs/Met_whole/Combined_bubble_metabolites.RDS differ
diff --git a/runs/Met_whole/Heatmap_all.jpg b/runs/Met_whole/Heatmap_all.jpg
index cb786e9687adcdec7f20090ccc4c097f15aff446..fdba6da94b86efd6746252d9f57882d37a8dbca3 100644
Binary files a/runs/Met_whole/Heatmap_all.jpg and b/runs/Met_whole/Heatmap_all.jpg differ
diff --git a/runs/Met_whole/Heatmap_all_met.RDS b/runs/Met_whole/Heatmap_all_met.RDS
index 20d9f3df2d288a0a18e9bd98bbb5f945c2075d14..376214e087b91510c7b2885cd03a160ff1b84b78 100644
Binary files a/runs/Met_whole/Heatmap_all_met.RDS and b/runs/Met_whole/Heatmap_all_met.RDS differ
diff --git a/runs/Met_whole/Metaboanalyst_Bubble_WTglB_vs_MUgl.png b/runs/Met_whole/Metaboanalyst_Bubble_WTglB_vs_MUgl.png
new file mode 100644
index 0000000000000000000000000000000000000000..47c5d761b0d3aa4c5af663685b3877a88f26dfb9
Binary files /dev/null and b/runs/Met_whole/Metaboanalyst_Bubble_WTglB_vs_MUgl.png differ
diff --git a/runs/Met_whole/Metaboanalyst_Bubble_WTglB_vs_MUwl.png b/runs/Met_whole/Metaboanalyst_Bubble_WTglB_vs_MUwl.png
new file mode 100644
index 0000000000000000000000000000000000000000..5567d490c5fb8ccdd5a066c95aadefdb09f6a321
Binary files /dev/null and b/runs/Met_whole/Metaboanalyst_Bubble_WTglB_vs_MUwl.png differ
diff --git a/runs/Met_whole/Metaboanalyst_bubble_WTglB_vs_MUgl.xlsx b/runs/Met_whole/Metaboanalyst_bubble_WTglB_vs_MUgl.xlsx
new file mode 100644
index 0000000000000000000000000000000000000000..25856bfa40e525f1d855907af473d856b5446b35
Binary files /dev/null and b/runs/Met_whole/Metaboanalyst_bubble_WTglB_vs_MUgl.xlsx differ
diff --git a/runs/Met_whole/Metaboanalyst_bubble_WTglB_vs_MUwl.xlsx b/runs/Met_whole/Metaboanalyst_bubble_WTglB_vs_MUwl.xlsx
new file mode 100644
index 0000000000000000000000000000000000000000..af9f1ac2114d47a4b76ba33cdec1e5e94c2c4879
Binary files /dev/null and b/runs/Met_whole/Metaboanalyst_bubble_WTglB_vs_MUwl.xlsx differ
diff --git a/runs/Met_whole/asparagine.png b/runs/Met_whole/asparagine.png
index 76c43aba1c1a811f22c7b74b41fa38e042145320..91f5ca98c1bac83b9da07f6111567e324de1d836 100644
Binary files a/runs/Met_whole/asparagine.png and b/runs/Met_whole/asparagine.png differ
diff --git a/runs/Met_whole/asparagine_leaves.RDS b/runs/Met_whole/asparagine_leaves.RDS
index 28beffc35f1fb4d2f7c087ad1bcadd1621e9a8a5..b7e16b7e058778a44aa6835448f5e1e3a2147af0 100644
Binary files a/runs/Met_whole/asparagine_leaves.RDS and b/runs/Met_whole/asparagine_leaves.RDS differ
diff --git a/runs/Met_whole/cana1_gl_pathway_results.csv b/runs/Met_whole/cana1_gl_pathway_results.csv
new file mode 100644
index 0000000000000000000000000000000000000000..cba40447cd40e9bdce5060ef4e7bfda893db1363
--- /dev/null
+++ b/runs/Met_whole/cana1_gl_pathway_results.csv
@@ -0,0 +1,59 @@
+,Total Cmpd,Hits,Raw p,#NAME?,Holm adjust,FDR,Impact
+Fructose and mannose metabolism,20,1,0.00029287,3.5333,0.016986,0.016986,0
+"Alanine, aspartate and glutamate metabolism",22,10,0.0010156,2.9933,0.05789,0.023646,0.84892
+Glutathione metabolism,26,7,0.0014774,2.8305,0.082733,0.023646,0.13362
+Cyanoamino acid metabolism,29,6,0.0016307,2.7876,0.089691,0.023646,0.35593
+Selenocompound metabolism,13,1,0.004896,2.3102,0.26438,0.056793,0
+Nitrogen metabolism,12,2,0.006858,2.1638,0.36348,0.057737,0
+Glyoxylate and dicarboxylate metabolism,29,11,0.007169,2.1445,0.37279,0.057737,0.50324
+Aminoacyl-tRNA biosynthesis,46,19,0.0079637,2.0989,0.40615,0.057737,0.11111
+Arginine biosynthesis,18,8,0.0093094,2.0311,0.46547,0.059994,0.34079
+Pyrimidine metabolism,38,3,0.010468,1.9801,0.51293,0.060714,0.05766
+Purine metabolism,63,4,0.011995,1.921,0.57578,0.063248,0.07925
+Arginine and proline metabolism,34,8,0.0155,1.8097,0.72852,0.074919,0.5247
+"Glycine, serine and threonine metabolism",33,8,0.018263,1.7384,0.84011,0.078969,0.60216
+Cysteine and methionine metabolism,46,6,0.019406,1.7121,0.87326,0.078969,0.16135
+Sulfur metabolism,15,3,0.020423,1.6899,0.89862,0.078969,0.10774
+Sphingolipid metabolism,17,1,0.026431,1.5779,1,0.090635,0
+"Valine, leucine and isoleucine biosynthesis",22,6,0.026565,1.5757,1,0.090635,0.10727
+Galactose metabolism,27,7,0.029971,1.5233,1,0.095385,0.29394
+Porphyrin and chlorophyll metabolism,48,1,0.031247,1.5052,1,0.095385,0
+Phosphatidylinositol signaling system,26,2,0.039039,1.4085,1,0.11321,0.03285
+Inositol phosphate metabolism,28,3,0.044463,1.352,1,0.1228,0.12552
+Carbon fixation in photosynthetic organisms,21,5,0.050027,1.3008,1,0.13189,0.10075
+Butanoate metabolism,17,5,0.05963,1.2245,1,0.15037,0.13636
+"Valine, leucine and isoleucine degradation",37,3,0.071051,1.1484,1,0.17171,0
+Glycerophospholipid metabolism,37,2,0.08671,1.0619,1,0.19859,0.09832
+Glycerolipid metabolism,21,4,0.089437,1.0485,1,0.19859,0.02235
+"Phenylalanine, tyrosine and tryptophan biosynthesis",22,4,0.092521,1.0338,1,0.19859,0.1016
+Pantothenate and CoA biosynthesis,23,4,0.095883,1.0183,1,0.19859,0.08423
+Nicotinate and nicotinamide metabolism,13,2,0.099296,1.0031,1,0.19859,0
+Glucosinolate biosynthesis,65,6,0.10452,0.98082,1,0.20206,0
+beta-Alanine metabolism,18,4,0.11245,0.94905,1,0.21039,0.32937
+Monobactam biosynthesis,8,2,0.12424,0.90576,1,0.22518,0
+Ascorbate and aldarate metabolism,18,5,0.1425,0.84619,1,0.2476,0.40299
+Citrate cycle (TCA cycle),20,7,0.14515,0.83819,1,0.2476,0.32413
+Pyruvate metabolism,22,4,0.16038,0.79485,1,0.24955,0.32193
+Phenylalanine metabolism,11,1,0.16098,0.79322,1,0.24955,0.47059
+"Tropane, piperidine and pyridine alkaloid biosynthesis",8,1,0.16098,0.79322,1,0.24955,0
+Lysine biosynthesis,9,3,0.16758,0.77578,1,0.24955,0
+Isoquinoline alkaloid biosynthesis,6,2,0.1678,0.77521,1,0.24955,0.5
+Propanoate metabolism,20,1,0.17295,0.76207,1,0.25078,0
+Starch and sucrose metabolism,22,7,0.33099,0.48019,1,0.46823,0.72431
+Tyrosine metabolism,16,4,0.34913,0.45701,1,0.48214,0.27703
+Lysine degradation,18,3,0.36048,0.44312,1,0.48622,0
+Biosynthesis of secondary metabolites - unclassified,5,1,0.40582,0.39166,1,0.52567,1
+Ubiquinone and other terpenoid-quinone biosynthesis,38,2,0.40785,0.3895,1,0.52567,0.00097
+Histidine metabolism,15,1,0.41804,0.37878,1,0.52709,0.04264
+Phenylpropanoid biosynthesis,46,3,0.46851,0.32928,1,0.57816,0.05295
+Amino sugar and nucleotide sugar metabolism,50,5,0.51068,0.29185,1,0.60855,0.10791
+Fatty acid biosynthesis,56,2,0.51412,0.28894,1,0.60855,0
+Glycolysis / Gluconeogenesis,26,3,0.6327,0.1988,1,0.73393,0.12036
+Pentose and glucuronate interconversions,16,3,0.657,0.18243,1,0.74718,0
+Pentose phosphate pathway,19,1,0.72613,0.13899,1,0.80062,0.11621
+Terpenoid backbone biosynthesis,30,1,0.7316,0.13573,1,0.80062,0
+Thiamine metabolism,22,2,0.7646,0.11657,1,0.82123,0
+Tryptophan metabolism,28,1,0.81944,0.086485,1,0.8487,0.12037
+Indole alkaloid biosynthesis,4,1,0.81944,0.086485,1,0.8487,0
+C5-Branched dibasic acid metabolism,6,2,0.93199,0.030589,1,0.94834,0
+Zeatin biosynthesis,21,2,0.97936,0.0090588,1,0.97936,0
diff --git a/runs/Met_whole/cana1_wl_pathway_results.csv b/runs/Met_whole/cana1_wl_pathway_results.csv
new file mode 100644
index 0000000000000000000000000000000000000000..84111004e8bfcae37b3e4f3c4927437b44b05c84
--- /dev/null
+++ b/runs/Met_whole/cana1_wl_pathway_results.csv
@@ -0,0 +1,59 @@
+,Total Cmpd,Hits,Raw p,#NAME?,Holm adjust,FDR,Impact
+Fructose and mannose metabolism,20,1,1.30E-08,7.8867,7.53E-07,7.53E-07,0
+Pantothenate and CoA biosynthesis,23,4,7.69E-08,7.1142,4.38E-06,1.69E-06,0.08423
+"Valine, leucine and isoleucine biosynthesis",22,6,1.18E-07,6.9289,6.60E-06,1.69E-06,0.10727
+Terpenoid backbone biosynthesis,30,1,1.63E-07,6.7886,8.95E-06,1.69E-06,0
+"Valine, leucine and isoleucine degradation",37,3,1.66E-07,6.7811,8.95E-06,1.69E-06,0
+Glucosinolate biosynthesis,65,6,1.75E-07,6.7567,9.28E-06,1.69E-06,0
+Tryptophan metabolism,28,1,2.89E-07,6.539,1.50E-05,2.10E-06,0.12037
+Indole alkaloid biosynthesis,4,1,2.89E-07,6.539,1.50E-05,2.10E-06,0
+Glutathione metabolism,26,7,5.86E-07,6.2319,2.93E-05,3.78E-06,0.13362
+Phosphatidylinositol signaling system,26,2,1.15E-06,5.9387,5.64E-05,6.68E-06,0.03285
+Inositol phosphate metabolism,28,3,2.03E-06,5.6935,9.72E-05,1.07E-05,0.12552
+Sulfur metabolism,15,3,1.08E-05,4.9681,0.00050581,5.20E-05,0.10774
+C5-Branched dibasic acid metabolism,6,2,2.11E-05,4.6765,0.00096885,9.40E-05,0
+Arginine and proline metabolism,34,8,3.60E-05,4.4434,0.0016212,0.00014926,0.5247
+Pyrimidine metabolism,38,3,4.58E-05,4.3393,0.0020144,0.00016923,0.05766
+Glyoxylate and dicarboxylate metabolism,29,11,4.67E-05,4.3308,0.0020144,0.00016923,0.50324
+Nitrogen metabolism,12,2,5.26E-05,4.2794,0.0022072,0.00017929,0
+Purine metabolism,63,4,5.74E-05,4.2409,0.0023544,0.00018504,0.07925
+Tyrosine metabolism,16,4,9.40E-05,4.027,0.0037586,0.00026937,0.27703
+Isoquinoline alkaloid biosynthesis,6,2,9.56E-05,4.0196,0.0037586,0.00026937,0.5
+Ubiquinone and other terpenoid-quinone biosynthesis,38,2,9.75E-05,4.0109,0.0037586,0.00026937,0.00097
+"Phenylalanine, tyrosine and tryptophan biosynthesis",22,4,0.0001319,3.8798,0.0048803,0.00034773,0.1016
+Starch and sucrose metabolism,22,7,0.00033327,3.4772,0.011998,0.00084042,0.72431
+Glycolysis / Gluconeogenesis,26,3,0.00047536,3.323,0.016638,0.0011052,0.12036
+Citrate cycle (TCA cycle),20,7,0.00047637,3.3221,0.016638,0.0011052,0.32413
+Phenylpropanoid biosynthesis,46,3,0.00088182,3.0546,0.0291,0.0018302,0.05295
+Phenylalanine metabolism,11,1,0.00088356,3.0538,0.0291,0.0018302,0.47059
+"Tropane, piperidine and pyridine alkaloid biosynthesis",8,1,0.00088356,3.0538,0.0291,0.0018302,0
+Glycerolipid metabolism,21,4,0.0018387,2.7355,0.05516,0.0036773,0.02235
+Sphingolipid metabolism,17,1,0.0019075,2.7195,0.055316,0.0036878,0
+Propanoate metabolism,20,1,0.0029375,2.532,0.082249,0.0054959,0
+Arginine biosynthesis,18,8,0.003661,2.4364,0.098847,0.0066355,0.34079
+Aminoacyl-tRNA biosynthesis,46,19,0.0047869,2.3199,0.12446,0.0084134,0.11111
+Galactose metabolism,27,7,0.0091016,2.0409,0.22754,0.015526,0.29394
+Butanoate metabolism,17,5,0.0098481,2.0066,0.23635,0.01632,0.13636
+Porphyrin and chlorophyll metabolism,48,1,0.011279,1.9477,0.25942,0.018172,0
+"Alanine, aspartate and glutamate metabolism",22,10,0.017776,1.7502,0.39107,0.027865,0.84892
+Cyanoamino acid metabolism,29,6,0.023478,1.6293,0.49304,0.035835,0.35593
+Lysine degradation,18,3,0.026983,1.5689,0.53967,0.039744,0
+Histidine metabolism,15,1,0.02741,1.5621,0.53967,0.039744,0.04264
+Pentose phosphate pathway,19,1,0.032153,1.4928,0.57875,0.045484,0.11621
+Thiamine metabolism,22,2,0.047944,1.3193,0.81504,0.066208,0
+Selenocompound metabolism,13,1,0.05842,1.2334,0.93472,0.078799,0
+Glycerophospholipid metabolism,37,2,0.06874,1.1628,1,0.090203,0.09832
+Lysine biosynthesis,9,3,0.069985,1.155,1,0.090203,0
+Pyruvate metabolism,22,4,0.078343,1.106,1,0.09878,0.32193
+Amino sugar and nucleotide sugar metabolism,50,5,0.080745,1.0929,1,0.099642,0.10791
+Ascorbate and aldarate metabolism,18,5,0.088924,1.051,1,0.10745,0.40299
+Pentose and glucuronate interconversions,16,3,0.12346,0.90847,1,0.14614,0
+Carbon fixation in photosynthetic organisms,21,5,0.16395,0.7853,1,0.19018,0.10075
+Cysteine and methionine metabolism,46,6,0.21104,0.67563,1,0.23597,0.16135
+beta-Alanine metabolism,18,4,0.21156,0.67456,1,0.23597,0.32937
+"Glycine, serine and threonine metabolism",33,8,0.22696,0.64406,1,0.24837,0.60216
+Nicotinate and nicotinamide metabolism,13,2,0.25348,0.59605,1,0.26781,0
+Monobactam biosynthesis,8,2,0.25396,0.59524,1,0.26781,0
+Biosynthesis of secondary metabolites - unclassified,5,1,0.44773,0.34898,1,0.46372,1
+Fatty acid biosynthesis,56,2,0.50236,0.29899,1,0.51117,0
+Zeatin biosynthesis,21,2,0.81428,0.089227,1,0.81428,0
diff --git a/runs/Met_whole/complete_norm.qs b/runs/Met_whole/complete_norm.qs
new file mode 100644
index 0000000000000000000000000000000000000000..cfbff313d9aa44244057a5960528e7f3fda9a2a1
Binary files /dev/null and b/runs/Met_whole/complete_norm.qs differ
diff --git a/runs/Met_whole/compound_db.qs b/runs/Met_whole/compound_db.qs
new file mode 100644
index 0000000000000000000000000000000000000000..21b880d7bd04e101d01a9527d18079a5550b8220
Binary files /dev/null and b/runs/Met_whole/compound_db.qs differ
diff --git a/runs/Met_whole/data_orig.qs b/runs/Met_whole/data_orig.qs
new file mode 100644
index 0000000000000000000000000000000000000000..eac70b4d67751c5f54deca0e8768817f3fb548f1
Binary files /dev/null and b/runs/Met_whole/data_orig.qs differ
diff --git a/runs/Met_whole/data_proc.qs b/runs/Met_whole/data_proc.qs
new file mode 100644
index 0000000000000000000000000000000000000000..cfbff313d9aa44244057a5960528e7f3fda9a2a1
Binary files /dev/null and b/runs/Met_whole/data_proc.qs differ
diff --git a/runs/Met_whole/name_map.csv b/runs/Met_whole/name_map.csv
new file mode 100644
index 0000000000000000000000000000000000000000..8d8ea26d1ed97308347a6fb9bcf6663b14b76585
--- /dev/null
+++ b/runs/Met_whole/name_map.csv
@@ -0,0 +1,101 @@
+"Query","Match","HMDB","PubChem","KEGG","SMILES","Comment"
+"group","NA","NA","NA","NA","NA","0"
+"C00147","Adenine","HMDB0000034","190","C00147","NC1=C2NC=NC2=NC=N1","1"
+"C00020","Adenosine monophosphate","HMDB0000045","6083","C00020","NC1=C2N=CN([C@@H]3O[C@H](COP(O)(O)=O)[C@@H](O)[C@H]3O)C2=NC=N1","1"
+"C00956","Aminoadipic acid","HMDB0000510","92136","C00956","N[C@@H](CCCC(O)=O)C(O)=O","1"
+"C06104","Adipic acid","HMDB0000448","196","C06104","OC(=O)CCCCC(O)=O","1"
+"C00041","L-Alanine","HMDB0000161","5950","C00041","C[C@H](N)C(O)=O","1"
+"C02512","3-Cyano-L-alanine","METPA0300",NA,"C02512",NA,"1"
+"C00099","beta-Alanine","HMDB0000056","239","C00099","NCCC(O)=O","1"
+"C00062","L-Arginine","HMDB0000517","6322","C00062","N[C@@H](CCCNC(N)=N)C(O)=O","1"
+"C00072","Ascorbic acid","HMDB0000044","54670067","C00072","[H][C@@]1(OC(=O)C(O)=C1O)[C@@H](O)CO","1"
+"C00152","L-Asparagine","HMDB0000168","6267","C00152","N[C@@H](CC(N)=O)C(O)=O","1"
+"C00049","L-Aspartic acid","HMDB0000191","5960","C00049","N[C@@H](CC(O)=O)C(O)=O","1"
+"C00156","4-Hydroxybenzoic acid","HMDB0000500","135","C00156","OC(=O)C1=CC=C(O)C=C1","1"
+"C00180","Benzoic acid","HMDB0001870","243","C00180","OC(=O)C1=CC=CC=C1","1"
+"C00556","Benzyl alcohol","HMDB0003119","244","C00556","OCC1=CC=CC=C1","1"
+"C00334","gamma-Aminobutyric acid","HMDB0000112","119","C00334","NCCCC(O)=O","1"
+"C10850","NA","NA","NA","NA","NA","0"
+"C10851","Calystegine B2","HMDB0038594","3693124","C10851","OC1C2CCC(O)(N2)C(O)C1O","1"
+"C00811","4-Hydroxycinnamic acid","HMDB0002035","637542","C00811","OC(=O)\C=C\C1=CC=C(O)C=C1","1"
+"C00158","Citric acid","HMDB0000094","311","C00158","OC(=O)CC(O)(CC(O)=O)C(O)=O","1"
+"C01571","Capric acid","HMDB0000511","2969","C01571","CCCCCCCCCC(O)=O","1"
+"C00503","Erythritol","HMDB0002994","222285","C00503","OC[C@H](O)[C@H](O)CO","1"
+"C21593","NA","NA","NA","NA","NA","0"
+"C00189","Ethanolamine","HMDB0000149","700","C00189","NCCO","1"
+"C00085","Fructose 6-phosphate","HMDB0000124","69507","C00085","OCC(=O)[C@@H](O)[C@H](O)[C@H](O)COP(O)(O)=O","1"
+"C00095","D-Fructose","HMDB0000660","439709","C00095","OC[C@H]1O[C@](O)(CO)[C@@H](O)[C@@H]1O","1"
+"C01019","L-Fucose","HMDB0000174","17106","C01019","C[C@@H]1O[C@@H](O)[C@@H](O)[C@H](O)[C@@H]1O","1"
+"C00122","Fumaric acid","HMDB0000134","444972","C00122","OC(=O)\C=C\C(O)=O","1"
+"C01235","beta-Cortol","HMDB0005821","439451","C01235","[H][C@@]12CC[C@](O)([C@H](O)CO)[C@@]1(C)C[C@H](O)[C@@]1([H])[C@@]2([H])CC[C@]2([H])C[C@H](O)CC[C@]12C","1"
+"C01115","NA","NA","NA","NA","NA","0"
+"C00333","Galacturonic acid","HMDB0002545","84740","C00333","O[C@@H](C=O)[C@@H](O)[C@@H](O)[C@H](O)C(O)=O","1"
+"C00345","6-Phosphogluconic acid","HMDB0001316","91493","C00345","O[C@H](COP(O)(O)=O)[C@@H](O)[C@H](O)[C@@H](O)C(O)=O","1"
+"C00103","Glucose 1-phosphate","HMDB0001586","65533","C00103","OC[C@H]1O[C@H](OP(O)(O)=O)[C@H](O)[C@@H](O)[C@@H]1O","1"
+"C00092","Glucose 6-phosphate","HMDB0001401","5958","C00092","OC1O[C@H](COP(O)(O)=O)[C@@H](O)[C@H](O)[C@H]1O","1"
+"C22350","NA","NA","NA","NA","NA","0"
+"C00031","D-Glucose","HMDB0000122","5793","C00031","OC[C@H]1O[C@@H](O)[C@H](O)[C@@H](O)[C@@H]1O","1"
+"C00025","Glutamic acid","HMDB0000148","33032","C00025","N[C@@H](CCC(O)=O)C(O)=O","1"
+"C00064","Glutamine","HMDB0000641","5961","C00064","N[C@@H](CCC(N)=O)C(O)=O","1"
+"C00026","Oxoglutaric acid","HMDB0000208","51","C00026","OC(=O)CCC(=O)C(O)=O","1"
+"C00489","Glutaric acid","HMDB0000661","743","C00489","OC(=O)CCCC(O)=O","1"
+"C00258","Glyceric acid","HMDB0000139","439194","C00258","OC[C@@H](O)C(O)=O","1"
+"C00093","Glycerol 3-phosphate","HMDB0000126","439162","C00093","OC[C@@H](O)COP(O)(O)=O","1"
+"C00116","Glycerol","HMDB0000131","753","C00116","OCC(O)CO","1"
+"C00037","Glycine","HMDB0000123","750","C00037","NCC(O)=O","1"
+"C00160","Glycolate",NA,"3460","C00160",NA,"1"
+"C17349","Guanidine","HMDB0001842","3520","C17349","NC(N)=N","1"
+"C01040","L-Gulonolactone","HMDB0003466","439373","C01040","[H][C@@]1(OC(=O)[C@@H](O)[C@H]1O)[C@@H](O)CO","1"
+"C00135","Histidine","HMDB0000177","6274","C00135","N[C@@H](CC1=CN=CN1)C(O)=O","1"
+"C00263","L-Homoserine","HMDB0000719","12647","C00263","N[C@@H](CCO)C(O)=O","1"
+"C04006","myo-Inositol 1-phosphate","HMDB0000213",NA,"C04006","O[C@H]1[C@H](O)[C@@H](O)[C@H](OP(O)(O)=O)[C@H](O)[C@@H]1O","1"
+"C00137","myo-Inositol","HMDB0000211",NA,"C00137","O[C@H]1[C@H](O)[C@@H](O)[C@H](O)[C@H](O)[C@@H]1O","1"
+"C00311","Isocitric acid","HMDB0000193","1198","C00311","OC(C(CC(O)=O)C(O)=O)C(O)=O","1"
+"C00407","Isoleucine","HMDB0000172","6306","C00407","CC[C@H](C)[C@H](N)C(O)=O","1"
+"C00490","Itaconic acid","HMDB0002092","811","C00490","OC(=O)CC(=C)C(O)=O","1"
+"C00186","Lactic acid","HMDB0000190","107689","C00186","C[C@H](O)C(O)=O","1"
+"C00123","Leucine","HMDB0000687","6106","C00123","CC(C)C[C@H](N)C(O)=O","1"
+"C00047","Lysine","HMDB0000182","5962","C00047","NCCCC[C@H](N)C(O)=O","1"
+"C00476","Isocarboxazid","HMDB0015377","3759","C00476","CC1=CC(=NO1)C(=O)NNCC1=CC=CC=C1","1"
+"C02612","(R)-2-Methylmalate","METPA0308",NA,"C02612",NA,"1"
+"C00149","Malic acid","HMDB0000156","222656","C00149","O[C@@H](CC(O)=O)C(O)=O","1"
+"C00185","Cellobiose","HMDB0000055","10712","C00185","OC[C@H]1O[C@@H](O[C@H]2[C@H](O)[C@@H](O)[C@H](O)O[C@@H]2CO)[C@H](O)[C@@H](O)[C@@H]1O","1"
+"C00208","D-Maltose","HMDB0000163","10991489","C00208","OC[C@H]1O[C@H](O[C@H]2[C@H](O)[C@H](O)[C@@H](O)O[C@@H]2CO)[C@H](O)[C@@H](O)[C@@H]1O","1"
+"C00073","Methionine","HMDB0000696","6137","C00073","CSCC[C@H](N)C(O)=O","1"
+"C00253","Nicotinic acid","HMDB0001488","938","C00253","OC(=O)C1=CN=CC=C1","1"
+"C01020","6-Hydroxynicotinic acid","HMDB0002658","72924","C01020","OC(=O)C1=CN=C(O)C=C1","1"
+"C01601","Pelargonic acid","HMDB0000847","8158","C01601","CCCCCCCCC(O)=O","1"
+"C00077","Ornithine","HMDB0000214","6262","C00077","NCCC[C@H](N)C(O)=O","1"
+"C01742","NA","NA","NA","NA","NA","0"
+"C00079","Phenylalanine","HMDB0000159","6140","C00079","N[C@@H](CC1=CC=CC=C1)C(O)=O","1"
+"C00009","Phosphate","HMDB0001429","1004","C00009","OP(O)(O)=O","1"
+"C00408","Pipecolic acid","HMDB0000070","849","C00408","OC(=O)C1CCCCN1","1"
+"C00148","Proline","HMDB0000162","145742","C00148","OC(=O)[C@@H]1CCCN1","1"
+"C01157","4-Hydroxyproline","HMDB0000725","5810","C01157","O[C@H]1CN[C@@H](C1)C(O)=O","1"
+"C06468","NA","NA","NA","NA","NA","0"
+"C00134","Putrescine","HMDB0001414","1045","C00134","NCCCCN","1"
+"C02502","2-Hydroxypyridine","HMDB0013751","8871","C02502","OC1=CC=CC=N1","1"
+"C01879","Pyroglutamic acid","HMDB0000267","7405","C01879","OC(=O)[C@@H]1CCC(=O)N1","1"
+"C00022","Pyruvic acid","HMDB0000243","1060","C00022","CC(=O)C(O)=O","1"
+"C00296","Quinic acid","HMDB0003072","6508","C00296","O[C@@H]1C[C@@](O)(C[C@@H](O)[C@H]1O)C(O)=O","1"
+"C00492","Raffinose","HMDB0003213","439242","C00492","OC[C@H]1O[C@@](CO)(O[C@H]2O[C@H](CO[C@H]3O[C@H](CO)[C@H](O)[C@H](O)[C@H]3O)[C@@H](O)[C@H](O)[C@H]2O)[C@@H](O)[C@@H]1O","1"
+"C00507","Rhamnose","HMDB0000849","25310","C00507","C[C@@H]1OC(O)[C@H](O)[C@H](O)[C@H]1O","1"
+"C00979","O-Acetylserine","HMDB0003011","99478","C00979","CC(=O)OC[C@H](N)C(O)=O","1"
+"C00065","Serine","HMDB0000187","5951","C00065","N[C@@H](CO)C(O)=O","1"
+"C00493","Shikimic acid","HMDB0003070","8742","C00493","O[C@@H]1CC(=C[C@@H](O)[C@H]1O)C(O)=O","1"
+"C00315","Spermidine","HMDB0001257","1102","C00315","NCCCCNCCCN","1"
+"C00042","Succinic acid","HMDB0000254","1110","C00042","OC(=O)CCC(O)=O","1"
+"C00089","Sucrose","HMDB0000258","5988","C00089","OC[C@H]1O[C@@](CO)(O[C@H]2O[C@H](CO)[C@@H](O)[C@H](O)[C@H]2O)[C@@H](O)[C@@H]1O","1"
+"C06424","Myristic acid","HMDB0000806","11005","C06424","CCCCCCCCCCCCCC(O)=O","1"
+"C01620","Threonic acid","HMDB0000943","5460407","C01620","OC[C@H](O)[C@@H](O)C(O)=O","1"
+"C00188","L-Threonine","HMDB0000167","6288","C00188","C[C@@H](O)[C@H](N)C(O)=O","1"
+"C01083","Trehalose","HMDB0000975","7427","C01083","OC[C@H]1O[C@H](O[C@H]2O[C@H](CO)[C@@H](O)[C@H](O)[C@H]2O)[C@H](O)[C@@H](O)[C@@H]1O","1"
+"C00078","L-Tryptophan","HMDB0000929","6305","C00078","N[C@@H](CC1=CNC2=C1C=CC=C2)C(O)=O","1"
+"C00483","Tyramine","HMDB0000306","5610","C00483","NCCC1=CC=C(O)C=C1","1"
+"C00082","L-Tyrosine","HMDB0000158","6057","C00082","N[C@@H](CC1=CC=C(O)C=C1)C(O)=O","1"
+"C05422","Dehydroascorbic acid","HMDB0001264","440667","C05422","[H][C@@]1(OC(=O)C(=O)C1=O)[C@@H](O)CO","1"
+"C00482","Sinapic acid","HMDB0032616","637775","C00482","COC1=CC(\C=C\C(O)=O)=CC(OC)=C1O","1"
+"C00106","Uracil","HMDB0000300","1174","C00106","O=C1NC=CC(=O)N1","1"
+"C00086","Urea","HMDB0000294","1176","C00086","NC(N)=O","1"
+"C00183","L-Valine","HMDB0000883","6287","C00183","CC(C)[C@H](N)C(O)=O","1"
+"C00181","D-Xylose","HMDB0000098","135191","C00181","O[C@@H]1COC(O)[C@H](O)[C@H]1O","1"
diff --git a/runs/Met_whole/prenorm.qs b/runs/Met_whole/prenorm.qs
new file mode 100644
index 0000000000000000000000000000000000000000..cfbff313d9aa44244057a5960528e7f3fda9a2a1
Binary files /dev/null and b/runs/Met_whole/prenorm.qs differ
diff --git a/runs/Met_whole/preproc.qs b/runs/Met_whole/preproc.qs
new file mode 100644
index 0000000000000000000000000000000000000000..cfbff313d9aa44244057a5960528e7f3fda9a2a1
Binary files /dev/null and b/runs/Met_whole/preproc.qs differ
diff --git a/runs/Met_whole/raw_dataview.csv b/runs/Met_whole/raw_dataview.csv
new file mode 100644
index 0000000000000000000000000000000000000000..d0b87c9ef0a4ec5fbebe3b27cb1d8dfacb5fe932
--- /dev/null
+++ b/runs/Met_whole/raw_dataview.csv
@@ -0,0 +1,51 @@
+"","aliquot_ID","genotype","group","C00147","C00020","C00956","C06104","C00041","C02512","C00099"
+"1","B17c","M82","green_leaf_M82",0.517550328266783,0.116297314385273,0.521565529098574,0.753867050446336,0.912389041443362,0.498391615702464,0.587732847033946
+"2","I05b","M82","green_leaf_M82",0.873197488141036,0.673488910735222,0.65133049828935,0.91454982019302,1.09606248301126,2.02933958337809,0.67000572421169
+"3","B05b","M82","green_leaf_M82",0.893863818037196,1.26603937004449,0.881668671004367,0.871399471536548,0.501880672476331,0.356120883168139,1.00876931959325
+"4","B17a","M82","green_leaf_M82",0.772333326603507,1.33401993755097,0.882331492905826,1.0129234705862,1.1120760489314,0.222983955208823,0.40236443457727
+"5","I06c","M82","green_leaf_M82",0.895739751006669,0.482208152992893,1.78263474524841,1.1734645281022,1.34804718755421,0.339719755337666,1.01419731819396
+"6","B06c","M82","green_leaf_M82",0.971279446981424,1.01449164952299,1.03253289474338,0.97692408074514,1.52166119301564,1.99717591134928,1.10362364221421
+"7","B06b","M82","green_leaf_M82",1.22139206948076,1.23024382217404,0.958739470886406,0.91222580911007,1.60465544060946,0.611845666274938,1.58021319970939
+"8","B11d","CANA1","green_leaf_CANA1",0.69605674587893,0.499725225926165,1.9392905628835,1.31336198336717,4.43005543951608,0.542494456871184,2.10304399387154
+"9","B14d","CANA1","green_leaf_CANA1",1.28296341670762,1.28447407502476,2.01194828743908,0.732881032419114,2.3918031379216,0.615742450303569,1.37280044279856
+"10","B12d","CANA1","green_leaf_CANA1",1.27000451234034,1.062606338981,2.37949215854601,1.03092785678314,3.37909302787288,1.04125619572873,1.56514262657594
+"11","B14a","CANA1","green_leaf_CANA1",0.707724435475351,0.834295885197712,0.935511810800953,1.45020749741446,2.91525725632371,1.22873072521149,1.90914374971176
+"12","B11c","CANA1","green_leaf_CANA1",0.885808267268499,0.632382949783313,0.958181099496651,1.56461542956252,4.04787074457592,0.721210335251307,1.61250812810363
+"13","B13d","CANA1","green_leaf_CANA1",1.10797805549254,0.73313059148844,2.05892840765958,0.897724636559396,2.19203314759952,1.29258322434074,1.53890222647633
+"14","I01c","CANA1","green_leaf_CANA1",1.19801089748378,1.95408109445323,1.60908994717651,0.923904120926657,3.38442487050826,0.547369766234931,1.59852928268639
+"15","2131598","M82","green_leaf_M82",0.285637304125361,0.198976739417061,0.680199804224252,0.726288802830615,0.357249250884837,0.441014925745586,0.414670336248597
+"16","2131600","M82","green_leaf_M82",0.130548405094935,0.62439029004948,0.489224559664597,1.40388809622427,0.317752149705409,0.730727539481527,2.62531432381745
+"17","2131596","M82","green_leaf_M82",0.285647485279867,0.290496753697845,0.73623828931482,1.03700941560192,1.34171635242358,1.4851255463651,1.06821983145755
+"18","2131599","M82","green_leaf_M82",0.611379867760303,2.67912911011042,3.5391936696865,2.76068612318482,1.27003553082329,3.99001036002593,3.21339813019689
+"19","2131597","M82","green_leaf_M82",2.95496142899016,0.618283091911722,1.36648501101274,1.12802001065491,1.31168014093896,1.5403867522058,0.916333574610609
+"20","2131609","CANA1","green_leaf_CANA1",1.23971093429929,0.962340463192019,0.270940678046416,4.55384429180132,0.149917378814419,13.637467964092,3.00960613589722
+"21","2131614","CANA1","green_leaf_CANA1",0.0160465378708477,0.153582005764058,1.39822933290782,1.18287684472437,2.41878578201321,8.70931463188714,2.23837410131661
+"22","2131610","CANA1","green_leaf_CANA1",0.337514679110327,0.28434147074013,0.620751590997593,0.891344145132462,2.47888276439447,2.92002389072043,1.77142589683725
+"23","2131620","CANA1","green_leaf_CANA1",0.269696886394736,0.0798384367918847,2.82141992314034,4.24096210827602,3.63923751210769,3.61708186204113,3.16054735156469
+"24","2131618","CANA1","green_leaf_CANA1",0.641097779732855,0.178517174009216,2.01573317415472,0.330825212877157,5.04533063598459,8.71872015415802,3.86739251117153
+"25","2131619","CANA1","green_leaf_CANA1",0.299403600988385,0.584013953826422,1.83722033682475,0.308217062211498,1.69128502457248,9.31406527453619,3.01036334920338
+"26","2131607","CANA1","green_leaf_CANA1",0.861752946359618,0.13264883768304,1.33301681447253,1.01399799342866,3.66971622704291,2.16735074348903,2.37462396264922
+"27","2131624","CANA1","green_leaf_CANA1",0.161479461188589,0.303197838781484,0.0494786874564523,1.66166242312318,1.90863447956696,5.02030626829994,1.46910366810027
+"28","2131621","CANA1","green_leaf_CANA1",0.341383738538429,1.15479348780799,0.0465358498044521,1.92228113434499,2.16165667343194,5.67219953334222,1.83255550575552
+"29","2131616","CANA1","green_leaf_CANA1",3.39902261006722,0.892248748225491,0.204342593894141,1.21995187357204,3.243161607097,4.76404302931545,3.15072465316875
+"30","2131617","CANA1","green_leaf_CANA1",0.263855394129024,1.28771446821961,1.05349043309771,1.135646426507,2.95926918273782,17.1181473364575,1.37131793645265
+"31","2131611","CANA1","green_leaf_CANA1",0.325341221692356,2.77728081929329,0.169648111394139,0.14225567818867,13.0955698535466,17.4140762294029,3.74631938400485
+"32","2131606","CANA1","green_leaf_CANA1",0.633935763645125,0.353007053368199,0.509979519583682,0.769123526007729,4.86113426568836,2.82100310837478,1.39908764264347
+"33","2131608","CANA1","green_leaf_CANA1",0.328514751201363,1.69641401464218,0.0828447874047525,1.9150328790548,0.0253197502828897,4.1606140176022,3.11020599745003
+"34","2131623","CANA1","green_leaf_CANA1",1.10646621332473,0.938525493665886,1.66346926998438,0.589695370850272,5.19163027209168,2.42444724353202,4.0559949613617
+"35","2131625","CANA1","green_leaf_CANA1",0.676219131267936,0.589074911245405,0.174692367740808,0.879796031517418,6.40757742019373,5.66335368751813,2.70077573781844
+"36","2131612","CANA1","green_leaf_CANA1",13.1490392620276,4.21086969571993,11.304734322256,9.71642593916786,10.6960638778318,12.4766938300322,7.23398217597849
+"37","2131622","CANA1","green_leaf_CANA1",0.0968303306794077,0.143726006204678,2.89475075180171,0.518999108606748,8.25328572236853,8.06846260105493,1.9740110123636
+"38","2227165","M82","green_leaf_M82",4.24810412745876,4.66381083257834,0.420773085670529,0.835830141796583,2.48362968134026,0.510146603295017,0.403557441755876
+"39","2227708","M82","green_leaf_M82",0.267587704444878,0.451080821857624,0.527385260568216,0.415260770286742,0.139730796228183,0.327851193472368,0.138028743315106
+"40","2227164","M82","green_leaf_M82",0.917113129179661,0.242945736173811,1.10498036191397,0.748032694502448,0.258577334092478,0.559139740866802,0.413715565228398
+"41","2227167","M82","green_leaf_M82",0.1536643191487,0.11409746679781,0.424716655768071,0.329629714198175,0.422856696521339,0.360019968122469,0.439855567835805
+"42","2227183","CANA1","green_leaf_CANA1",0.0573009532640558,1.93196007013482,1.21909963607767,0.474821598936868,0.561854615410136,0.724678386257927,1.30179090650046
+"43","2227172","CANA1","green_leaf_CANA1",1.11732576854196,0.58691021378005,1.48000155219887,0.020562833255018,0.393395933317842,11.1687133715111,1.02423836254977
+"44","2227177","CANA1","green_leaf_CANA1",0.245890718834525,1.0428653701969,0.221088431094138,0.714938518760093,0.227330674470048,1.88120187453601,0.39155791040523
+"45","2227176","CANA1","green_leaf_CANA1",0.492568923853148,0.821789617621,0.734091982758946,1.82614546685926,0.270591871737609,0.865329016682951,0.439786909372301
+"46","2227182","CANA1","green_leaf_CANA1",0.21390231442955,0.807063482511505,1.96661289283662,0.37342217500698,1.53171960407979,0.479772852860241,0.464766391973201
+"47","2227173","CANA1","green_leaf_CANA1",0.0172011140676406,0.414995789948616,0.0130756958741895,0.195816199292027,0.0508198103733565,0.404325773381279,0.262885116308655
+"48","2227174","CANA1","green_leaf_CANA1",0.474732332334937,0.595592179962751,2.40424548845806,0.209176946775961,5.39796343410305,0.0399933376787281,0.827465649307349
+"49","2227186","CANA1","green_leaf_CANA1",0.105494683300382,0.27956095951952,0.496406637076318,0.144465643272226,0.244777500815178,0.340212771603715,0.451476803166359
+"50","2227170","CANA1","green_leaf_CANA1",0.291795158345364,1.23892356532495,0.811418298964994,0.803567641282244,0.775606859823324,0.201954017803061,0.676292780152962
diff --git a/runs/Met_whole/row_norm.qs b/runs/Met_whole/row_norm.qs
new file mode 100644
index 0000000000000000000000000000000000000000..cfbff313d9aa44244057a5960528e7f3fda9a2a1
Binary files /dev/null and b/runs/Met_whole/row_norm.qs differ
diff --git a/runs/Met_whole/tosend.rds b/runs/Met_whole/tosend.rds
new file mode 100644
index 0000000000000000000000000000000000000000..504fabeb64efe845a4f5502e51691eb5bad7ab83
Binary files /dev/null and b/runs/Met_whole/tosend.rds differ
diff --git a/studies/canal-1_2017_GH/resources/Labels.xlsx b/studies/canal-1_2017_GH/resources/Labels.xlsx
new file mode 100644
index 0000000000000000000000000000000000000000..1c0731f37dbd77f73703aef46766b8cb37c3fe6a
Binary files /dev/null and b/studies/canal-1_2017_GH/resources/Labels.xlsx differ
diff --git a/workflows/Met_whole/Met_whole.R b/workflows/Met_whole/Met_whole.R
index 8d451cffa6853a33f6a98dee92f386cb23ebb6b5..be4ff2132d616baba4388260955ca334a7c47b99 100644
--- a/workflows/Met_whole/Met_whole.R
+++ b/workflows/Met_whole/Met_whole.R
@@ -11,6 +11,8 @@ library(viridis)
library(car)
library(broom)
library(here)
+library(ggrepel)
+library(openxlsx)
# Directory setting -------------------------------------------------------
@@ -49,7 +51,7 @@ metdat_GC_class <- readxl::read_xlsx(here(dat_in,
isa_tidy <- read_csv(here(prim_source, "isa_tidy.csv")) %>%
distinct(source_name_study, tissue_group, year, leaf_col, genotype, treatment,
- genotype_alt, tissue, group)
+ genotype_alt, tissue, group, sample_name_study)
fc_1_prim <- read_csv(here(prim_source, "mean_values_se_n.csv")) %>%
select(-Compound_Class) %>%
@@ -95,7 +97,8 @@ fc_1_ind_sec <- read_csv(here(sec_source, "individual_values.csv")) %>%
"tissue_group",
"sample_name_study")) %>%
left_join(plot_label, by = c("genotype_alt" = "geno",
- "tissue" = "tissue_label"))
+ "tissue" = "tissue_label")) %>%
+ rename(aliquot_ID = source_name_study)
sig_sec <- read_csv(here(sec_source, "p_values.csv"))
@@ -109,7 +112,15 @@ lip_met <- read_csv(here(lip_source, "metabolitedata_lipids.csv")) %>%
fc_1_lip <- read_csv(here(lip_source, "mean_values_se_n.csv")) %>%
left_join(lip_met)
-fc_1_ind_lip <- read_csv(here(lip_source, "individual_values.csv"))
+fc_1_ind_lip <- read_csv(here(lip_source, "individual_values.csv")) %>%
+ select(-group) %>%
+ left_join(isa_tidy_lc, by = c("tissue",
+ "treatment",
+ "genotype",
+ "tissue_group",
+ "aliquot_ID" = "source_name_study")) %>%
+ left_join(plot_label, by = c("genotype_alt" = "geno",
+ "tissue" = "tissue_label"))
sig_lip <- read_csv(here(lip_source, "p_values.csv"))
@@ -287,7 +298,7 @@ ggsave(plot = heat_all,
units = "cm",
dpi = 300)
-if(!is.null(dev.cur())){
+if(!dev.cur() == 1){
dev.off()
}
@@ -409,165 +420,176 @@ ggsave(plot = p1, here(out, "asparagine.png"),
saveRDS(p1, here(out, "asparagine_leaves.RDS"))
-# MetaboanalystR ----------------------------------------------------------
-
-fc_1 <- read_csv(here(analysis, "mean_values_se_n.csv")) %>%
- mutate(genotype = fct_relevel(genotype, wt))
-
-genotypes <- fc_1 %>%
- distinct(alias, genotype)
-
-fc_1_ind <- read_csv(here(analysis,"individual_values.csv")) %>%
- mutate(genotype = fct_relevel(genotype, wt))
-
-sig_GC <- read_csv(here(analysis,"p_values.csv")) #%>%
-
-maf <- readxl::read_xlsx("H:/3. cmQTL mapping/Shared_source_files/Primary_metabolites_classification_MWA.xlsx")
-
-met_dat <- fc_1 %>%
- distinct(met, Compound_Name, Compound_Class, component)
-
-genotypes <- fc_1 %>%
- distinct(alias, genotype)
-
-#think of way to select best features for analysis
-skip = T
-if(skip == T) {
-} else{
-
-}
-
-pathqea <- function(wt = "M82",
- geno_path,
- col_path,
- tiss_path,
- group_path){
-
- comp <- str_c(geno_path, col_path, tiss_path, sep = "_")
-
- sig_kegg <- fc_1 %>%
- separate(tissue, into = c("leaf_col", "tissue"), sep = 5) %>%
- mutate(group = as_factor(str_c(leaf_col, tissue, genotype, sep = "_"))) %>%
- ungroup() %>%
- filter(group %in% group_path, !is.na(component), treatment == "control") %>%
- #left_join(metdat_GC_class, by = c("component" = "Xcal_name_xreport")) #%>%
- pivot_wider(id_cols = c(KEGG_ID_ChEBI_mapped, met),
- names_from = group,
- values_from = mean_fc) %>%
- filter(!is.na(KEGG_ID_ChEBI_mapped)) %>%
- rowwise() %>%
- mutate(total_change = sum(abs(c_across(where(is.numeric))))) %>%
- group_by(KEGG_ID_ChEBI_mapped) %>%
- mutate(rank = rank(-total_change, ties.method = "first")) %>%
- filter(rank == 1)
-
- conc_tab <- fc_1_ind %>%
- separate(tissue, into = c("leaf_col", "tissue"), sep = 5) %>%
- select(-Compound_Name) %>%
- left_join(met_dat) %>%
- mutate(group = as_factor(str_c(leaf_col, tissue, genotype, sep = "_"))) %>%
- ungroup() %>%
- filter(group %in% group_path, !is.na(component), treatment == "control") %>%
- left_join(metdat_GC_class, by = c("component" = "Xcal_name_xreport")) %>%
- filter(!is.na(KEGG_ID_ChEBI_mapped), met %in% sig_kegg$met) %>%
- #group_by(KEGG_ID_ChEBI_mapped) %>%
- #mutate(rank_fc = -rank(abs(fc)),
- # rank_p = rank(p.value),
- # rank_c = rank_fc + rank_p,
- # rank_rank = rank(rank_c, ties.method = "average"),
- # rank_rank = if_else(rank_rank%%1 != 0, rank_p, rank_rank)) %>%
- #ungroup() %>%
- #filter(rank_rank == 1) %>%
- #distinct(KEGG_ID_ChEBI_mapped, LIMS_ID, tissue, treatment, genotype, .keep_all = T) %>%
- pivot_wider(id_cols = c(aliquot_ID, genotype, group),
- names_from = KEGG_ID_ChEBI_mapped,
- values_from = fc) %>%
- as.data.frame()
-
- write_csv(conc_tab, here(out, "ma_concentration_table.csv"))
-
- # Enrichment --------------------------------------------------------------
-
- mSet <- NULL
- # Create vector consisting of compounds for enrichment analysis
- #tmp.vec <- c("Acetoacetic acid", "Beta-Alanine", "Creatine", "Dimethylglycine", "Fumaric acid", "Glycine", "Homocysteine", "L-Cysteine", "L-Isolucine", "L-Phenylalanine", "L-Serine", "L-Threonine", "L-Tyrosine", "L-Valine", "Phenylpyruvic acid", "Propionic acid", "Pyruvic acid", "Sarcosine")
- # Create mSetObj for storing objects created during your analysis
- anal.type <- "pathqea"
- mSet<-InitDataObjects(data.type = "conc", anal.type = "pathqea", paired = FALSE)
- mSet <- Read.TextData(mSetObj = mSet, filePath = here(out, "ma_concentration_table.csv"), format = "rowu", lbl.type = "disc")
- # Set up mSetObj with the list of compounds
- #mSet<-Setup.MapData(mSet, conc_tab);
- # Cross reference list of compounds against libraries (hmdb, pubchem, chebi, kegg, metlin)
- mSet<-CrossReferencing(mSetObj = mSet, q.type = "kegg", lipid = F);
- # Creates a mapping result table; shows HMDB, KEGG, PubChem, etc. IDs
- # Saved as "name_map.csv" or can be found in mSet$dataSet$map.table
- # Compounds with no hits will contain NAs across the columns
- mSet<-CreateMappingResultTable(mSet);
- # From the mapping result table, L-Isolucine has no matches
- # Now, perform potential matching with our database against this compound
- #mSet<-PerformDetailMatch(mSet, "L-Isolucine");
- # Get list of candidates for matching
- # Results are found in mSet$name.map$hits.candidate.list
- #mSet<-GetCandidateList(mSet);
- # Replace L-Isolucine with selected compound (L-Isoleucine)
- #mSet<-SetCandidate(mSet, "L-Isolucine", "L-Isoleucine");
- # Select the pathway library, ranging from mammals to prokaryotes
- # Note the third parameter, where users need to input the KEGG pathway version.
-
- 1
- mSet<-SanityCheckData(mSet)
- mSet<-ReplaceMin(mSet)
- mSet<-PreparePrenormData(mSet)
- mSet<-Normalization(mSet, "NULL", "NULL", "NULL")
-
- # Use "current" for the latest KEGG pathway library or "v2018" for the KEGG pathway library version prior to November 2019.
- mSet<-SetKEGG.PathLib(mSet, "ath", "current")
- # Set the metabolite filter
- # Default set to false
- mSet<-SetMetabolomeFilter(mSet, T);
- # Calculate the over representation analysis score, here we selected to use the hypergeometric test (alternative is Fisher's exact test)
- # A results table "pathway_results.csv" will be created and found within your working directory
- api.base <- "http://api.xialab.ca"
- mSet<-CalculateQeaScore(mSetObj = mSet, nodeImp = "rbc", "gt")
- # Plot of the Pathway Analysis Overview
- #mSet<-PlotPathSummary(mSet, show.grid = T, "path_view_0_", "png", dpi = 72, width=NA)
- # Plot a specific metabolic pathway, in this case "Glycine, serine and threonine metabolism"
- #mSet<-PlotKEGGPath(mSet, "Glycine, serine and threonine metabolism",528, 480, "png")
-
- paths <- tibble("pathway ID" = mSet$analSet$path.ids,
- pathway = mSet$analSet$path.nms)
-
- mapped_cpds <- tibble("pathway ID" = names(mSet$analSet$qea.hits),
- hits = map(mSet$analSet$qea.hits,
- .f = ~tibble(Compound = names(.x),
- KEGG = .x))) %>%
- unnest(cols = c("pathway ID", hits)) %>%
- left_join(paths)
-
- path_res <- read_csv("pathway_results.csv",
- col_names = c("pathway",
- "total compounds",
- "hits", "p-values",
- "-log10(p)",
- "holm",
- "FDR",
- "impact"),
- skip = 1) %>%
- left_join(paths) %>%
- mutate(ratio = hits/`total compounds`,
- `-log10(p)` = -log10(`p-values`)) %>%
- select("pathway ID", "pathway", ratio, "total compounds", "hits", "p-values", "-log10(p)", "holm", "FDR", "impact")
-
- wb <- createWorkbook()
- addWorksheet(wb, "Pathway results")
- addWorksheet(wb, "Mapped Compounds")
- writeDataTable(wb, sheet = 1, x = path_res)
- writeDataTable(wb, sheet = 2, x = mapped_cpds)
- write.xlsx(list(path_res, mapped_cpds),
- here(out, str_c(comp, "pathway_results.xlsx", sep = "_")))
-}
-
-
+# # MetaboanalystR ----------------------------------------------------------
+# #currently defunct/browser utility used
+# #not recommended for further usage until stable
+# library(MetaboAnalystR)
+#
+# # fc_1 <- read_csv(here(analysis, "mean_values_se_n.csv")) %>%
+# # mutate(genotype = fct_relevel(genotype, wt))
+# #
+# # genotypes <- fc_1 %>%
+# # distinct(alias, genotype)
+# #
+# # fc_1_ind <- read_csv(here(analysis,"individual_values.csv")) %>%
+# # mutate(genotype = fct_relevel(genotype, wt))
+# #
+# # sig_GC <- read_csv(here(analysis,"p_values.csv")) #%>%
+# #
+# # maf <- readxl::read_xlsx("H:/3. cmQTL mapping/Shared_source_files/Primary_metabolites_classification_MWA.xlsx")
+# #
+# # met_dat <- fc_1 %>%
+# # distinct(met, Compound_Name, Compound_Class, component)
+# #
+# # genotypes <- fc_1 %>%
+# # distinct(alias, genotype)
+#
+# #think of way to select best features for analysis
+# skip = T
+# if(skip == T) {
+# } else{
+#
+# }
+#
+# #for test
+# geno_path <- "CANA1"
+# col_path <- c("green", "white")
+# tiss_path <- "leaf"
+# group_path <- c("green_leaf_CANA1", "green_leaf_M82")
+#
+# pathqea <- function(wt = "M82",
+# geno_path,
+# col_path,
+# tiss_path,
+# group_path){
+#
+# comp <- str_c(geno_path, col_path, tiss_path, sep = "_")
+#
+# sig_kegg <- fc_1 %>%
+# separate(tissue, into = c("leaf_col", "tissue"), sep = 5) %>%
+# mutate(group = as_factor(str_c(leaf_col, tissue, genotype, sep = "_"))) %>%
+# ungroup() %>%
+# filter(group %in% group_path, !is.na(component), treatment == "control") %>%
+# #left_join(metdat_GC_class, by = c("component" = "Xcal_name_xreport")) #%>%
+# pivot_wider(id_cols = c(KEGG_ID_ChEBI_mapped, met),
+# names_from = group,
+# values_from = mean_fc) %>%
+# filter(!is.na(KEGG_ID_ChEBI_mapped)) %>%
+# rowwise() %>%
+# mutate(total_change = sum(abs(c_across(where(is.numeric))))) %>%
+# group_by(KEGG_ID_ChEBI_mapped) %>%
+# mutate(rank = rank(-total_change, ties.method = "first")) %>%
+# filter(rank == 1)
+#
+# conc_tab <- fc_1_ind %>%
+# separate(tissue, into = c("leaf_col", "tissue"), sep = 5) %>%
+# select(-Compound_Name) %>%
+# left_join(met_dat) %>%
+# mutate(group = as_factor(str_c(leaf_col, tissue, genotype, sep = "_"))) %>%
+# ungroup() %>%
+# filter(group %in% group_path, !is.na(component), treatment == "control") %>%
+# left_join(metdat_GC_class, by = c("component" = "Xcal_name_xreport")) %>%
+# filter(!is.na(KEGG_ID_ChEBI_mapped), met %in% sig_kegg$met) %>%
+# #group_by(KEGG_ID_ChEBI_mapped) %>%
+# #mutate(rank_fc = -rank(abs(fc)),
+# # rank_p = rank(p.value),
+# # rank_c = rank_fc + rank_p,
+# # rank_rank = rank(rank_c, ties.method = "average"),
+# # rank_rank = if_else(rank_rank%%1 != 0, rank_p, rank_rank)) %>%
+# #ungroup() %>%
+# #filter(rank_rank == 1) %>%
+# #distinct(KEGG_ID_ChEBI_mapped, LIMS_ID, tissue, treatment, genotype, .keep_all = T) %>%
+# pivot_wider(id_cols = c(aliquot_ID, genotype, group),
+# names_from = KEGG_ID_ChEBI_mapped,
+# values_from = fc) %>%
+# as.data.frame()
+#
+# write_csv(conc_tab, here(out, "ma_concentration_table.csv"))
+#
+# # Enrichment --------------------------------------------------------------
+#
+# #setwd(here(out))
+#
+# mSet <- NULL
+# # Create vector consisting of compounds for enrichment analysis
+# #tmp.vec <- c("Acetoacetic acid", "Beta-Alanine", "Creatine", "Dimethylglycine", "Fumaric acid", "Glycine", "Homocysteine", "L-Cysteine", "L-Isolucine", "L-Phenylalanine", "L-Serine", "L-Threonine", "L-Tyrosine", "L-Valine", "Phenylpyruvic acid", "Propionic acid", "Pyruvic acid", "Sarcosine")
+# # Create mSetObj for storing objects created during your analysis
+# anal.type <- "pathqea"
+# mSet<-InitDataObjects(data.type = "conc", anal.type = "pathqea", paired = FALSE)
+# mSet <- Read.TextData(mSetObj = mSet, filePath = here(out, "ma_concentration_table.csv"), format = "rowu", lbl.type = "disc")
+# # Set up mSetObj with the list of compounds
+# #mSet<-Setup.MapData(mSet, conc_tab);
+# # Cross reference list of compounds against libraries (hmdb, pubchem, chebi, kegg, metlin)
+# mSet<-CrossReferencing(mSetObj = mSet, q.type = "kegg", lipid = F);
+# # Creates a mapping result table; shows HMDB, KEGG, PubChem, etc. IDs
+# # Saved as "name_map.csv" or can be found in mSet$dataSet$map.table
+# # Compounds with no hits will contain NAs across the columns
+# mSet<-CreateMappingResultTable(mSet);
+# # From the mapping result table, L-Isolucine has no matches
+# # Now, perform potential matching with our database against this compound
+# #mSet<-PerformDetailMatch(mSet, "L-Isolucine");
+# # Get list of candidates for matching
+# # Results are found in mSet$name.map$hits.candidate.list
+# #mSet<-GetCandidateList(mSet);
+# # Replace L-Isolucine with selected compound (L-Isoleucine)
+# #mSet<-SetCandidate(mSet, "L-Isolucine", "L-Isoleucine");
+# # Select the pathway library, ranging from mammals to prokaryotes
+# # Note the third parameter, where users need to input the KEGG pathway version.
+#
+# 1
+# mSet<-SanityCheckData(mSet)
+# mSet<-ReplaceMin(mSet)
+# mSet<-PreparePrenormData(mSet)
+# mSet<-Normalization(mSet, "NULL", "NULL", "NULL")
+#
+# # Use "current" for the latest KEGG pathway library or "v2018" for the KEGG pathway library version prior to November 2019.
+# mSet<-SetKEGG.PathLib(mSet, "ath", "current")
+# # Set the metabolite filter
+# # Default set to false
+# mSet<-SetMetabolomeFilter(mSet, T);
+# # Calculate the over representation analysis score, here we selected to use the hypergeometric test (alternative is Fisher's exact test)
+# # A results table "pathway_results.csv" will be created and found within your working directory
+# api.base <- "http://api.xialab.ca"
+# mSet<-CalculateQeaScore(mSetObj = mSet, nodeImp = "rbc", "gt")
+# # Plot of the Pathway Analysis Overview
+# #mSet<-PlotPathSummary(mSet, show.grid = T, "path_view_0_", "png", dpi = 72, width=NA)
+# # Plot a specific metabolic pathway, in this case "Glycine, serine and threonine metabolism"
+# #mSet<-PlotKEGGPath(mSet, "Glycine, serine and threonine metabolism",528, 480, "png")
+#
+# paths <- tibble("pathway ID" = mSet$analSet$path.ids,
+# pathway = mSet$analSet$path.nms)
+#
+# mapped_cpds <- tibble("pathway ID" = names(mSet$analSet$qea.hits),
+# hits = map(mSet$analSet$qea.hits,
+# .f = ~tibble(Compound = names(.x),
+# KEGG = .x))) %>%
+# unnest(cols = c("pathway ID", hits)) %>%
+# left_join(paths)
+#
+# path_res <- read_csv("pathway_results.csv",
+# col_names = c("pathway",
+# "total compounds",
+# "hits", "p-values",
+# "-log10(p)",
+# "holm",
+# "FDR",
+# "impact"),
+# skip = 1) %>%
+# left_join(paths) %>%
+# mutate(ratio = hits/`total compounds`,
+# `-log10(p)` = -log10(`p-values`)) %>%
+# select("pathway ID", "pathway", ratio, "total compounds", "hits", "p-values", "-log10(p)", "holm", "FDR", "impact")
+#
+# wb <- createWorkbook()
+# addWorksheet(wb, "Pathway results")
+# addWorksheet(wb, "Mapped Compounds")
+# writeDataTable(wb, sheet = 1, x = path_res)
+# writeDataTable(wb, sheet = 2, x = mapped_cpds)
+# write.xlsx(list(path_res, mapped_cpds),
+# here(out, str_c(comp, "pathway_results.xlsx", sep = "_")))
+# }
+#
+#
# Metaboanalyst bubble plots ---------------------------------------------------
path_tidy <- readxl::read_xlsx(here(dat_in, "ma_pathway_abbreviations.xlsx"))
@@ -601,12 +623,10 @@ bub_theme <- theme(axis.text.x = element_markdown(vjust = 2, hjust = 0.5, size =
# WTglB_vs_MUwl -----------------------------------------------------------
cont <- "WTglB_vs_MUwl"
-bub <- read_csv("H:/2. CANA1, SCO2/CANA1_GC_analysis/230130_analysis/cana1_wl_pathway_results.csv",
+bub <- read_csv(here(out, "cana1_wl_pathway_results.csv"),
col_names = c("pathway", "total compounds", "hits", "p-values", "-log10(p)", "holm", "FDR", "impact"),
skip = 1)
-setwd(out_dir)
-
bub_tidy <- bub %>%
left_join(path_tidy) %>%
rename(path_old = pathway,
@@ -630,23 +650,18 @@ bub_tidy %>%
#theme(legend.position = "")
scale_fill_continuous(type = "viridis")
-ggsave(str_c("Metaboanalyst_Bubble_", cont, ".png"),
+ggsave(here(out, str_c("Metaboanalyst_Bubble_", cont, ".png")),
width = 5.5, height = 10, units = "cm")
-write.xlsx(bub_tidy, str_c("H:/2. CANA1, SCO2/Manuscript_CANA1_2022/Tables/Metaboanalyst_bubble_", cont, ".xlsx"))
+write.xlsx(bub_tidy, here(out, str_c("Metaboanalyst_bubble_", cont, ".xlsx")))
# WTglB_vs_MUgl -----------------------------------------------------------
-setwd(current)
-
-readxl::read_xlsx("221026_MUglB_vs_WTgl_ChemRICH_results.xlsx")
cont <- "WTglB_vs_MUgl"
-bub <- read_csv("H:/2. CANA1, SCO2/CANA1_GC_analysis/230130_analysis/cana1_gl_pathway_results.csv",
+bub <- read_csv(here(out,"cana1_gl_pathway_results.csv"),
col_names = c("pathway", "total compounds", "hits", "p-values", "-log10(p)", "holm", "FDR", "impact"),
skip = 1)
-setwd(out_dir)
-
bub_tidy <- bub %>%
left_join(path_tidy) %>%
rename(path_old = pathway,
@@ -670,10 +685,10 @@ bub_tidy %>%
#theme(legend.position = "")
scale_fill_continuous(type = "viridis")
-ggsave(str_c("Metaboanalyst_Bubble_", cont, ".png"),
+ggsave(here(out, str_c("Metaboanalyst_Bubble_", cont, ".png")),
width = 5.5, height = 10, units = "cm")
-write.xlsx(bub_tidy, str_c("H:/2. CANA1, SCO2/Manuscript_CANA1_2022/Tables/Metaboanalyst_bubble_", cont, ".xlsx"))
+write.xlsx(bub_tidy, here(out, str_c("Metaboanalyst_bubble_", cont, ".xlsx")))
## Combined ----------------------------------------------------------------
@@ -700,383 +715,246 @@ bub_comb %>%
scale_fill_continuous(type = "viridis") +
ylim(c(0,7))
-saveRDS(last_plot(),"H:/2. CANA1, SCO2/Manuscript_CANA1_2022/Figures/RDS_files/Combined_bubble_metabolites.RDS")
+saveRDS(last_plot(), here(out,
+ "Combined_bubble_metabolites.RDS"))
-ggsave(str_c("Bubble_", "comb", ".png"),
+ggsave(here(out, str_c("Bubble_", "comb", ".png")),
width = 16.5, height = 10, units = "cm")
-# Pathview ----------------------------------------------------------------
-
-
-######### --------------- try double pathview ####
-library(pathview)
-setwd("H:/2. CANA1, SCO2/CANA1_RNAseq_experiment/CANA1_RNAseq_analysis")
-allfiles_allgenes <- list.files(pattern = "allgenes.csv")
-dates_allfiles <- as.numeric(str_extract(allfiles_allgenes, pattern = "^\\d{6}"))
-cur_date <- as.numeric(str_replace_all(Sys.Date(),"^.{2}|-",""))
-latest_allgenes <- allfiles_allgenes[[which.min(cur_date - dates_allfiles)]]
-
-all_genes <- read_csv(latest_allgenes)
-sampleTable <- read.csv("CANA1_2017_RNAseq_samplelist.csv")
-
-all_genes_long <- all_genes %>%
- dplyr::select(-contains("log2"),
- -contains("padj"),) %>%
- pivot_longer(cols = c(contains("Mu"),
- contains ("Wt")),
- names_to = "geno",
- values_to = "exp") %>%
- filter(is.na(exp)==F) %>%
- mutate(geno = as_factor(geno)) %>%
- separate(col = geno,
- into = c("geno", "rep"),
- sep = -1) %>%
- mutate(geno_tissue = as_factor(geno)) %>%
- separate(col = geno,
- into = c("geno","tissue"),
- sep=2) %>%
- mutate(tissue = as_factor(tissue),
- geno = as_factor(geno)) %>%
- dplyr::rename(gene = Row.names)
-
-entrez_solyc_ID_paper <-read_csv("solyc_vs_entrez_unique_original.csv")
-
-degs <- all_genes %>%
- dplyr::select(Row.names, contains("padj"), contains("Log2FC")) %>%
- mutate(DEG_WTglB_vs_MUwl = if_else(padj_WTglB_vs_MUwl <= 0.1 & abs(log2FC_WTglB_vs_MUwl) >= 1, T,F),
- DEG_WTglB_vs_MUgl = if_else(padj_WTglB_vs_MUgl <= 0.1 & abs(log2FC_WTglB_vs_MUgl) >= 1, T,F),
- DEG_MUglB_vs_MUwl = if_else(padj_MUglB_vs_MUwl <= 0.1 & abs(log2FC_MUglB_vs_MUwl) >= 1, T,F),
- DEG_WTsB_vs_MUs = if_else(padj_WTsB_vs_MUs <= 0.1 & abs(log2FC_WTsB_vs_MUs) >= 1, T,F),
- Row.names = str_replace_all(Row.names, pattern = "(\\.\\d){2,3}", replacement = "")) %>%
- rename(SolycID = Row.names) %>%
- left_join(entrez_solyc_ID_paper)
-
-
-pv_leaf_genes <- degs %>%
- select(log2FC_WTglB_vs_MUwl, DEG_WTglB_vs_MUwl,log2FC_WTglB_vs_MUgl, DEG_WTglB_vs_MUgl,
- SolycID, EntrezID) %>%
- filter(!is.na(EntrezID)) %>%
- #mutate(log2FC_WTglB_vs_MUwl = if_else(DEG_WTglB_vs_MUwl == T, log2FC_WTglB_vs_MUwl, 0),
- # log2FC_WTglB_vs_MUgl = if_else(DEG_WTglB_vs_MUgl == T, log2FC_WTglB_vs_MUgl, 0)) %>%
- select(EntrezID, log2FC_WTglB_vs_MUgl, log2FC_WTglB_vs_MUwl) %>%
- as.data.frame()
-row.names(pv_leaf_genes) <- pv_leaf_genes$EntrezID
-pv_leaf_genes_plot <- pv_leaf_genes[,-1]
+# Heatmap SAD----------------------------------------------------------------
+
+library(pheatmap)
+plot_label <- read.xlsx("studies/canal-1_2017_GH/resources/Labels.xlsx")
+red_met_class <- read_csv(here(dat_in, "reduced_met_classes.csv")) %>%
+ select(-n)
-setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
-getwd()
-#test pathview, need distinct occurences of KEGG_IDs
-setwd(out_dir)
-#metdat_GC_class <- readxl::read_xlsx("H:/3. cmQTL mapping/Shared_source_files/210118_primary_metabolites_classification.xlsx")
+heat_base <- fc_1_ind %>%
+ filter(!str_detect(genotype, "on"),
+ !str_detect(treatment, "drought")) %>%
+ select(- Compound_Name,
+ - geno_tissue,
+ - geno_label,
+ - geno_tissue_label,
+ - genotype_alt) %>%
+ left_join(met_dat) %>%
+ left_join(plot_label, by = c("genotype" = "geno_old",
+ "tissue" = "tissue_label")) %>%
+ group_by(met,
+ Compound_Name,
+ Compound_Class,
+ tissue,
+ treatment,
+ genotype,
+ geno_label,
+ tissue_group,
+ geno_tissue_label) %>%
+ mutate(sad = (abs(fc - median(fc)))/median(fc)) %>%
+ summarise(sad = mean(sad)) %>%
+ ungroup()
-sig_kegg <- fc_1 %>%
- filter(tissue_group == "leaves", genotype %in% c("M82", "CANA1")) %>%
- mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_")),
- lfc = log2(mean_fc)) %>%
+sad_M82 <- heat_base %>%
+ filter(genotype == "M82") %>%
+ group_by(met,
+ Compound_Name,
+ Compound_Class,
+ tissue,
+ treatment,
+ genotype,
+ geno_label,
+ tissue_group,
+ geno_tissue_label) %>%
+ summarise(sad_M82 = mean(sad)) %>%
ungroup() %>%
- filter(treatment == "control") %>%
- #left_join(metdat_GC_class, by = c("component" = "Xcal_name_xreport")) %>%
- pivot_wider(id_cols = c(KEGG_ID_ChEBI_mapped, met),
- names_from = group,
- values_from = lfc) %>%
- filter(!is.na(KEGG_ID_ChEBI_mapped)) %>%
- rowwise() %>%
- mutate(total_change = sum(abs(c_across(where(is.numeric))))) %>%
- group_by(KEGG_ID_ChEBI_mapped) %>%
- mutate(rank = rank(-total_change, ties.method = "first")) %>%
- filter(rank == 1)
-
-
-pv_leaf_compound <- fc_1 %>%
- mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_")),
- lfc = log2(mean_fc)) %>%
- #mutate(lfc = if_else(signif == "*", lfc, 0)) %>%
- filter(tissue_group == "leaves", genotype == "CANA1") %>%
- #group_by(tissue, met,stress) %>%
+ select(-genotype, -geno_label, -geno_tissue_label, -tissue)
+
+heat_nom <- heat_base %>%
+ left_join(sad_M82) %>%
+ mutate(mean_fc = sad/sad_M82) %>%
+ filter(!is.na(geno_tissue_label)) %>%
+ group_by(Compound_Name, met) %>%
+ mutate(log_norm = log2(mean_fc),
+ log_norm = if_else(is.infinite(log_norm), 0, log_norm),
+ log_norm_level = (log_norm - mean(log_norm))/(max(log_norm)-min(log_norm)),
+ z_score = (mean_fc - mean(mean_fc))/(max(mean_fc)-min(mean_fc))) %>%
ungroup() %>%
- filter(treatment == "control"
- , met %in% sig_kegg$met
- ) %>%
- #left_join(metdat_GC_class, by = c("component" = "Xcal_name_xreport")) %>%
- pivot_wider(id_cols = c(KEGG_ID_ChEBI_mapped, met, genotype),
- names_from = group,
- values_from = lfc) %>%
- filter(!is.na(KEGG_ID_ChEBI_mapped)) %>%
- #rowwise() %>%
- #mutate(total_change = sum(abs(c_across(where(is.numeric))))) %>%
- #group_by(KEGG_ID_ChEBI_mapped) %>%
- #mutate(rank = rank(-total_change, ties.method = "first")) %>%
- #filter(rank == 1) %>%
- #select(-met, - genotype, -total_change, -rank) %>%
+ mutate(group = as_factor(str_c(tissue, genotype, treatment, sep = "_"))) %>%
+ pivot_wider(id_cols = c(Compound_Name, Compound_Class, met),
+ names_from = geno_tissue_label,
+ values_from = log_norm) %>%
+ left_join(red_met_class) %>%
+ select(-Compound_Class, Compound_Class = Compound_Class_dense) %>%
+ arrange(Compound_Class, Compound_Name) %>%
as.data.frame()
-row.names(pv_leaf_compound) <- pv_leaf_compound$KEGG_ID_ChEBI_mapped
-pv_leaf_compound_plot <- pv_leaf_compound[,-(1:3)]
-
-
-pathview(gene.data = pv_leaf_genes_plot[,1:2],
- cpd.data = pv_leaf_compound_plot[,1:2],
- species = "sly",
- pathway.id = "00290",
- kegg.native = T,
- low = list(gene = "magenta", cpd = "blue"),
- mid = list(gene = "gray", cpd ="gray"),
- high = list(gene = "green", cpd = "yellow"),
- limit = list(gene = 2, cpd = 4),
- bins = list(gene = 10, cpd = 10),
- out.suffix = "val_leu_ile_leaves")
-
-pathview(gene.data = pv_leaf_genes_plot[,1:2],
- cpd.data = pv_leaf_compound_plot[,1:2],
- species = "sly",
- pathway.id = "00250",
- kegg.native = T,
- low = list(gene = "magenta", cpd = "blue"),
- mid = list(gene = "gray", cpd ="gray"),
- high = list(gene = "green", cpd = "yellow"),
- limit = list(gene = 2, cpd = 4),
- bins = list(gene = 10, cpd = 10),
- out.suffix = "ala_asp_asn_leaves")
-
-pathview(gene.data = pv_leaf_genes_plot[,1:2],
- cpd.data = pv_leaf_compound_plot[,1:2],
- species = "sly",
- pathway.id = "00195",
- kegg.native = T,
- low = list(gene = "magenta", cpd = "blue"),
- mid = list(gene = "gray", cpd ="gray"),
- high = list(gene = "green", cpd = "yellow"),
- limit = list(gene = 2, cpd = 4),
- bins = list(gene = 10, cpd = 10),
- out.suffix = "photosynthesis_leaves")
-
-
-pathview(gene.data = pv_leaf_genes_plot[,1:2],
- cpd.data = pv_leaf_compound_plot[,1:2],
- species = "sly",
- pathway.id = "00196",
- kegg.native = T,
- low = list(gene = "magenta", cpd = "blue"),
- mid = list(gene = "gray", cpd ="gray"),
- high = list(gene = "green", cpd = "yellow"),
- limit = list(gene = 2, cpd = 4),
- bins = list(gene = 10, cpd = 10),
- out.suffix = "photosynthesis_antenna_leaves")
-
-pathview(gene.data = pv_leaf_genes_plot[,1:2],
- cpd.data = pv_leaf_compound_plot[,1:2],
- species = "sly",
- pathway.id = "00500",
- kegg.native = T,
- low = list(gene = "magenta", cpd = "blue"),
- mid = list(gene = "gray", cpd ="gray"),
- high = list(gene = "green", cpd = "yellow"),
- limit = list(gene = 2, cpd = 4),
- bins = list(gene = 10, cpd = 10),
- out.suffix = "starch_suc_leaves")
-
-pathview(gene.data = pv_leaf_genes_plot[,1:2],
- cpd.data = pv_leaf_compound_plot[,1:2],
- species = "sly",
- pathway.id = "00710",
- kegg.native = T,
- low = list(gene = "magenta", cpd = "blue"),
- mid = list(gene = "gray", cpd ="gray"),
- high = list(gene = "green", cpd = "yellow"),
- limit = list(gene = 2, cpd = 4),
- bins = list(gene = 10, cpd = 10),
- out.suffix = "met_path_leaves")
-
-pathview(gene.data = pv_leaf_genes_plot[,1:2],
- cpd.data = pv_leaf_compound_plot[,1:2],
- species = "sly",
- pathway.id = "00592",
- kegg.native = T,
- low = list(gene = "magenta", cpd = "blue"),
- mid = list(gene = "gray", cpd ="gray"),
- high = list(gene = "green", cpd = "yellow"),
- limit = list(gene = 2, cpd = 4),
- bins = list(gene = 10, cpd = 10),
- out.suffix = "alpha_lin_acid_leaves")
-
-pathview(gene.data = pv_leaf_genes_plot[,1:2],
- cpd.data = pv_leaf_compound_plot[,1:2],
- species = "sly",
- pathway.id = "00620",
- kegg.native = T,
- low = list(gene = "magenta", cpd = "blue"),
- mid = list(gene = "gray", cpd ="gray"),
- high = list(gene = "green", cpd = "yellow"),
- limit = list(gene = 2, cpd = 4),
- bins = list(gene = 10, cpd = 10),
- out.suffix = "pyruvate_leaves")
-
-pathview(gene.data = pv_leaf_genes_plot[,1:2],
- cpd.data = pv_leaf_compound_plot[,1:2],
- species = "sly",
- pathway.id = "00020",
- kegg.native = T,
- low = list(gene = "magenta", cpd = "blue"),
- mid = list(gene = "gray", cpd ="gray"),
- high = list(gene = "green", cpd = "yellow"),
- limit = list(gene = 2, cpd = 4),
- bins = list(gene = 10, cpd = 10),
- out.suffix = "tca_leaves")
-
-pathview(gene.data = pv_leaf_genes_plot[,1:2],
- cpd.data = pv_leaf_compound_plot[,1:2],
- species = "sly",
- pathway.id = "00010",
- kegg.native = T,
- low = list(gene = "magenta", cpd = "blue"),
- mid = list(gene = "gray", cpd ="gray"),
- high = list(gene = "green", cpd = "yellow"),
- limit = list(gene = 2, cpd = 4),
- bins = list(gene = 10, cpd = 10),
- out.suffix = "glycolysis_leaves")
-
-pathview(gene.data = pv_leaf_genes_plot[,1:2],
- cpd.data = pv_leaf_compound_plot[,1:2],
- species = "sly",
- pathway.id = "00280",
- kegg.native = T,
- low = list(gene = "magenta", cpd = "blue"),
- mid = list(gene = "gray", cpd ="gray"),
- high = list(gene = "green", cpd = "yellow"),
- limit = list(gene = 2, cpd = 4),
- bins = list(gene = 10, cpd = 10),
- out.suffix = "val_ile_leu_deg")
-
-pathview(gene.data = pv_leaf_genes_plot[,1:2],
- cpd.data = pv_leaf_compound_plot[,1:2],
- species = "sly",
- pathway.id = "00630",
- kegg.native = T,
- low = list(gene = "magenta", cpd = "blue"),
- mid = list(gene = "gray", cpd ="gray"),
- high = list(gene = "green", cpd = "yellow"),
- limit = list(gene = 2, cpd = 4),
- bins = list(gene = 10, cpd = 10),
- out.suffix = "glyoxylate")
-
-pathview(gene.data = pv_leaf_genes_plot[,1:2],
- cpd.data = pv_leaf_compound_plot[,1:2],
- species = "sly",
- pathway.id = "00030",
- kegg.native = T,
- low = list(gene = "magenta", cpd = "blue"),
- mid = list(gene = "gray", cpd ="gray"),
- high = list(gene = "green", cpd = "yellow"),
- limit = list(gene = 2, cpd = 4),
- bins = list(gene = 10, cpd = 10),
- out.suffix = "PPP")
-
-
-# SBGN View ---------------------------------------------------------------
-library(SBGNview)
-data("sbgn.xmls")
-data("pathways.info", "pathways.stats")
-data("mapped.ids")
-
-gene.data <- as.matrix(pv_leaf_genes_plot)
-cpd.data <- as.matrix(pv_leaf_compound_plot)
-
-SBGN_gene <- changeDataId(data.input.id = gene.data,
- input.type = "entrez",
- output.type = "KO",
- mol.type = "gene",
- sum.method = "sum",
- org = "sly")
-
-SBGNview(gene.data = SBGN_gene[,1:2],
- cpd.data = cpd.data[,1:2],
- org = "sly",
- input.sbgn = "PWY-1042",
- output.file = "Glyc_test_MetaCyc",
- gene.id.type = "KO",
- cpd.id.type = "kegg",
- output.formats = "svg")
-
-# Test own layout ---------------------------------------------------------
-
-setwd(current)
-
-#invertase <- load("Invertase_test_file.sbgn")
-
-gene <- tribble(~gene, ~mugl, ~muwl,
- "Inv", -1, 1,
- "GAPDH", 3,4) %>%
- mutate(mugl = as.numeric(mugl),
- muwl = as.numeric(muwl)) %>%
- as.data.frame()
+rownames(heat_nom) <- heat_nom$met
-rownames(gene) <- gene$gene
-gene <- gene[,-1]
-gene <- as.matrix(gene)
+mat.heat_nom <- heat_nom %>%
+ select(#contains("flowers"), contains("roots"),
+ #contains("*canal-1*"),
+ #contains("M82")
+ contains("fr")
+ ) %>% as.matrix()
+annotation_row <- heat_nom %>%
+ select(Compound_Class)
-cpd <- tribble(~cpd, ~mugl, ~muwl,
- "Sucrose", "2", "1",
- "Glucose", "1", "3",
- "Fructose", "4", "0") %>%
+rownames(annotation_row) <- heat_nom$met
+
+annotation_col <- heat_base %>%
+ distinct(tissue, treatment, geno_label, geno_tissue_label) %>%
+ #mutate(group = as_factor(str_c(tissue, genotype, treatment, sep = "_"))) %>%
+ filter(geno_tissue_label %in% colnames(mat.heat_nom)) %>%
+ select(tissue, genotype = geno_label, geno_tissue_label) %>%
as.data.frame()
-rownames(cpd) <- cpd$cpd
-
-cpd <- as.matrix(cpd)
-
-#id.mapping.gene_tbl <- tribble(~SYMBOL, ~ENTREZID,
-# "Inv", "Invertase",
-# "prevent", "character_vector") %>%
-# as.data.frame()
-
-id.mapping.gene <- as.matrix(id.mapping.gene_tbl)
-
-SBGNview(gene.data = gene[,1:2],
- #cpd.data = cpd[,2:3],
- #org = "sly",
- input.sbgn = "Invertase_test_file_manually_modified.sbgn",
- sbgn.id.attr = "id",
- sbgn.gene.id.type = "test",
- #id.mapping.gene = id.mapping.gene,
- output.file = "Invertase_test",
- gene.id.type = "test",
- #cpd.id.type = "kegg",
- output.formats = "svg",)
-
-#load("230221_Figures/SBGNview.tmp.data/ath_ENTREZID_KO.RData")
-#load("230221_Figures/SBGNview.tmp.data/kegg_metacyc.SBGN.RData")
-
-id.mapping.gene.table <- entrez_solyc_ID_paper %>%
- select(entrez_in = EntrezID) %>%
- mutate(entrez_out = entrez_in,
- #entrez_in = str_c(entrez_in, "...")
- ) %>%
- as.data.frame() %>%
- as.matrix()
-
-#Vanted SBGN
-setwd(current)
-
-pv_leaf_genes_matrix <- pv_leaf_genes_plot %>% as.matrix()
-
-SBGNview(gene.data = pv_leaf_genes_matrix[,1:2],
- #cpd.data = pv_leaf_compound_plot[,1:2],
- org = "sly",
- input.sbgn = "Carbon fixation in photosynthetic organisms.sbgn",
- sbgn.id.attr = "label",
- sbgn.gene.id.type = "entrez_in",
- #id.mapping.gene = id.mapping.gene.table,
- #sbgn.cpd.id.type = "KEGG",
- output.file = "VANTED_transformed_KEGG",
- gene.id.type = "entrez_in",
- #cpd.id.type = "kegg",
- output.formats = "svg")
-
-load("230221_Figures/SBGNview.tmp.data/ath_ENTREZID_KO.RData")
-load("230221_Figures/SBGNview.tmp.data/kegg_metacyc.SBGN.RData")
+rownames(annotation_col) <- annotation_col$geno_tissue_label
+
+annotation_col <- annotation_col %>%
+ select(-geno_tissue_label, tissue, genotype)
+
+heat_cols <- colnames(mat.heat_nom)
+
+ann_colors = list(
+ tissue = c('fruit' = "red",
+ '#greenleaf' = "darkgreen", 'whiteleaf' = "beige",
+ #'roots' = "grey",
+ 'stem' = "lightgreen"
+ #,flowers = "yellow"
+ ),
+ genotype = c(M82 = cb_scale[1], `*canal-1*` = cb_scale[2]),
+ #treatment = c("control" = vir_scale[1], "drought" = vir_scale[3]),
+ Compound_Class = c("Amino Acid or derivative" = "#0072B2",
+ "Carboxylic Acid" = "#D55E00",
+ "Carbohydrate or derivative" = "#CC79A7",
+ "Cinnamic Acid" = "#484E37",
+ "Dipeptide" = "#4AEE2F",
+ "Flavonoid (glycosides)" = "#BAD23A",
+ "Galactolipid" = "#EF000B",
+ "Phospholipid" = "#19605B",
+ "Steroidal Glycoalkaloids" = "#7A5DF0",
+ "Triacylglyceride" = "#6E3455",
+ "Other" = "#000000")
+)
+
+pheatmap.all <- pheatmap(mat.heat_nom,
+ colorRampPalette(c("#440154FF", "white", "#FDE725FF"))(25),
+ #cellwidth = 16,
+ #cellheight = 4,
+ breaks = seq(-6.25, 6.25, 0.5),
+ #clustering_distance_rows = dist((mat.heat_nom), method = "euclidean"),
+ cluster_rows = F,
+ cluster_cols = F,
+ annotation_names_row = F,
+ show_rownames = F,
+ annotation_row = annotation_row,
+ annotation_col = annotation_col,
+ #display_numbers = mat.heat_nom_signif,
+ number_color = "black",
+ fontsize_number = 6,
+ fontsize = 6,
+ annotation_colors = ann_colors,
+ #filename = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+ # "cmQTL_val_1_heatmap_rel_tissue_wt.jpg",
+ # sep = "_")
+)
+
+dev.off()
+heat_all <- ggplotify::as.ggplot(pheatmap.all)
+saveRDS(heat_all, here(out, "Heatmap_all_met.RDS"))
+
+ggsave(plot = heat_all,
+ here(out, "Heatmap_fruit_sad.jpg"),
+ width = 16,
+ height = 10,
+ units = "cm",
+ dpi = 300)
+
+if(!dev.cur() == 1){
+ dev.off()
+}
+
+
+# Correlation to yield data? ----------------------------------------------
+
+yield <- read.xlsx(here("assays/canal-1_Phenotyping_2019/dataset/canal-1_yield_data.xlsx"))
+
+isa_fruit <- isa_tidy %>%
+ filter(!str_detect(genotype, "on"), ! year == 2017) %>%
+ distinct(source_name_study, genotype, genotype_alt, treatment)
+
+yield_sad <- yield %>%
+ filter(fruit_yield != 0) %>%
+ mutate(aliquot_ID = as.character(aliquot_ID)) %>%
+ left_join(isa_fruit, by = c("aliquot_ID" = "source_name_study",
+ "treatment")) %>%
+ group_by(treatment, genotype, genotype_alt) %>%
+ mutate(yield_sad = (abs(fruit_yield - median(fruit_yield))/median(fruit_yield)))
+
+fruit_met_sad <- fc_1_ind %>%
+ filter(!str_detect(genotype, "on"),
+ treatment == "control", tissue == "fruit") %>%
+ mutate(genotype = as_factor(genotype)) %>%
+ group_by(tissue, treatment, genotype, met) %>%
+ mutate(sad = (abs(fc - median(fc)))/median(fc)) %>%
+ ungroup()
+
+fruit_met_sad_wide <- fruit_met_sad %>%
+ pivot_wider(id_cols = c(tissue, treatment, genotype, aliquot_ID, sample_name_study),
+ names_from = met,
+ values_from = sad)
+
+# Test metabolite variation -----------------------------------------------
+
+library(rstatix)
+
+sad <- fc_1_ind %>%
+ filter(!str_detect(genotype, "on"),
+ treatment == "control", tissue == "fruit") %>%
+ mutate(genotype = as_factor(genotype)) %>%
+ group_by(tissue, treatment, genotype, met) %>%
+ mutate(sad = (abs(fc - median(fc)))/median(fc)) %>%
+ ungroup()
+
+signif <- sad %>%
+ group_by(tissue, treatment, met) %>%
+ wilcox_test(formula = sad ~ genotype,
+ ref.group = "M82",
+ p.adjust.method = "fdr")
+
+p1 <- sad %>%
+ filter(met == "m_86", str_detect(tissue, "fruit"), treatment == "control", !str_detect(genotype, "on")) %>%
+ #mutate(genotype = as_factor(geno_tissue_label),
+ # genotype = fct_relevel(genotype, c("gl M82", "gl *canal-1*", "wl *canal-1*"))) %>%
+ ggplot(aes(x = genotype, y = sad)) +
+ geom_boxplot() +
+ geom_dotplot()#+
+ #geom_errorbar(aes(ymin = (mean_fc-se), ymax = (mean_fc + se)), position = position_dodge(0.9), width = 0.25, linewidth = 0.75)+
+ stat_pvalue_manual(step.group.by = c("treatment", "tissue_group"), sig_comb_bar, label = "p.signif", y.position = "y.position",
+ step.increase = 0.07, tip.length = 0.01,
+ hide.ns = F) +
+ # theme(axis.text.x = element_markdown(angle = 45, hjust = 1),
+ # panel.background = element_rect(fill = "white"),
+ # panel.border = element_rect(color = "black",fill = NA),
+ # text = element_text(size = 14),
+ # legend.title = element_blank(),
+ # legend.text = element_markdown()) +
+ com_theme +
+ xlab("") +
+ ylab("Mean fold-change") +
+ #facet_grid(rows = vars(treatment), cols = vars(tissue_group)) +
+ ggtitle(label = "asparagine") +
+ scale_fill_manual(values = c("darkgreen", "beige"), aesthetics = "fill") +
+ scale_y_log10()
+
+p1
+
+ggsave(plot = p1, here(out, "asparagine.png"),
+ width = 16, height = 10, units = "cm")
+
+
+saveRDS(p1, here(out, "asparagine_leaves.RDS"))
diff --git a/workflows/Pigment_SpecPhot_analysis/Pigment_SpecPhot_Yield_analysis.R b/workflows/Pigment_SpecPhot_analysis/Pigment_SpecPhot_Yield_analysis.R
new file mode 100644
index 0000000000000000000000000000000000000000..04be5ce1993cd25909efee669a2e2c0eb1cfa7aa
--- /dev/null
+++ b/workflows/Pigment_SpecPhot_analysis/Pigment_SpecPhot_Yield_analysis.R
@@ -0,0 +1,731 @@
+library(tidyverse)
+library(car)
+library(broom)
+library(ggtext)
+library(ggpubr)
+library(openxlsx)
+library(ggbeeswarm)
+
+sampledata <- readxl::read_excel("H:/2. CANA1, SCO2/CANA1_2018_2019_samplelist.xlsx",
+ sheet = 2)
+extraction <- readxl::read_excel("H:/2. CANA1, SCO2/CANA1_2018_2019_samplelist.xlsx",
+ sheet = 3)
+treatment <- readxl::read_excel("H:/2. CANA1, SCO2/CANA1_2018_2019_samplelist.xlsx",
+ sheet = 4)
+
+av_yield <- treatment %>%
+ filter(Name != "no_LIMS") %>%
+ mutate(fruit_yield = as.numeric(fruit_yield)) %>%
+ group_by(genotype,trtmnt) %>%
+ summarise(av_yield = mean(fruit_yield, na.rm = T),
+ sd_yield = sd(fruit_yield, na.rm = T),
+ obs = n()) %>%
+ mutate(cv_yield = sd_yield/av_yield)
+
+iCV <- av_yield %>%
+ pivot_wider(id_cols = c(genotype),
+ names_from = trtmnt,
+ values_from = cv_yield) %>%
+ mutate(iCV = drought/control)
+
+treatment <- treatment %>%
+ filter(Name != "no_LIMS") %>%
+ mutate(genotype = as_factor(genotype),
+ trtmnt = as_factor(trtmnt))
+
+yield_control <- treatment %>%
+ filter(trtmnt == "control")
+
+yield_drought <- treatment %>%
+ filter(trtmnt == "drought")
+
+plot_label <- tribble(~genotype, ~group,
+ "M82", "M82",
+ "CANA1", "*canal-1*")
+
+y_for_p <- treatment %>%
+ group_by(trtmnt) %>%
+ summarise(y.pos = 1.1*max(fruit_yield))
+
+p_ctrl <- tibble(p_value = kruskal.test(fruit_yield ~ genotype , data = yield_control)$p.value) %>%
+ mutate(trtmnt = "control")
+
+p_drought <- tibble(p_value = kruskal.test(fruit_yield ~ genotype , data = yield_drought)$p.value) %>%
+ mutate(trtmnt = "drought")
+
+yield_p <- bind_rows(p_ctrl, p_drought)
+
+sig_yield_plot <- y_for_p %>%
+ left_join(yield_p) %>%
+ mutate(group1 = "M82",
+ group2 = "*canal-1*",
+ group = group1,
+ p.signif = if_else(p_value <= 0.0005, "***",
+ if_else(p_value <= 0.005, "**",
+ if_else(p_value <= 0.05, "*", "ns")))) %>%
+ mutate(y = "fruit_yield")
+
+
+# Theme -------------------------------------------------------------------
+
+
+com_theme <- theme(axis.text.x = element_markdown(angle = 45, hjust = 1, size = 6, family = "sans"),
+ axis.text.y = element_text(size = 6, family = "sans"),
+ axis.title.x = element_blank(),
+ axis.title.y = element_text(size = 6, family = "sans"),
+ panel.background = element_rect(fill = "white"),
+ panel.border = element_rect(color = "black",fill = NA),
+ strip.text = element_text(size = 8, family = "sans", margin = margin(t = 1, r = 1, b = 1, l = 1 , unit = "pt")),
+ text = element_text(size = 6, family = "sans"),
+ legend.title = element_blank(),
+ legend.text = element_markdown(size = 6),
+ plot.margin = unit(c(1,1,1,1), "mm"),
+ legend.margin = margin(t = 0, r = 2, b = 0, l = 2 , unit = "mm"))
+
+
+# Yield plot --------------------------------------------------------------
+
+
+treatment <- treatment %>%
+ left_join(plot_label)
+
+binwidth <- treatment %>%
+ summarise(binwidth = (max(fruit_yield) - min(fruit_yield))/50)
+
+treatment %>%
+ mutate(group = as_factor(group),
+ trtmnt = as_factor(trtmnt)) %>%
+ ggplot(aes(x = group, y = fruit_yield, fill = group)) +
+ geom_boxplot(size = 0.2, outlier.size = 0.2) +
+ geom_quasirandom(method = "quasirandom", shape = 21, size = 0.5) +
+ #geom_dotplot(binaxis = "y", stackdir = "center", binwidth = binwidth$binwidth) +
+ stat_pvalue_manual(data = sig_yield_plot, y.position = "y.pos", label = "p.signif", xmin = "group1", xmax = "group2") +
+ facet_grid(cols = vars(trtmnt), scales = "free") +
+ scale_fill_manual(values = c("darkgreen", "beige")) +
+ labs(x = "Genotype", y = "Fruit yield/ g/plant") +
+ com_theme
+
+ggsave(str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"_fruit_yield.jpg"),
+ width = 89,
+ height = 60,
+ units = "mm",
+ dpi = 300)
+
+saveRDS(last_plot(), "H:/2. CANA1, SCO2/Manuscript_CANA1_2022/Figures/RDS_files/canal1_yield.RDS")
+
+sum_stat <- treatment %>%
+ group_by(trtmnt, group) %>%
+ summarise(mean = mean(fruit_yield),
+ sd = sd(fruit_yield),
+ cv = sd/mean,
+ n = n()) %>%
+ ungroup() %>%
+ mutate(group = str_remove_all(group, "\\*"),
+ mean = round(mean, digits = 4),
+ sd = round(sd, digits = 4),
+ cv = round(cv, digits = 4)) %>%
+ select(Treatment = trtmnt, Genotype = group, Mean = mean, SD = sd, CV = cv, N = n) %>%
+ mutate(Genotype = as_factor(Genotype),
+ Genotype = fct_relevel(Genotype, c("M82", "canal-1"))) %>%
+ arrange(Treatment, Genotype)
+
+write.xlsx(sum_stat, "H:/2. CANA1, SCO2/Manuscript_CANA1_2022/Tables/Yield_sum_stat.xlsx")
+
+# Test MAD ------------------------------------------------------
+
+mad <- treatment %>%
+ group_by(trtmnt, genotype) %>%
+ mutate(mad = abs(fruit_yield - median(fruit_yield))/median(fruit_yield)) %>%
+ ungroup()
+
+mad_nest <- mad %>%
+ group_by(trtmnt) %>%
+ nest()
+
+mad_p <- mad_nest %>%
+ mutate(krusk = map(data, .f = ~kruskal.test(.x$mad ~ .x$genotype)$p.value)) %>%
+ select(trtmnt, krusk)
+
+y_for_p <- mad %>%
+ group_by(trtmnt) %>%
+ summarise(y.pos = 1.1*max(mad))
+
+sig_yield_plot <- mad_p %>%
+ left_join(y_for_p) %>%
+ mutate(group1 = "M82",
+ group2 = "*canal-1*",
+ group = group1,
+ p.signif = if_else(krusk <= 0.0005, "***",
+ if_else(krusk <= 0.005, "**",
+ if_else(krusk <= 0.05, "*", "ns")))) %>%
+ mutate(y = "mad")
+
+binwidth <- mad %>%
+ summarise(binwidth = (max(mad) - min(mad))/50)
+
+mad %>%
+ mutate(group = as_factor(group),
+ trtmnt = as_factor(trtmnt)) %>%
+ ggplot(aes(x = group, y = mad, fill = group)) +
+ geom_boxplot(size = 0.2, outlier.size = 0.2) +
+ geom_quasirandom(method = "quasirandom", shape = 21, size = 0.5) +
+ #geom_dotplot(binaxis = "y", stackdir = "center", binwidth = binwidth$binwidth) +
+ stat_pvalue_manual(data = sig_yield_plot, hide.ns = T, y.position = "y.pos", label = "p.signif") +
+ facet_grid(cols = vars(trtmnt), scales = "free") +
+ scale_fill_manual(values = c("darkgreen", "beige")) +
+ labs(x = "Genotype", y = "Scaled MAD Fruit yield/ g/plant") +
+ com_theme
+
+ggsave(str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"_fruit_yield_MAD.jpg"),
+ width = 89,
+ height = 60,
+ units = "mm",
+ dpi = 300)
+
+saveRDS(last_plot(), "H:/2. CANA1, SCO2/Manuscript_CANA1_2022/Figures/RDS_files/canal1_yield_mad_scaled.RDS")
+
+sum_stat <- mad %>%
+ group_by(trtmnt, group) %>%
+ summarise(mean = mean(mad),
+ sd = sd(mad),
+ cv = sd/mean,
+ n = n()) %>%
+ ungroup() %>%
+ mutate(group = str_remove_all(group, "\\*"),
+ mean = round(mean, digits = 4),
+ sd = round(sd, digits = 4),
+ cv = round(cv, digits = 4)) %>%
+ select(Treatment = trtmnt, Genotype = group, Mean = mean, SD = sd, N = n) %>%
+ mutate(Genotype = as_factor(Genotype),
+ Genotype = fct_relevel(Genotype, c("M82", "canal-1"))) %>%
+ arrange(Treatment, Genotype)
+
+write.xlsx(sum_stat, "H:/2. CANA1, SCO2/Manuscript_CANA1_2022/Tables/Yield_MAD_sum_stat.xlsx")
+
+
+# Test Levene ------------------------------------------------------
+skip = T
+if(skip == T){
+ print("skipped")
+} else{
+lev <- treatment %>%
+ group_by(trtmnt, genotype) %>%
+ mutate(lev = abs(log10(fruit_yield) - log10(median(fruit_yield)))) %>%
+ ungroup()
+
+lev_nest <- lev %>%
+ group_by(trtmnt) %>%
+ nest()
+
+lev_p <- lev_nest %>%
+ mutate(krusk = map(data, .f = ~kruskal.test(.x$lev ~ .x$genotype)$p.value)) %>%
+ select(trtmnt, krusk)
+
+y_for_p <- lev %>%
+ group_by(trtmnt) %>%
+ summarise(y.pos = 1.1*max(lev))
+
+sig_yield_plot <- lev_p %>%
+ left_join(y_for_p) %>%
+ mutate(group1 = "M82",
+ group2 = "*canal-1*",
+ group = group1,
+ p.signif = if_else(krusk <= 0.0005, "***",
+ if_else(krusk <= 0.005, "**",
+ if_else(krusk <= 0.05, "*", "ns")))) %>%
+ mutate(y = "lev")
+
+binwidth <- lev %>%
+ summarise(binwidth = (max(lev) - min(lev))/50)
+
+lev %>%
+ mutate(group = as_factor(group),
+ trtmnt = as_factor(trtmnt)) %>%
+ ggplot(aes(x = group, y = lev, fill = group)) +
+ geom_boxplot(size = 0.2, outlier.size = 0.2) +
+ geom_dotplot(binaxis = "y", stackdir = "center", binwidth = binwidth$binwidth) +
+ stat_pvalue_manual(data = sig_yield_plot, hide.ns = T, y.position = "y.pos", label = "p.signif") +
+ facet_grid(cols = vars(trtmnt), scales = "free") +
+ scale_fill_manual(values = c("darkgreen", "beige")) +
+ labs(x = "Genotype", y = "Scaled lev Fruit yield/ g/plant") +
+ com_theme
+
+ggsave(str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"_fruit_yield_lev.jpg"),
+ width = 89,
+ height = 60,
+ units = "mm",
+ dpi = 300)
+
+saveRDS(last_plot(), "H:/2. CANA1, SCO2/Manuscript_CANA1_2022/Figures/RDS_files/canal1_yield_lev_scaled.RDS")
+
+sum_stat <- lev %>%
+ group_by(trtmnt, group) %>%
+ summarise(mean = mean(lev),
+ sd = sd(lev),
+ cv = sd/mean,
+ n = n()) %>%
+ ungroup() %>%
+ mutate(group = str_remove_all(group, "\\*"),
+ mean = round(mean, digits = 4),
+ sd = round(sd, digits = 4),
+ cv = round(cv, digits = 4)) %>%
+ select(Treatment = trtmnt, Genotype = group, Mean = mean, SD = sd, N = n) %>%
+ mutate(Genotype = as_factor(Genotype),
+ Genotype = fct_relevel(Genotype, c("M82", "canal-1"))) %>%
+ arrange(Treatment, Genotype)
+
+write.xlsx(sum_stat, "H:/2. CANA1, SCO2/Manuscript_CANA1_2022/Tables/Yield_lev_sum_stat.xlsx")
+}
+
+# Chlorophyll -------------------------------------------------------------
+
+chl1 <- readxl::read_excel("210128_Chl_CANA1_1_to_10.xlsx",
+ sheet = 2)
+
+chl1_reformat <- chl1 %>%
+ pivot_longer(cols = matches("^\\d{1,2}$"),
+ names_to = "column",
+ values_to = "value") %>%
+ pivot_wider(id_cols = c(row,column),
+ names_from = identity,
+ values_from = value)
+
+chl2 <- readxl::read_excel("210130_Chl_CANA1_1_to_10.xlsx",
+ sheet = 2)
+
+chl2_reformat <- chl2 %>%
+ pivot_longer(cols = matches("^\\d{1,2}$"),
+ names_to = "column",
+ values_to = "value") %>%
+ pivot_wider(id_cols = c(row,column),
+ names_from = identity,
+ values_from = value)
+
+chl3 <- readxl::read_excel("210131_Chl_CANA1_1_to_10.xlsx",
+ sheet = 2)
+
+chl3_reformat <- chl3 %>%
+ pivot_longer(cols = matches("^\\d{1,2}$"),
+ names_to = "column",
+ values_to = "value") %>%
+ pivot_wider(id_cols = c(row,column),
+ names_from = identity,
+ values_from = value)
+
+chl_all <- chl1_reformat %>%
+ bind_rows(chl2_reformat,chl3_reformat) %>%
+ filter(!is.na(extraction_num))
+
+blank <- chl_all %>%
+ filter(extraction_num == 0) %>%
+ summarise(blank_652 = mean(`652`),
+ blank_665 = mean(`665`),
+ blank_470 = mean(`470`))
+
+sampledata_tidy <- sampledata %>%
+ left_join(extraction) %>%
+ left_join(treatment, by = "Aliquot ID") %>%
+ mutate(Genotype = as_factor(Genotype))
+
+write_csv(sampledata_tidy, str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"sampledata_tidy.csv"))
+
+levels(sampledata_tidy$Genotype)
+sampledata_tidy$Genotype <- ordered(sampledata_tidy$Genotype,
+ levels = c("M82", "CANA1", "CRISPR_8619", "CRISPR_7472"))
+
+chl_calc <- chl_all %>%
+ left_join(sampledata_tidy) %>%
+ filter(trtmnt == "control") %>%
+ mutate(path_652 = (`652`- blank$blank_652)/0.51,
+ path_665 = (`665`- blank$blank_665)/0.51,
+ path_470 = (`470`- blank$blank_470)/0.51,
+ chl_a = ((16.72 * path_665 - 9.16 * path_652)),
+ chl_b = ((34.09 * path_652 - 15.28 * path_665)),
+ car_tot = (((1000 * path_470 - 1.63 * chl_a - 104.96 * chl_b) / 221)),
+ chla_ug_g = if_else(chl_a < 0 , 0, (chl_a*2.5/fw) *1000),
+ chlb_ug_g = if_else(chl_b < 0 , 0, (chl_b*2.5/fw) *1000),
+ car_tot_ug_g = if_else(car_tot < 0, 0, (car_tot*2.5/fw) *1000),
+ ratio_chla_chlb = chla_ug_g/chlb_ug_g,
+ #ratio_chla_chlb = if_else(is.finite(ratio_chla_chlb), ratio_chla_chlb, 1),
+ Tissue = as_factor(Tissue),
+ tissue_group = as_factor(str_extract(Tissue, "(?<!\\w)\\w+$")),
+ tissue_short = abbreviate(Tissue, minlength = 2),
+ group_anova = as_factor(str_c(tissue_short,group, sep =" "))) %>%
+ filter(is.finite(ratio_chla_chlb))
+
+levels(chl_calc$group_anova)
+
+chl_calc$group_anova <- fct_relevel(chl_calc$group_anova,
+ c("gl M82",
+ "gl *canal-1*",
+ "wl *canal-1*",
+ "st M82",
+ "st *canal-1*"))
+
+chl_calc %>%
+ filter(tissue_group == "stems") %>%
+ ggplot(aes(x = Tissue, y = chla_ug_g, fill = genotype)) +
+ geom_dotplot(binaxis = "y", stackdir = "center", position = "dodge")
+
+chl_plot <- chl_calc %>%
+ pivot_longer(cols = c(contains("ug_g"), ratio_chla_chlb),
+ names_to = "pigment",
+ values_to = "concentration") %>%
+ mutate(pigment = as_factor(pigment))
+
+sum_stats <- chl_plot %>%
+ group_by(pigment, Tissue, genotype) %>%
+ summarise(norm_dist = shapiro.test(concentration)$p.value,
+ mean = mean(concentration),
+ sd = sd(concentration),
+ n = n()) %>%
+ mutate(se = sd/sqrt(n))
+
+data_leaves <- chl_calc %>%
+ filter(tissue_group == "leaves") %>%
+ select(group_anova, chla_ug_g, chlb_ug_g,car_tot_ug_g,ratio_chla_chlb)
+
+numeric_data_leaves <- data_leaves %>%
+ select(where(is.numeric))
+
+map_dfc(.x = numeric_data_leaves, .f = ~leveneTest(.x ~ data_leaves$group_anova))
+map_dfc(.x = numeric_data_leaves, .f = ~kruskal.test(.x ~ data_leaves$group_anova)$p.value)
+leaves_wilcox <- map(.x = numeric_data_leaves, .f = ~pairwise.wilcox.test(.x, data_leaves$group_anova))
+
+sig_leaves <- leaves_wilcox %>%
+ map(.f = tidy) %>%
+ map2(.y = names(leaves_wilcox), .f = ~.x %>% mutate(var = .y)) %>%
+ purrr::reduce(bind_rows)
+
+data_stems <- chl_calc %>%
+ filter(tissue_group == "stems") %>%
+ select(group_anova, chla_ug_g, chlb_ug_g,car_tot_ug_g, ratio_chla_chlb)
+
+numeric_data_stems <- data_stems %>%
+ select(where(is.numeric))
+
+map_dfc(.x = numeric_data_stems, .f = ~leveneTest(.x ~ data_stems$group_anova))
+map_dfc(.x = numeric_data_stems, .f = ~kruskal.test(.x ~ data_stems$group_anova)$p.value)
+stems_t.test <- map(.x = numeric_data_stems, .f = ~pairwise.t.test(.x, data_stems$group_anova))
+
+sig_stems <- stems_t.test %>%
+ map(.f = tidy) %>%
+ map2(.y = names(stems_t.test), .f = ~.x %>% mutate(var = .y)) %>%
+ purrr::reduce(bind_rows)
+
+#single tests
+#leveneTest(chla_ug_g ~ group_anova, data = data_chl_a_leaves)
+#kruskal.test(chla_ug_g ~ group_anova, data = data_chl_a_leaves)$p.value
+#pairwise.wilcox.test(data_chl_a_leaves$chla_ug_g,
+ # data_chl_a_leaves$group_anova)
+
+pigments <- tibble(
+ chla_ug_g = "chlorophyll a",
+ chlb_ug_g = "chlorophyll b",
+ car_tot_ug_g = "total carotenoids",
+ ratio_chla_chlb = "chlorophyll a/b ratio") %>%
+ pivot_longer(cols = everything(),
+ names_to = "pigment",
+ values_to = "pigment_label")
+
+
+# Leaves pigment plot -----------------------------------------------------
+
+chl_plot_leaves <- chl_plot %>%
+ filter(pigment != "ratio_chla_chlb",
+ tissue_group == "leaves") %>%
+ mutate(group = group_anova) %>%
+ left_join(pigments)
+
+tissues <- chl_plot_leaves %>%
+ distinct(group, Tissue)
+
+y_for_p <- chl_plot_leaves %>%
+ group_by(pigment) %>%
+ summarise(y.position = 1.1*max(concentration))
+
+sig_leaves_plot <- sig_leaves %>%
+ left_join(tissues, by = c("group2" = "group")) %>%
+ mutate(p.signif = if_else(p.value <= 0.0005, "***",
+ if_else(p.value <= 0.005, "**",
+ if_else(p.value <= 0.05, "*", "ns")))) %>%
+ left_join(y_for_p, by = c("var" = "pigment")) %>%
+ left_join(pigments, by = c("var" = "pigment")) %>%
+ filter(var != "ratio_chla_chlb")
+
+binwidth <- chl_plot_leaves %>%
+ summarise(binwidth = (max(concentration) - min(concentration))/100) %>%
+ ungroup()
+
+chl_plot_leaves %>%
+ ggplot(aes(x = group, y = concentration, fill = Tissue)) +
+ geom_boxplot(size = 0.2, outlier.size = 0.2) +
+ #geom_beeswarm(aes(x = group, y = concentration, fill = Tissue),
+ # shape = 21) +
+ geom_quasirandom(method = "quasirandom", shape = 21, size = 0.5) +
+ stat_pvalue_manual(sig_leaves_plot, y.position = "y.position",
+ label = "p.signif", step.increase = 0.07,
+ hide.ns = T) +
+ facet_grid(cols = vars(pigment_label), scales = "free") +
+ scale_fill_manual(values = c("darkgreen", "beige")) +
+ labs(x = "Tissue and genotype", y = "concentration/ ?g/g FW") +
+ com_theme
+
+ggsave(str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"_leaf_pigments.jpg"),
+ width = 183,
+ height = 100,
+ units = "mm",
+ dpi = 300)
+
+saveRDS(last_plot(), "H:/2. CANA1, SCO2/Manuscript_CANA1_2022/Figures/RDS_files/canal1_leaf_pigments.RDS")
+
+sum_stat <- chl_plot_leaves %>%
+ group_by(group, pigment_label) %>%
+ summarise(mean = mean(concentration),
+ sd = sd(concentration),
+ cv = sd/mean,
+ n = n()) %>%
+ ungroup() %>%
+ mutate(group = str_remove_all(group, "\\*"),
+ mean = round(mean, digits = 4),
+ sd = round(sd, digits = 4),
+ cv = round(cv, digits = 4)) %>%
+ select(Pigment = pigment_label, `Genotype and Tissue` = group, Mean = mean, SD = sd, N = n) %>%
+ mutate(`Genotype and Tissue` = as_factor(`Genotype and Tissue`),
+ `Genotype and Tissue` = fct_relevel(`Genotype and Tissue`, c("gl M82", "gl canal-1", "wl canal-1"))) %>%
+ arrange(Pigment, `Genotype and Tissue`)
+
+write.xlsx(sum_stat, "H:/2. CANA1, SCO2/Manuscript_CANA1_2022/Tables/Leaf_pigments_sum_stat.xlsx")
+
+
+# stems pigment plot -----------------------------------------------------
+
+chl_plot_stems <- chl_plot %>%
+ filter(pigment != "ratio_chla_chlb",
+ tissue_group == "stems") %>%
+ mutate(group = group_anova) %>%
+ left_join(pigments)
+
+tissues <- chl_plot_stems %>%
+ distinct(group, Tissue)
+
+y_for_p <- chl_plot_stems %>%
+ group_by(pigment) %>%
+ summarise(y.position = 1.1*max(concentration))
+
+sig_stems_plot <- sig_stems %>%
+ left_join(tissues, by = c("group2" = "group")) %>%
+ mutate(p.signif = if_else(p.value <= 0.0005, "***",
+ if_else(p.value <= 0.005, "**",
+ if_else(p.value <= 0.05, "*", "ns")))) %>%
+ left_join(y_for_p, by = c("var" = "pigment")) %>%
+ left_join(pigments, by = c("var" = "pigment")) %>%
+ filter(var != "ratio_chla_chlb") %>%
+ mutate(group = group1)
+
+binwidth <- chl_plot_stems %>%
+ summarise(binwidth = (max(concentration) - min(concentration))/100) %>%
+ ungroup()
+
+chl_plot_stems %>%
+ ggplot(aes(x = group, y = concentration, fill = group)) +
+ geom_boxplot(size = 0.2, outlier.size = NA, outlier.shape = NA) +
+ #geom_beeswarm(aes(x = group, y = concentration, fill = group),
+ # shape = 21) +
+ geom_quasirandom(method = "quasirandom", shape = 21, size = 0.5) +
+ stat_pvalue_manual(sig_stems_plot, y.position = "y.position",
+ label = "p.signif", step.increase = 0.07,
+ hide.ns = T) +
+ facet_grid(cols = vars(pigment_label), scales = "free") +
+ scale_fill_manual(values = c("darkgreen", "#a0c45a63")) +
+ labs(x = "Tissue and genotype", y = "concentration/ ?g/g FW") +
+ com_theme
+
+ggsave(str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"_stem_pigments.jpg"),
+ width = 183,
+ height = 100,
+ units = "mm",
+ dpi = 300)
+
+saveRDS(last_plot(), "H:/2. CANA1, SCO2/Manuscript_CANA1_2022/Figures/RDS_files/canal1_stem_pigments.RDS")
+
+sum_stat <- chl_plot_stems %>%
+ group_by(group, pigment_label) %>%
+ summarise(mean = mean(concentration),
+ sd = sd(concentration),
+ cv = sd/mean,
+ n = n()) %>%
+ ungroup() %>%
+ mutate(group = str_remove_all(group, "\\*"),
+ mean = round(mean, digits = 4),
+ sd = round(sd, digits = 4),
+ cv = round(cv, digits = 4)) %>%
+ select(Pigment = pigment_label, `Genotype and Tissue` = group, Mean = mean, SD = sd, N = n) %>%
+ mutate(`Genotype and Tissue` = as_factor(`Genotype and Tissue`),
+ `Genotype and Tissue` = fct_relevel(`Genotype and Tissue`, c("st M82", "st canal-1"))) %>%
+ arrange(Pigment, `Genotype and Tissue`)
+
+write.xlsx(sum_stat, "H:/2. CANA1, SCO2/Manuscript_CANA1_2022/Tables/Stems_pigments_sum_stat.xlsx")
+
+
+# Comparison pigments apolar LC -------------------------------------------
+
+pig_LC <- read_csv("pigments_apolar_LC.csv") %>%
+ left_join(plot_label) %>%
+ mutate(tissue_group = as_factor(str_extract(tissue, "(?<!\\w)\\w+$")),
+ tissue_short = abbreviate(tissue, minlength = 2),
+ group_anova = as_factor(str_c(tissue_short,group, sep =" ")))
+
+pigments_LC <- tibble(
+ Chlorophyll_a = "chlorophyll a",
+ Chlorophyll_b = "chlorophyll b",
+ car_tot_ug_g = "total carotenoids",
+ ratio_chla_chlb = "chlorophyll a/b ratio") %>%
+ pivot_longer(cols = everything(),
+ names_to = "Compound_Name",
+ values_to = "pigment_label")
+
+chl_calc_LC <- pig_LC %>%
+ filter(str_detect(Compound_Name, "Chlorophyll"),
+ tissue_group == "leaves" | tissue_group == "stems",
+ year == "2018" | year == "2019") %>%
+ mutate(fw_norm = area/fw_in_mg) %>%
+ select(Compound_Name, genotype, Tissue = tissue, aliquot_ID, tissue_group, fw_norm, group_anova) %>%
+ left_join(pigments_LC) %>%
+ pivot_wider(id_cols = c(genotype, Tissue, aliquot_ID, tissue_group, group_anova),
+ names_from = pigment_label,
+ values_from = fw_norm)
+
+levels(chl_calc_LC$group_anova)
+
+chl_calc_LC$group_anova <- fct_relevel(chl_calc_LC$group_anova,
+ c("gl M82",
+ "gl *canal-1*",
+ "wl *canal-1*",
+ "st M82",
+ "st *canal-1*"))
+
+chl_calc_LC %>%
+ filter(tissue_group == "stems") %>%
+ ggplot(aes(x = Tissue, y = `chlorophyll a`, fill = genotype)) +
+ geom_dotplot(binaxis = "y", stackdir = "center", position = "dodge")
+
+chl_plot_LC <- chl_calc_LC %>%
+ pivot_longer(cols = c(contains("chloro")),
+ names_to = "pigment",
+ values_to = "concentration") %>%
+ mutate(pigment = as_factor(pigment))
+
+sum_stats_LC <- chl_plot_LC %>%
+ group_by(pigment, Tissue, genotype) %>%
+ summarise(norm_dist = shapiro.test(concentration)$p.value,
+ mean = mean(concentration),
+ sd = sd(concentration),
+ n = n()) %>%
+ mutate(se = sd/sqrt(n))
+
+data_leaves_LC <- chl_calc_LC %>%
+ filter(tissue_group == "leaves") %>%
+ select(group_anova, `chlorophyll a`, `chlorophyll b`)
+
+numeric_data_leaves <- data_leaves %>%
+ select(where(is.numeric))
+
+map_dfc(.x = numeric_data_leaves, .f = ~leveneTest(.x ~ data_leaves$group_anova))
+map_dfc(.x = numeric_data_leaves, .f = ~kruskal.test(.x ~ data_leaves$group_anova)$p.value)
+leaves_wilcox_LC <- map(.x = numeric_data_leaves, .f = ~pairwise.wilcox.test(.x, data_leaves_LC$group_anova))
+
+sig_leaves_LC <- leaves_wilcox_LC %>%
+ map(.f = tidy) %>%
+ map2(.y = names(leaves_wilcox), .f = ~.x %>% mutate(var = .y)) %>%
+ purrr::reduce(bind_rows)
+
+data_stems_LC <- chl_calc_LC %>%
+ filter(tissue_group == "stems") %>%
+ select(group_anova, `chlorophyll a`, `chlorophyll b`)
+
+numeric_data_stems_LC <- data_stems_LC %>%
+ select(where(is.numeric))
+
+map_dfc(.x = numeric_data_stems_LC, .f = ~leveneTest(.x ~ data_stems_LC$group_anova))
+map_dfc(.x = numeric_data_stems_LC, .f = ~kruskal.test(.x ~ data_stems_LC$group_anova)$p.value)
+stems_t.test <- map(.x = numeric_data_stems_LC, .f = ~pairwise.t.test(.x, data_stems_LC$group_anova))
+
+sig_stems <- stems_t.test %>%
+ map(.f = tidy) %>%
+ map2(.y = names(stems_t.test), .f = ~.x %>% mutate(var = .y)) %>%
+ purrr::reduce(bind_rows)
+
+#single tests
+#leveneTest(chla_ug_g ~ group_anova, data = data_chl_a_leaves)
+#kruskal.test(chla_ug_g ~ group_anova, data = data_chl_a_leaves)$p.value
+#pairwise.wilcox.test(data_chl_a_leaves$chla_ug_g,
+# data_chl_a_leaves$group_anova)
+
+
+# Scale and Combine both datasets ---------------------------------------------------
+
+
+chl_plot_leaves_LC <- chl_plot_LC %>%
+ filter(pigment != "ratio_chla_chlb",
+ tissue_group == "leaves") %>%
+ mutate(group = group_anova) %>%
+ rename(pigment_label = pigment)
+
+chl_leaves_phot <- chl_plot_leaves %>%
+ filter(!str_detect(pigment_label, "carotenoid")) %>%
+ select(genotype, Tissue, aliquot_ID = `Aliquot ID`, tissue_group, group_anova, pigment_label, concentration) %>%
+ mutate(method = "Photometer")
+
+chl_leaves_LC <- chl_plot_leaves_LC %>%
+ mutate(method = "apolar LC")
+
+id_in_both <- chl_leaves_phot %>%
+ full_join(chl_leaves_LC) %>%
+ distinct(group_anova, aliquot_ID, method) %>%
+ group_by(group_anova, aliquot_ID) %>%
+ summarise(n = n()) %>%
+ ungroup()
+
+chl_leaves_comb <- chl_leaves_phot %>%
+ bind_rows(chl_leaves_LC) %>%
+ left_join(id_in_both) %>%
+ filter(n == 2)
+
+wt_leaves <- chl_leaves_comb %>%
+ filter(group_anova == "gl M82") %>%
+ group_by(pigment_label, method) %>%
+ summarise(wt_mean = mean(concentration))
+
+chl_leaves_scaled <- chl_leaves_comb %>%
+ left_join(wt_leaves) %>%
+ mutate(level = concentration/wt_mean)
+
+chl_leaves_scaled %>%
+ ggplot(aes(x = group_anova, y = level, fill = method)) +
+ geom_boxplot() +
+ geom_point(shape = 21, position = position_dodge(0.75)) +
+ facet_wrap(facets = vars(pigment_label)) +
+ com_theme +
+ theme(legend.position = "bottom")
+
+ggsave("comparison_chlorophyll.jpg", width = 16, height = 20, units = "cm")
+
+chl_leaves_scaled %>%
+ ggplot(aes(x = method, y = level, fill = method)) +
+ geom_boxplot() +
+ geom_point(shape = 21, position = position_dodge(0.75)) +
+ geom_line(aes(x = method, y = level, group = aliquot_ID)) +
+ facet_grid(cols = vars(pigment_label), rows = vars(group_anova), scales = "free")
+
+ggsave("comparison_chlorophyll_ids_outlier.jpg", width = 16, height = 20, units = "cm")
+
+chl_leaves_scaled %>%
+ filter(aliquot_ID != "2131611") %>%
+ ggplot(aes(x = method, y = level, fill = method)) +
+ geom_boxplot() +
+ geom_point(shape = 21, position = position_dodge(0.75)) +
+ geom_line(aes(x = method, y = level, group = aliquot_ID)) +
+ facet_grid(cols = vars(pigment_label), rows = vars(group_anova), scales = "free")
+
+ggsave("comparison_chlorophyll_ids_outlier_removed.jpg", width = 16, height = 20, units = "cm")
diff --git a/workflows/Pigment_SpecPhot_analysis/Pigment_SpecPhot_analysis.R b/workflows/Pigment_SpecPhot_analysis/Pigment_SpecPhot_analysis.R
index 05be75e479ac717e180e647a28e1ecf34359fdd4..04be5ce1993cd25909efee669a2e2c0eb1cfa7aa 100644
--- a/workflows/Pigment_SpecPhot_analysis/Pigment_SpecPhot_analysis.R
+++ b/workflows/Pigment_SpecPhot_analysis/Pigment_SpecPhot_analysis.R
@@ -1,6 +1,3 @@
-rm(list = ls())
-setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
-getwd()
library(tidyverse)
library(car)
library(broom)
@@ -471,7 +468,7 @@ chl_plot_leaves %>%
hide.ns = T) +
facet_grid(cols = vars(pigment_label), scales = "free") +
scale_fill_manual(values = c("darkgreen", "beige")) +
- labs(x = "Tissue and genotype", y = "concentration/ µg/g FW") +
+ labs(x = "Tissue and genotype", y = "concentration/ ?g/g FW") +
com_theme
ggsave(str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"_leaf_pigments.jpg"),
@@ -541,7 +538,7 @@ chl_plot_stems %>%
hide.ns = T) +
facet_grid(cols = vars(pigment_label), scales = "free") +
scale_fill_manual(values = c("darkgreen", "#a0c45a63")) +
- labs(x = "Tissue and genotype", y = "concentration/ µg/g FW") +
+ labs(x = "Tissue and genotype", y = "concentration/ ?g/g FW") +
com_theme
ggsave(str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"_stem_pigments.jpg"),