One SNP per credible set, is it normal?

[Ruiqi-CUB]

Hi there! Thank you so much for developing this awesome software!

I am new to fine mapping and I find it is very useful! However, I was wondering if you could help me with you experience to see if the results I got is biologically meaningful, or there is something going wrong.

I followed the tutorial here to do fine mapping on my dataset, after running susieR, I got 6 credible sets, and each of them contains one SNP as shown on the figure:

I was just wondering if this is "normal". Would you mind helping me have a look?

Here is a breif describtion of my data and my code:
Input for susieR:

1. LD matrix: SNPs in the +-1MB of the target set was extracted, and LD matrix were calculated by the plink command below:

plink --bfile OZ034921.1_20700459_1MB_GWAS_region \
    --keep-allele-order \
    --r square \
    --extract OZ034921.1_20700459_1MB_GWAS_snps.txt \
    --out OZ034921.1_20700459_1MB_GWAS_region_LD \
    --allow-extra-chr \
    --noweb

2. gwas results: GWAS was run on a lightly filtered SNP set, and then results were extracted for SNPs in that +-1Mb window using the same file as above OZ034921.1_20700459_1MB_GWAS_snps.txt

I made sure that the order of snps match in those two files.

SusieR was run with the following command:

# read gwas
gwas <- fread("OZ034921.1_20700459_1MB_GWAS_region_GWAS.txt", header = FALSE)
setnames(gwas, c(
  "chr", "rs", "ps", "n_miss", "allele1", "allele0", "af",
  "beta", "se", "logl_H1", "l_remle", "l_mle",
  "p_wald", "p_lrt", "p_score"
))

# Read LD matrix without assigning column names
ld_raw <- fread("OZ034921.1_20700459_1MB_GWAS_region_LD.ld", header = FALSE)
# Convert to numeric matrix and remove dimnames
ld.mat <- as.matrix(ld_raw)
storage.mode(ld.mat) <- "numeric"
dimnames(ld.mat) <- NULL

fitted_rss <- susie_rss(
  bhat = gwas$beta,
  shat = gwas$se,
  R = ld.mat,
  n = 151,
  max_iter = 500  # increase iteration cap
)

The only warning I saw is:
"Failed to query server: Transport endpoint is not connected
HINT: For large R or large XtX, consider installing the Rfast package for better performance.
Warning message:
In system("timedatectl", intern = TRUE) :
running command 'timedatectl' had status 1" but it seems unrelevant

Here is the PIP plot if it is helpful

Thank you very much!
Ruiqi

[pcarbo]

@Ruiqi-CUB I see that the sample size is n = 151, which is quite small. Do you have access to the full (genotype and phenotype) data? If so, you could try running susie() on the full data, and susie_rss() on the summary data, and compare the results.

[Ruiqi-CUB]

Thank you Peter! I work on non-model organisms so 151 is the best that I can get.
I am not too familiar with susie() and the only tutorial I found is here. I was wondering how I should convert my gwas data to use susie(). Would you mind checking if my understanding is correct?

1. X: I do have the genotype data in vcf format. How should I convert it to the matrix that I can use for susie()?
2. Y: the trait I used in GWAS is a discrete trait (0 or 1) (could this also be the issue?). My understanding is that X should be a vector with contains those trait value. Is this correct?

[pcarbo]

@Ruiqi-CUB Answering your questions:

1. In the past I have used PLINK to convert a VCF file to a .raw file, then this can be easily read into R using, say, fread from the data.table package. There are surely better ways to do this using Bioconductor packages, but this is one approach that I am familiar with. What you want in the end is a numeric matrix in which the genotypes are encoded as allele counts (0, 1 or 2).

2. One can reasonably approximate a logistic regression with a linear regression as long as the individual SNP effects are not "too large". This is discussed for example in Pirinen et al 2013.

Hope this helps.

[Ruiqi-CUB]

Thank you! Just did it with susie():
X: genotype in the 1Mb window converted by plink

str(X)
 num [1:151, 1:3024] 0 0 0 0 1 0 0 0 0 0 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:3024] "OZ034921.1_19701047_A_C_C" "OZ034921.1_19701405_G_A_A" "OZ034921.1_19702877_A_T_A" "OZ034921.1_19703085_C_G_G" ...

Y: binary trait

str(Y)
 int [1:150] 1 0 1 1 1 1 0 0 1 1 ...

I did susier:

fitted <- susie(X, Y,
                L = 10,
                estimate_prior_variance = TRUE,
                estimate_residual_variance = TRUE,
                verbose = TRUE)

but the results look terrible:

# Number of credible sets
length(cs$cs)
[1] 10
> sapply(cs$cs, length)
  L1   L2   L3   L4   L5   L6   L7   L8   L9  L10
 794 2849 2839 2828 2820 2818 2822 2829 2838 2846
> summary(fitted$pip)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
0.002357 0.002420 0.002618 0.003300 0.003164 0.136920

Would you mind giving me some suggestions? Should I try other bayesian methods? or should I just use heuristic fine mapping, or even just finding genes in proximity?

[pcarbo]

@Ruiqi-CUB Before running SuSiE, did you try running an association analysis (e.g., with PLINK) to determine if any of the SNPs have strong (and perhaps significant) associations with the binary phenotype?

[Ruiqi-CUB]

No I did not.

I was planning to go with GWAS route with susie_rss(). I did find SNPs with pvalues above the cutoff from GWAS.
Should I run an association analysis to see if the results overlapping with significant SNPs from GWAS?

I am pretty new to find mapping. Any further suggestions would be much appreciate! Thank you!

[pcarbo]

@Ruiqi-CUB Typically the process is to first run an association analysis to find regions of the genome that are significantly associated with the phenotype ("QTLs"), then the second step is to finemap the regions containing these significant associations.

[Ruiqi-CUB]

Oh I think used GWAS to indentify QTLs, and I am fine mapping 1MB window around SNPs that have p values below 5 x 10-8. Do you think that I should try other methods for identifying QTLs and fine mapping?

[pcarbo]

@Ruiqi-CUB Okay, then I would start by fine-mapping some QTLs that contain the strongest associations. It is hard to say why your credible sets contain so many SNPs, but one possible explanation is that there is very strong and/or long-range LD (correlation) between many SNPs. Hope this helps. With some luck, some of the QTLs you have identified have less strong LD and therefore offer the potential to better finemap the SNPs.