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

From ZhangLabWiki
Jump to navigation Jump to search
>Dinh
No edit summary
>Dinh
mNo edit summary
 
Line 69: Line 69:
  9      1.000  1.000  0.941  0.957  0.941  0.966  0.971  0.960  1.000  0.0241
  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
  10      1.000  1.000  0.882  1.000  0.956  0.983  0.971  1.000  1.000  0.0390
  pvalue = 0.1
  pvalue = 0.2


* '''We should reject fewer CpG sites now.'''
* '''We should reject fewer CpG sites now.'''

Latest revision as of 00:24, 9 February 2012

NMF for CpG related to GA[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.