I wonder if we could use a warnings() (or even a stop()) when people try to feed summary data from MHC region to the program.

The main reason is that LD pattern of MHC is complicated, and may be hard to estimated reliably from 1000 Genomes data.

Related notes: https://crerecombinase.github.io/PolygenicRSS/New_WholeGenome.html#mhc

I wonder if we could add a check to verify that the orders of quh (rotated single-SNP summary statistics) and D (eigenvalues of LD matrix) are the same.

Related notes: https://crerecombinase.github.io/PolygenicRSS/New_WholeGenome.html#bug_fixes

@CreRecombinase Hi Nick,

Below are a few links that accompany my notes posted on Slack — i spent some time after our meeting to clean up these codes; hopefully they are not too hard to read and use:

• fit_rssp.m: fit a RSSp model by maximizing its marginal likelihood
• estimate_pve.m: obtain posterior mean and simulate posterior sample for PVE, given a fitted RSSp

You should have access to these scripts since i have asked Peter to added you to my developer repo. Please let me know if this is not the case.

Before including these scripts into your work, please first make sure they work on your side. To this end I have created a simple example script single_test_nick.m for you to run (you need add input data files there). I only tested these scripts in Matlab 2017b: module load matlab/2017b.

A few more remarks:

• The baseline model that you are mainly interested in corresponds to my RSSp model with only c0 term; that is,

k = 1;
[cvec_1, obj_1] = fit_rssp(vhat, dvec, zeros(k,1));

• fit_rssp.m is built on the fmincon function in Matlab, and so you can use it to confirm your R/C++ implementation of baseline model.

• Can you please apply estimate_pve.m to re-estimate PVE for Framingham eQTL data? It should be not that hard, since all the RSSp models have been fitted and all EVDs have been precomputed.

• Can you please also apply my RSSp models with confounding correction (i.e. with c1, c2, c3, ... terms) to some genes with extremely large PVE estimates in Framingham data? Ideally we want to see reduced PVE estimates for these "problematic" genes.

I wonder if we could let the program compute and display the genomic control factor for the input summary dataset.

Genomic control factor paper: https://www.ncbi.nlm.nih.gov/pubmed/11315092

The following code may be useful.

# SOURCE: http://genometoolbox.blogspot.com/2014/08/how-to-calculate-genomic-inflation.html
# compute genomic control factor
calc_gc <- function(pval) {
  if (all(is.na(pval))) {
    lambda_gc <- NA
    return(lambda_gc)
  } else {
    pval <- pval[!is.na(pval)]
    chisq <- qchisq(1-pval, 1)
    lambda_gc <- median(chisq)/qchisq(0.5,1)
    return(lambda_gc)
  }
}

Displaying this factor can help us understand why sometimes the estimated PVE is large when the true PVE is close to zero (maybe the input summary data were not corrected for population stratification?)