Dinh/Dinh 2012/NOTES/2012-2-6: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Dinh
mNo edit summary
>Dinh
No edit summary
Line 51: Line 51:
* The problem line was this:
* The problem line was this:
   $methylStat->add_data($mC/1000); #add in methylation value
   $methylStat->add_data($mC/1000); #add in methylation value
* I changed it to this:
* The correct code should be this:
     $methylStat->add_data(sprintf("%4.3f", $mC/$num_obs)); #add in methylation value
     $methylStat->add_data(sprintf("%4.3f", $mC/$num_obs)); #add in methylation value
* New simulation result:
* New simulation result:
Line 71: Line 71:
  pvalue = 0.1
  pvalue = 0.1


* We should be able to reject fewer CpG sites now.
* '''We should reject fewer CpG sites now.'''

Revision as of 00:23, 9 February 2012

NMF for CpG related to GA

  • 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

  • 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!

  • 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.1
  • We should reject fewer CpG sites now.