Last updated: 2018-08-31

workflowr checks: (Click a bullet for more information)

  • ✔ R Markdown file: up-to-date

    Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

  • ✔ Repository version: a177f20

    Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.

    Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated: 
    
Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    analysis/.RData
    Ignored:    analysis/.Rhistory
    Ignored:    analysis/ALStruct_cache/
    Ignored:    data/.Rhistory
    Ignored:    docs/figure/.DS_Store

Untracked files:
    Untracked:  .dropbox
    Untracked:  Icon
    Untracked:  _workflowr.yml
    Untracked:  analysis/GTEX-cogaps.Rmd
    Untracked:  analysis/SPCAvRP.rmd
    Untracked:  analysis/eQTL.perm.rand.pdf
    Untracked:  analysis/ieQTL.perm.rand.pdf
    Untracked:  analysis/mash_bhat_z.Rmd
    Untracked:  analysis/mash_ieqtl_permutations.Rmd
    Untracked:  analysis/sc_bimodal.Rmd
    Untracked:  analysis/susie_example.Rmd
    Untracked:  analysis/test_sparse.Rmd
    Untracked:  analysis/z.txt
    Untracked:  code/multivariate_testfuncs.R
    Untracked:  data/4matthew/
    Untracked:  data/4matthew2/
    Untracked:  data/E-MTAB-2805.processed.1/
    Untracked:  data/ENSG00000156738.Sim_Y2.RDS
    Untracked:  data/GDS5363_full.soft.gz
    Untracked:  data/GSE41265_allGenesTPM.txt
    Untracked:  data/Muscle_Skeletal.ACTN3.pm1Mb.RDS
    Untracked:  data/Thyroid.FMO2.pm1Mb.RDS
    Untracked:  data/bmass.HaemgenRBC2016.MAF01.Vs2.MergedDataSources.200kRanSubset.ChrBPMAFMarkerZScores.vs1.txt.gz
    Untracked:  data/bmass.HaemgenRBC2016.Vs2.NewSNPs.ZScores.hclust.vs1.txt
    Untracked:  data/bmass.HaemgenRBC2016.Vs2.PreviousSNPs.ZScores.hclust.vs1.txt
    Untracked:  data/finemap_data/fmo2.sim/b.txt
    Untracked:  data/finemap_data/fmo2.sim/dap_out.txt
    Untracked:  data/finemap_data/fmo2.sim/dap_out2.txt
    Untracked:  data/finemap_data/fmo2.sim/dap_out2_snp.txt
    Untracked:  data/finemap_data/fmo2.sim/dap_out_snp.txt
    Untracked:  data/finemap_data/fmo2.sim/data
    Untracked:  data/finemap_data/fmo2.sim/fmo2.sim.config
    Untracked:  data/finemap_data/fmo2.sim/fmo2.sim.k
    Untracked:  data/finemap_data/fmo2.sim/fmo2.sim.k4.config
    Untracked:  data/finemap_data/fmo2.sim/fmo2.sim.k4.snp
    Untracked:  data/finemap_data/fmo2.sim/fmo2.sim.ld
    Untracked:  data/finemap_data/fmo2.sim/fmo2.sim.snp
    Untracked:  data/finemap_data/fmo2.sim/fmo2.sim.z
    Untracked:  data/finemap_data/fmo2.sim/pos.txt
    Untracked:  data/logm.csv
    Untracked:  data/m.cd.RDS
    Untracked:  data/m.cdu.old.RDS
    Untracked:  data/m.new.cd.RDS
    Untracked:  data/m.old.cd.RDS
    Untracked:  data/mainbib.bib.old
    Untracked:  data/mat.csv
    Untracked:  data/mat.txt
    Untracked:  data/mat_new.csv
    Untracked:  data/paintor_data/
    Untracked:  data/temp.txt
    Untracked:  data/y.txt
    Untracked:  data/y_f.txt
    Untracked:  docs/figure/eigen.Rmd/
    Untracked:  docs/figure/fmo2.sim.Rmd/
    Untracked:  docs/figure/newVB.elbo.Rmd/
    Untracked:  docs/figure/rbc_zscore_mash2.Rmd/
    Untracked:  docs/figure/rbc_zscore_mash2_analysis.Rmd/
    Untracked:  docs/figure/rbc_zscores.Rmd/
    Untracked:  docs/trend_files/
    Untracked:  docs/z.txt
    Untracked:  explore_udi.R
    Untracked:  output/fit.varbvs.RDS
    Untracked:  output/glmnet.fit.RDS
    Untracked:  output/test.bv.txt
    Untracked:  output/test.gamma.txt
    Untracked:  output/test.hyp.txt
    Untracked:  output/test.log.txt
    Untracked:  output/test.param.txt
    Untracked:  output/test2.bv.txt
    Untracked:  output/test2.gamma.txt
    Untracked:  output/test2.hyp.txt
    Untracked:  output/test2.log.txt
    Untracked:  output/test2.param.txt
    Untracked:  output/test3.bv.txt
    Untracked:  output/test3.gamma.txt
    Untracked:  output/test3.hyp.txt
    Untracked:  output/test3.log.txt
    Untracked:  output/test3.param.txt
    Untracked:  output/test4.bv.txt
    Untracked:  output/test4.gamma.txt
    Untracked:  output/test4.hyp.txt
    Untracked:  output/test4.log.txt
    Untracked:  output/test4.param.txt
    Untracked:  output/test5.bv.txt
    Untracked:  output/test5.gamma.txt
    Untracked:  output/test5.hyp.txt
    Untracked:  output/test5.log.txt
    Untracked:  output/test5.param.txt

Unstaged changes:
    Modified:   analysis/_site.yml
    Deleted:    analysis/chunks.R
    Modified:   analysis/eigen.Rmd
    Modified:   analysis/fmo2.sim.Rmd
    Modified:   analysis/newVB.Rmd
    Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
Expand here to see past versions:
    File Version Author Date Message
    Rmd a177f20 stephens999 2018-08-31 workflowr::wflow_publish(“ash-depletion.Rmd”)

Introduction

Let \(J\) denote the number of classes (eg 4 for DNA data).

Suppose we have observed \(J\)-vectors of counts \(x =(x_1,\dots,x_J)\) and \(y=(y_1,\dots,y_J)\), with \[x|p \sim Mult(n_p,p)\] and \[y|q \sim Mult(n_q,q).\]

If \(n_p,n_q\) are large it is natural to use a Poisson approximation: \[x_j \sim Poi(n_p p_j); y_j \sim Poi(n_q q_j)\] from which we have: \[x_j | (x_j+y_j) \sim Bin(x_j+y_j, \rho_j) \quad [*]\] where \[\rho_j = n_p p_j / (n_p p_j + n_q q_j).\]

Now note that \[\log[\rho_j/(1-\rho_j)] = n_p p_j / n_q q_j = \log[n_p/n_q] + \log[p_j/q_j] \quad [**]\] So estimating \(\log(p_j/q_j)\) is effectively the same problem as estimating \(\log(\rho_j/(1-\rho_j))\).

Now a natural esimate of \(\log(\rho_j/(1-\rho_j))\) from [*] is \(\log(x_j/y_j)\), but that does not work when either \(x_j\) or \(y_j\) is 0. We had exactly this problem in smash (Xing and Stephens). In that paper (section B.1) we developed a solution, which gives an estimator for \(\log(\rho_j/(1-\rho_j))\) and its standard error. So the idea is we can use that estimator (subtracting \(\log[n_p/n_q]\) as in [**]) as an estimator of \(\log[p_j/q_j]\). We also have standard errors, and can thus shrink these using ashr (estimating the mode using mode="estimate"). This gives us shrunken estimates of \(\log[p_j/q_j]\), and note that the shrinkage will be strongest for those with large se, which is the ones with small counts (especially 0s!)

This reproducible R Markdown analysis was created with workflowr 1.1.1