WalhoutLab
commented
2 years ago Hi Kamoors,
When the peaks are not clear or the default initial points (expected stats) are very different from where the actual peaks are, the initial points do matter, because the fitting can be stuck in a local minimum and you want to control where the local fitting starts (i.e., around the visible peaks).
In addition, a negative peak does not provide a valid (physically interpretable) set of statistics.
You should therefore change your initial points according to where you think the peaks of subpopulations are. For this distribution, it is hard to judge if you have a right or left shoulder (a shoulder would indicate the second subpopulation under the main peak) in your main population, or even if you have one. It is just clear that you have a blip on the left (seems like a third small subpopulation). So, I would set mu1 (the first peak) at around 2.3 (as far as I can judge where the blip is from the plot) with sigma1 around 0.5, and try to fit the main subpopulations with trial and error. For example, if we believe there is a right shoulder because of a tight subpopulation there, then something like mu2=5+/-1.5, mu3=7+/-0.7 could do (the second numbers are sigma values).
If the negative peak persists in your trials, I would fit two curves instead of three and manually judge the rare expression threshold based on the blip (unless, the blip becomes the secondary subpopulation and the main peak fits the rest well, that also seems possible in this plot and should definitely be tried if you have not already done so).
The order of curves in the output does not matter in the sense that it is not a technical problem; you should just decide on where the high, low and middle peaks are. However, if the order does not follow your initial points, that proves they were not effective in directing the fit.
Finally, if you still cannot come up with a satisfactory fitting, I would be happy to give it a shot if you share the data. This seems like an interesting histogram and may help improve CatExp too.
Safak
-- Lutfu Safak Yilmaz, Ph.D. Associate Professor
Department of Systems Biology University of Massachusetts Chan Medical School 368 Plantation Street, AS5.1057 Worcester MA, 01605-2324 Tel: (508) 856-6711
From: kamoors @.> Sent: Friday, October 14, 2022 11:15 AM To: WalhoutLab/MERGE @.> Cc: Subscribed @.***> Subject: [WalhoutLab/MERGE] CatExp plots (Issue #5)
Hi! me again...
I am now trying to use iMat++ for a completely different model (not iCEL, not Recon), and I'm currently using the CatExp.py stuff to categorize my gene expression data. However, upon using the gaussian curve plotting to obtain the thresholding values, I encounter some weird stuff.. Firstly, I cannot really grasp how to define the expectedStats values, and how exactly they affect the fitting (thus I'm leaving them as in your walkthrough) Second, when I do trimodal, I get a good R^2 value, but the mu values are not in the correct order (why?) - the mu1 matches the mean of the curve most to the right, mu2 matches the one on the left (which is negative?), and mu3 matches the curve in the middle. Is this normal?
Thus, later on I defined the rare cutoff as mu2, low cutoff as mu3, and high cutoff as mu1 - correct? [trimodal]https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fuser-images.githubusercontent.com%2F35234605%2F195881190-19e61072-ef77-4fc6-a9b4-3018a68d54ed.png&data=05%7C01%7CLutfuSSafak.Yilmaz%40umassmed.edu%7Caaa86a6e7e0748c6656408daadf6e8b6%7Cee9155fe2da34378a6c44405faf57b2e%7C0%7C0%7C638013573224323381%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=YqHdN6Gdalf0WWoDs5zyRM8w4sdAHJHwMp1c23%2BozB8%3D&reserved=0 {'mu1': 5.631200169345438, 'sigma1': 1.9850268078063493, 'mu2': 2.585393029325549, 'sigma2': 0.20107748436559067, 'mu3': 4.411776037922199, 'sigma3': 1.0847997002620264, 'N1': 2296.8310628226477, 'N3': 1795.7250793209037, 'N2': -50.55614214355137, 'Rsquared': 0.9960556525960251, 'RSS': 15.82686593199582, 'n': 100}
Thanks in advance :)
— Reply to this email directly, view it on GitHubhttps://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2FWalhoutLab%2FMERGE%2Fissues%2F5&data=05%7C01%7CLutfuSSafak.Yilmaz%40umassmed.edu%7Caaa86a6e7e0748c6656408daadf6e8b6%7Cee9155fe2da34378a6c44405faf57b2e%7C0%7C0%7C638013573224323381%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=cc3yawTBlj%2FAW9JmYGzxTMyJ5k1hmWcAV%2Fb2MaqSxCo%3D&reserved=0, or unsubscribehttps://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAA3DYFIQBUPGDIPVMXWVV2DWDF2IPANCNFSM6AAAAAARFKV4VI&data=05%7C01%7CLutfuSSafak.Yilmaz%40umassmed.edu%7Caaa86a6e7e0748c6656408daadf6e8b6%7Cee9155fe2da34378a6c44405faf57b2e%7C0%7C0%7C638013573224479657%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=jMCpevZRDgmIYM5eN62KQ0HD1RNC5%2Bx9chm0GZFSLss%3D&reserved=0. You are receiving this because you are subscribed to this thread.Message ID: @.***>
kamoors
Yet another interesting case.. I cannot seem to get a good fit (probably because of the plateau) 
Hi! me again...
I am now trying to use iMat++ for a completely different model (not iCEL, not Recon), and I'm currently using the CatExp.py stuff to categorize my gene expression data. However, upon using the gaussian curve plotting to obtain the thresholding values, I encounter some weird stuff.. Firstly, I cannot really grasp how to define the expectedStats values, and how exactly they affect the fitting (thus I'm leaving them as in your walkthrough) Second, when I do trimodal, I get a good R^2 value, but the mu values are not in the correct order (why?) - the mu1 matches the mean of the curve most to the right, mu2 matches the one on the left (which is negative?), and mu3 matches the curve in the middle. Is this normal?
Thus, later on I defined the rare cutoff as mu2, low cutoff as mu3, and high cutoff as mu1 - correct?
{'mu1': 5.631200169345438,
'sigma1': 1.9850268078063493,
'mu2': 2.585393029325549,
'sigma2': 0.20107748436559067,
'mu3': 4.411776037922199,
'sigma3': 1.0847997002620264,
'N1': 2296.8310628226477,
'N3': 1795.7250793209037,
'N2': -50.55614214355137,
'Rsquared': 0.9960556525960251,
'RSS': 15.82686593199582,
'n': 100}
Thanks in advance :)