get credible sets from susieR

[mariesaitou]

Hi, I hope all is well. We are trying to get credible sets from susieR and have several questions.

(0) Basically, what we have is (A) SusieR result, and (B) alpha values from another method based on the (A) SusieR result. And we'd like to get credible sets based on (B), the external alpha values. What would be the most valid way?

(Q1) Susie and susie_get_cs give different results. I think the default parameters are different somewhere. I’m trying to put in something very simple and compare the results to find out the parameters.

For example, the following gives only one credible set,

fitted <-susie(mydata covtest[[1]][["values"]],
      L = 10, verbose = TRUE)
print(fitted$sets)

While the following gives three credible sets.

sets <- susie_get_cs(fitted,
                     coverage = 0.95,
                     min_abs_corr = 0.1)
print(sets)

(Q2) A very naive susie_get_cs(etermal_alpha_matrix) did not work;

> susie_get_cs(alpha, X = NULL,

+             Xcorr = NULL,
+             coverage = 0.95,
+             min_abs_corr = 0.5,
+             dedup = TRUE,
+             squared = FALSE,
+             check_symmetric = TRUE,
+             n_purity = 100
+ )
Error in rep(TRUE, nrow(res$alpha)) : invalid ‘times’ argument 

I appreciate your help.

[stephens999]

for Q1, the default in susie appears to be min_abs_cor=0.5. Maybe that is the difference?

For Q2 I guess you could try putting your alpha values into your fit and run susie_get_cs. ie new_fit = fitted; new_fit$alpha = my_alpha; susie_get_cs(new_fit) ? I should say i'm not 100% sure that will work as I don't recall all the details of susie_get_cs....

[mariesaitou]

Thank you very much, I'll try to investigate it.

[mariesaitou]

Hi, I tried the suggestions above but I got odd results... Could you help me to properly run it?

(1) Is it possible to include the exact same options for susie and susie_get_cs? I ran susie with the following setting, and would like to get cs with the same setting. L = 10, estimate_residual_variance = TRUE , estimate_prior_variance = FALSE, coverage = 0.95, scaled_prior_variance = 0.95, verbose = TRUE I naively put them after susie_get_cs(), but got the following error:

Error in susie_get_cs(new_fit, L = 10, estimate_residual_variance = TRUE,  : 
  unused arguments (L = 10, estimate_residual_variance = TRUE, estimate_prior_variance = FALSE, scaled_prior_variance = 0.95, verbose = TRUE)

(2) Then I ran susie_get_cs(new_fit) with the default setting, but I got oddly many credible sets with almost all the input SNPs. I suspect I have done something wrong.

cs_from_wb <- as.list(NULL)
## get credible set with overwritten alphas with various priors.
for(i in 1:nrow(genelist)){
  alpha_genesubset=subset(qingboYRI,qingboYRI$gene==genelist$gene[i]) # pull from genelist (from 1 to nrow(genelist))
  wb_alpha=alpha_genesubset[,20:29]# position fixed for the all genes  for wb
  new_fit = (fitted.test.YRI[[i]]) 
  new_fit[["alpha"]] = t(as.matrix(wb_alpha))
  cs_from_wb[[i]] <- susie_get_cs(new_fit)
}

And this is the actual R data.

Thank you very much for your help!!

[gaow]

@mariesaitou you can check ?susieR::susie_get_cs for accepted parameters for that function and the defaults. The minimum absolute correlation is 0.5 by default. But in order to perform filters on the CSs you need to provide the LD information through X or Xcorr. For example based on what i figure from the data you provided:

library(susieR)
gene = "ENSG00000272779"
i = which(genelist$gene==gene)
alpha_genesubset=subset(qingboYRI,qingboYRI$gene==genelist$gene[i])  # pull from genelist (from 1 to nrow(genelist))
wb_alpha=alpha_genesubset[,20:29]# position fixed for the all genes  for wb
new_fit = (fitted.test.YRI[[i]]) 
new_fit[["alpha"]] = t(as.matrix(wb_alpha))
# I dont know how you organized the genotype data. This is just a quick guess of what the gneotypes should be for this gene
m = ncol(new_fit[["alpha"]])
p = nrow(genotype.YRI.data)
X = t(genotype.YRI.data[seq(p-m,p),])
tmp_i = susie_get_cs(new_fit, X=X)
tmp_i$cs

Notice that I put in X there which is the key to your problem. The result is null for the gene ENSG00000272779.

[mariesaitou]

Thank you very much, I will look into it!