Dinh/Dinh 2012/NOTES/2012-2-6
Jump to navigation
Jump to search
[edit]
- Can we use NMF to cluster samples?
- A = methylation matrix (nxm),
- For k>1, find W (nxk) and H (kxm) such that A ~ WH
- Samples are clustered based on W and H
- The NMF problem is solved without requiring orthogonality (as in PCA) - such that there are many possible W and H.
- To solve NMF, the algorithm initializes W and H randomly, then reiterate toward an objective function (minimizing the distance between WH and A)
- Because the algorithm doesn't always converge to the same answer, we can calculate the robustness of the clustering by generating a consensus matrix from multiple runs.
- Install packages "NMF" and "RColorBrewer" on meangenemachine.
- Downloaded matrix for 35 sites from Dr. Zhang's wiki: http://genome-tech.ucsd.edu/LabNotes/index.php/Kun:LabNotes/CpgSeq/2012-2-6
Simulation p-values for CpG sites in 9 technical replicates[edit]
- I only finished the simulation for ~184K sites, but the matrix was cut-off at the bottom, there was ~400K sites to be tested.
- Out of the ~184K sites, 64K have STD >= 0.1
- Within the 64K that have STD >= 0.1, the distribution of p-values are as follows:
Mean = 0.00897, SD = 0.0922 Number of sites p < 0.05 = 63,667 (99%) Number of sites p < 0.01 = 63,622 (99%) Number of sites p < 0.001 = 63,578 (99%)
- Nearly all sites with STD >= 0.1 across the technical replicates are bad.
- ~120K sites have STD < 0.1 across the technical replicates
- Within the ~120K sites with STD < 0.1, the distribution of p-values are as follows:
Mean = 0.35, SD = 0.47 Number of sites p < 0.05 = 77110 (64%) Number of sites p < 0.01 = 76406 (63%) Number of sites p < 0.001 = 75016 (62%) However, a majority of these sites have very low standard deviation 50% have STD <= 0.02854 75% have STD <= 0.05926
- If we only reject CpGs with p<0.05 and STD>=0.1, we will reject ~36% of sites (64K)
Update simulation script![edit]
- According to Dr. Zhang's suggestion, I looked at the simulated values at 1 CpG site and saw that it was calculating the methylation value incorrectly.
Observed STD = 0.0584 Est_mC = 0.968 Round# mC#1 mC#2 mC#3 mC#4 mC#5 mC#6 mC#7 mC#8 mC#9 Simulated STD RD 14 12 69 83 68 58 107 91 85 0 0.005 0.003 0.016 0.021 0.066 0.058 0.033 0.024 0.027 0.0216 1 0.004 0.003 0.016 0.023 0.066 0.057 0.033 0.025 0.027 0.0215 2 0.005 0.003 0.017 0.023 0.064 0.058 0.032 0.024 0.026 0.0210 3 0.005 0.003 0.017 0.023 0.066 0.055 0.033 0.023 0.027 0.0210 4 0.005 0.003 0.017 0.023 0.064 0.055 0.033 0.024 0.026 0.0205 5 0.005 0.003 0.015 0.023 0.064 0.057 0.033 0.024 0.027 0.0210 6 0.004 0.003 0.016 0.021 0.067 0.057 0.033 0.023 0.027 0.0218 7 0.005 0.003 0.016 0.021 0.068 0.058 0.034 0.024 0.026 0.0221 8 0.004 0.003 0.017 0.022 0.067 0.057 0.031 0.025 0.027 0.0216 9 0.005 0.003 0.017 0.023 0.066 0.056 0.032 0.024 0.023 0.0212 10 0.005 0.003 0.017 0.023 0.066 0.056 0.034 0.023 0.025 0.0212 pvalue = 0.0
- The problem line was this:
$methylStat->add_data($mC/1000); #add in methylation value
- The correct code should be this:
$methylStat->add_data(sprintf("%4.3f", $mC/$num_obs)); #add in methylation value
- New simulation result:
Observed STD = 0.0584 Est_mC = 968 Round# mC#1 mC#2 mC#3 mC#4 mC#5 mC#6 mC#7 mC#8 mC#9 Simulated STD RD 14 12 69 83 68 58 107 91 85 0 0.800 1.000 1.000 0.913 0.985 0.948 0.912 0.960 1.000 0.0652 1 1.000 1.000 0.941 1.000 0.971 0.966 1.000 1.000 0.926 0.0289 2 0.800 1.000 0.941 0.957 1.000 0.983 0.971 1.000 0.963 0.0626 3 1.000 1.000 1.000 0.913 0.971 0.966 1.000 1.000 0.963 0.0297 4 1.000 1.000 1.000 1.000 0.956 0.966 1.000 0.960 1.000 0.0198 5 1.000 1.000 0.941 0.957 0.956 0.948 0.971 1.000 0.926 0.0279 6 1.000 1.000 1.000 1.000 0.956 0.931 0.971 1.000 0.889 0.0399 7 1.000 1.000 0.941 0.913 0.956 0.983 0.971 1.000 1.000 0.0315 8 1.000 1.000 1.000 0.957 0.985 1.000 0.941 0.960 0.963 0.0233 9 1.000 1.000 0.941 0.957 0.941 0.966 0.971 0.960 1.000 0.0241 10 1.000 1.000 0.882 1.000 0.956 0.983 0.971 1.000 1.000 0.0390 pvalue = 0.2
- We should reject fewer CpG sites now.