Run isoTWAS associations
stageassoc.RmdThis module will cover how to use functions from isotwas
to map isoform-level trait associations using a two-step testing
framework. Then, we will run probabilistic fine-mapping to account of
linkage disequilibrium (LD) in a locus with overlapping isoforms that
are associated with a trait. Along the way, we’ll be covering some
practical considerations for computation and feature selection.
For this, you will need:
a set of isoTWAS models,
an appropriate LD reference panel, and
GWAS summary statistics for your trait of interest.
Let’s talk a little about each of these three components.
Pre-computed models and LD matrix
You can always use your own set of models that you have trained on your own isoform expression QTL data (paired genotype and RNA-seq data). For convenience, we have trained and made publicly available isoTWAS models for 48 tissues from the Genotype-Tissue Expression Project (GTEx), adult brain cortex from PsychENCODE and AMP-AD, and developmental brain cortex from PsychENCODE. There are Zenodo repositories for the GTEx, adult cortex, and developmental cortex models. Please cite our manuscript and the associated Zenodo DOI if you use these models.
To download an entire Zenodo repository, you may use the Python
package Zenodo_get (https://zenodo.org/record/1261813). To install and use,
here are some sample commands to directly download all of the adult
cortex repository off Zenodo:
In this tutorial, we are only using isoTWAS models for 3 genes:
CDC42 (ENSG00000070831), WNT4
(ENSG00000162552), and CDC42-IT1
(ENSG00000230068). Let’s take a look at what the models look
like, which are saved as .RDS files.
set.seed(1789)
model_test = readRDS(system.file("extdata",
"ENSG00000070831_isoTWAS.RDS",
package = "isotwas"))
head(model_test)## Feature SNP Chromosome Position Build ALT REF Weight
## 2 ENST00000344548 rs10917023 1 21599110 hg38 A G 0.0257706583
## 3 ENST00000344548 rs10917145 1 22086152 hg38 A G -0.2503237794
## 4 ENST00000344548 rs11584687 1 22550747 hg38 A G 0.0200388736
## 5 ENST00000344548 rs12082914 1 21493245 hg38 G T 0.0167839153
## 6 ENST00000344548 rs12118362 1 21445504 hg38 G A 0.0133152139
## 7 ENST00000344548 rs12127343 1 21543280 hg38 A G 0.0007759457
## R2 Gene
## 2 0.09314035 ENSG00000070831
## 3 0.09314035 ENSG00000070831
## 4 0.09314035 ENSG00000070831
## 5 0.09314035 ENSG00000070831
## 6 0.09314035 ENSG00000070831
## 7 0.09314035 ENSG00000070831If you train your own models, ensure that your model files contain columns for the:
isoform name,
SNP name or rsID,
chromosome of the SNP,
position of the SNP,
alternative allele,
reference allele, and
the predictive weight for the isoform.
For reference, our models include the genome build on which the model was trained, so users can align their GWAS to the models based on genome build.
We also include in-sample LD matrices for each isoTWAS model. These
are simply .RDS files that contain a square matrix with the
LD between the SNPs in the correspondin row and column.
LD_test = readRDS(system.file("extdata",
"ENSG00000070831_LDMatrix.RDS",
package = "isotwas"))
head(as.matrix(LD_test)[,1:10])## Loading required package: Matrix## rs10916927 rs6658526 rs1354792 rs12567861 rs10916930 rs7521711
## rs10916927 1.0000000 0.9440432 0.8544395 -0.6393393 -0.6405133 0.8536513
## rs6658526 0.9440432 1.0000000 0.8005756 -0.5896286 -0.5924963 0.7998371
## rs1354792 0.8544395 0.8005756 1.0000000 -0.7649469 -0.7657869 0.9990775
## rs12567861 -0.6393393 -0.5896286 -0.7649469 1.0000000 0.9971883 -0.7642413
## rs10916930 -0.6405133 -0.5924963 -0.7657869 0.9971883 1.0000000 -0.7650804
## rs7521711 0.8536513 0.7998371 0.9990775 -0.7642413 -0.7650804 1.0000000
## rs4654899 rs960564 rs6663452 rs7532667
## rs10916927 0.9097877 0.9966544 0.9501309 0.9971286
## rs6658526 0.8503948 0.9427488 0.8973393 0.9458899
## rs1354792 0.9361329 0.8530981 0.8951509 0.8538566
## rs12567861 -0.7173327 -0.6390579 -0.6724022 -0.6393520
## rs10916930 -0.7179086 -0.6402319 -0.6737519 -0.6405219
## rs7521711 0.9371127 0.8532494 0.8952461 0.8530689Note that the column and row names for the matrix correspond to the rsIDs in the model.
Practical consideration: In many settings,
you will not have a corresponding file that contains the LD matrix for
the SNPs at the locus. In this case, we suggest either one of two
approaches. First, you can generate the LD matrix from the data where
the model was trained, if you have access to the individual-level
genotypes for this data. Alternatively, you can use an LD reference
panel that is well-aligned to the genetic ancestry of the training
dataset. This LD reference panel is usually 1000 Genomes, since it’s
publicly available. You can use bigsnpr::snp_cor or
bigsnpr::bed_cor to generate similar-looking LD matrices,
given you have the PLINK style files
(.bed/.bim/.fam) for the LD
reference.
Preparing the GWAS summary statistics
The last file you will need is the GWAS summary statistics for your trait of interest. In this tutorial, we are using the GWAS summary statistics for schizophrenia, made available by the Psychiatric Genomics Consortium. We are specifically using the GWAS summary statistics from European-ancestry individual. You can find these files here: https://figshare.com/articles/dataset/scz2022/19426775.
Since GWAS summary statistics often have varied column names, we use a munging script to clean the data. In addition, we subset the GWAS summary statistics to only SNPs annotated in HapMap3 to ensure that there is maximal overlap. This munging script is available here, courtesy of Michael Margolis.
In this tutorial, we have subsetted only to the SNPs that are within the locus of the example genes. Let’s take a look at this munged file.
gwas = read.table(system.file("extdata",
"tutorial_SCZ_2022.hm3_filter.sumstats.gz",
package = "isotwas"),
header=T)
head(gwas)## SNP A1 A2 Z N
## 1 rs10916927 T C 3.108 146347.8
## 2 rs6658526 C T 2.954 146347.8
## 3 rs1354792 T C 1.636 146347.8
## 4 rs12567861 C T -1.641 146347.8
## 5 rs10916930 T C -1.684 146347.8
## 6 rs7521711 A G 1.627 146347.8This munged summary statistics table has only 5 columns:
The SNP rsID,
the alternative allele,
the reference allele,
the Z-score (beta value divided by standard error) for the GWAS association, and
the effective sample size.
We advise you use the munging script to format your summary statistics files. At the very least, please ensure that your summary statistics file contains a column for the SNP identifier, the alternative allele, and either the Z-score or the beta value and standard error for each SNP.
Generating nominal Z-scores
The first step for isoTWAS trait-mapping is to generate the isoform-level Z-scores using the weighted burden test. In general, you’ll have a vector of gene names with corresponding isoTWAS model files and/or LD files. The following code snippet loops through this vector of genes, reads in the model file and LD file, and then computes a Z-score for each isoform of the gene that was well-predicted in cross-validation.
gene_names = c('ENSG00000070831',
'ENSG00000162552',
'ENSG00000230068')
out_df = data.frame(Gene = c(),
Feature = c(),
Z = c(),
P = c(),
permute.P = c(),
topSNP = c(),
topSNP.P = c())
for (gene in gene_names){
model = readRDS(system.file("extdata",
paste0(gene,"_isoTWAS.RDS"),
package = "isotwas"))
ld = readRDS(system.file("extdata",
paste0(gene,"_LDMatrix.RDS"),
package = "isotwas"))
for (tx in unique(model$Feature)){
sumstats.cur = subset(gwas,SNP %in% subset(model, Feature == tx)$SNP)
tx_df = isotwas::burdenTest(mod = subset(model, Feature == tx),
ld = ld,
gene = gene,
sumStats = sumstats.cur,
chr = 'Chromosome',
pos = 'Position',
a1 = 'A1',
a2 = 'A2',
a1_mod = 'ALT',
a2_mod = 'REF',
snpName = 'SNP',
Z = 'Z',
beta = NULL,
se = NULL,
featureName = 'Feature',
R2cutoff = .01,
alpha = 1e-3,
nperms = 1000,
usePos = F)
out_df = rbind(out_df,tx_df)
}
}## Registered S3 methods overwritten by 'huge':
## method from
## plot.roc pROC
## plot.sim lava
## print.roc pROC
## print.sim lava## Registered S3 methods overwritten by 'spls':
## method from
## predict.splsda caret
## print.splsda caret
head(out_df)## Gene Feature Z P permute.P topSNP
## 1 ENSG00000070831 ENST00000344548 -3.239210 0.0011986108 1.00000000 rs2143103
## 2 ENSG00000070831 ENST00000315554 1.268410 0.2046514827 1.00000000 rs11582542
## 3 ENSG00000070831 ENST00000656825 3.268385 0.0010816323 1.00000000 rs11582542
## 4 ENSG00000070831 ENST00000498236 -2.906460 0.0036554386 1.00000000 rs11582542
## 5 ENSG00000070831 ENST00000400259 3.344840 0.0008232993 0.01198801 rs11582542
## 6 ENSG00000070831 ENST00000667384 0.575199 0.5651567515 1.00000000 rs11582542
## topSNP.P
## 1 0.0004331767
## 2 0.0003396706
## 3 0.0003396706
## 4 0.0003396706
## 5 0.0003396706
## 6 0.0003396706The out_df variable contains a data frame with 7
columns:
the gene name,
the isoform name,
the isoform-level Z-score for the trait association,
the P-value corresponding to this Z-score,
the permutation P-value (more on this below),
the top SNP based on GWAS P-value within 1 Mb of the gene locus,
and the P-value corresponding to this GWAS SNP.
The permutation test test assesses how much signal is added by
isoform expression, given the GWAS architecture of the locus, and
controls for large LD blocks. We shuffle the SNP-to-isoform weights to
generate a null distribution and use this null to generate the
permutation P-value. We only run this permutation test if the nominal
P-value from the isoform-level Z-score is less than alpha,
as defined in the burdenTest function.
Stage-wise hypothesis test
isoTWAS employs a two-step testing framework. First, we
control from false discovery rate (FDR) across gene families, using
either Benjamini-Hochberg FDR control (recommended) or Bonferroni. Next,
for genes that pass FDR control (adjusted P-values for the genes are
lower than our defined significance threshold), we drop down to the
isoform level and run family-wise error rate (FWER) control using
Shaffer’s modified sequentially rejective Bonferroni procedure. Here, we
are controlling both the FDR (alpha1) and FWER
(alpha2) to 0.05.
Let’s run the first stage now: gene-level FDR control.
suppressPackageStartupMessages(library(tidyverse))
gene = out_df %>%
group_by(Gene) %>%
summarise(Screen.P = isotwas::p_screen(P))
gene = as.data.frame(gene)
head(gene)## Gene Screen.P
## 1 ENSG00000070831 0.0016636386
## 2 ENSG00000162552 0.0008302377
## 3 ENSG00000230068 0.0004539979
alpha1=.05
G = nrow(gene)
gene$Screen.P.Adjusted = p.adjust(gene$Screen.P,method = 'fdr')
R = length(unique(gene$Gene[gene$Screen.P.Adjusted < alpha1]))
alpha2 = (R*alpha1)/G
print(gene)## Gene Screen.P Screen.P.Adjusted
## 1 ENSG00000070831 0.0016636386 0.001663639
## 2 ENSG00000162552 0.0008302377 0.001245357
## 3 ENSG00000230068 0.0004539979 0.001245357We see that, on the gene-level, all 3 genes pass FDR control. That is, we see that the adjusted screening P-value is less than 0.05 for all 3 genes. Now, we have to drop down to the isoform level for all 3 genes.
Let’s run the second stage: isoform-level FWER control for these gene families.
isoform_new = as.data.frame(matrix(nrow = 0,
ncol = ncol(out_df)+2))
colnames(isoform_new) = c(colnames(out_df),'Screen.P','Confirmation.P')
gene = gene[order(gene$Screen.P),]
ttt = merge(out_df,
gene[,c('Gene','Screen.P',
'Screen.P.Adjusted')],
by = 'Gene')
isoform_new = ttt %>%
group_by(Gene) %>%
summarise(Feature = Feature,
Confirmation.P = isotwas::p_confirm(P,alpha = alpha2))## `summarise()` has grouped output by 'Gene'. You can override using the
## `.groups` argument.
isoform_new = merge(isoform_new,ttt,by=c('Gene','Feature'))
isoform_new$Confirmation.P = ifelse(isoform_new$Screen.P.Adjusted < 0.05,
isoform_new$Confirmation.P,
1)
isoform_new = isoform_new[,c('Gene','Feature','Z','P','permute.P',
'topSNP','topSNP.P',
'Screen.P','Screen.P.Adjusted','Confirmation.P')]
print(isoform_new)## Gene Feature Z P permute.P
## 1 ENSG00000070831 ENST00000315554 1.268410 0.2046514827 1.000000000
## 2 ENSG00000070831 ENST00000344548 -3.239210 0.0011986108 1.000000000
## 3 ENSG00000070831 ENST00000400259 3.344840 0.0008232993 0.011988012
## 4 ENSG00000070831 ENST00000498236 -2.906460 0.0036554386 1.000000000
## 5 ENSG00000070831 ENST00000656825 3.268385 0.0010816323 1.000000000
## 6 ENSG00000070831 ENST00000662562 -3.278162 0.0010448529 1.000000000
## 7 ENSG00000070831 ENST00000667384 0.575199 0.5651567515 1.000000000
## 8 ENSG00000162552 ENST00000290167 -3.441877 0.0005776926 0.041958042
## 9 ENSG00000162552 ENST00000415567 -3.179540 0.0014750883 1.000000000
## 10 ENSG00000230068 ENST00000431803 -3.506525 0.0004539979 0.007992008
## topSNP topSNP.P Screen.P Screen.P.Adjusted Confirmation.P
## 1 rs11582542 0.0003396706 0.0016636386 0.001663639 1.0000000000
## 2 rs2143103 0.0004331767 0.0016636386 0.001663639 0.0062691173
## 3 rs11582542 0.0003396706 0.0016636386 0.001663639 0.0049397958
## 4 rs11582542 0.0003396706 0.0016636386 0.001663639 0.0109663157
## 5 rs11582542 0.0003396706 0.0016636386 0.001663639 0.0062691173
## 6 rs2143103 0.0004331767 0.0016636386 0.001663639 0.0062691173
## 7 rs11582542 0.0003396706 0.0016636386 0.001663639 1.0000000000
## 8 rs11582542 0.0003396706 0.0008302377 0.001245357 0.0005776926
## 9 rs11582542 0.0003396706 0.0008302377 0.001245357 0.0014750883
## 10 rs11582542 0.0003396706 0.0004539979 0.001245357 0.0000000000## Gene Feature Z P permute.P
## 3 ENSG00000070831 ENST00000400259 3.344840 0.0008232993 0.011988012
## 8 ENSG00000162552 ENST00000290167 -3.441877 0.0005776926 0.041958042
## 10 ENSG00000230068 ENST00000431803 -3.506525 0.0004539979 0.007992008
## topSNP topSNP.P Screen.P Screen.P.Adjusted Confirmation.P
## 3 rs11582542 0.0003396706 0.0016636386 0.001663639 0.0049397958
## 8 rs11582542 0.0003396706 0.0008302377 0.001245357 0.0005776926
## 10 rs11582542 0.0003396706 0.0004539979 0.001245357 0.0000000000We have now generated screening gene-level FDR-adjusted P-values and confirmation isoform-level FWER controlled P-values. We see that there are 3 isoforms at this locus that pass all 3 levels of hypothesis testing: screening FDR-adjusted P < 0.05, confirmation FWER-adjusted P < 0.05, and nominal permutation P < 0.05. These isoforms are ENST00000400259 (gene: ENSG00000070831), ENST00000290167 (gene: ENSG00000162552), and ENST00000431803 (gene: ENSG00000230068).
Fine-mapping overlapping and trait-associated isoforms
Since these two isoforms overlap and likely share LD-linked SNPs that predict each isoform, we can’t be sure which of the three isoforms are likely to be carrying the trait association signal. In short, we account of the correlations between the genetically-predicted portions of the expression of these isoforms, enumerate the possible causal configurations, and use a Bayesian approach to generate a posterior inclusion probability (PIP) for each isoform. These PIPs can be used to generate a 90% credible set of isoforms that best explain the trait association at the locus. See Mancuso et al’s manuscript about gene-level fine-mapping for more of the intuition and mathematical details.
For fine-mapping, you no longer need the GWAS summary statistics. You will need:
TWAS summary statistics (something that looks like
isoform_new),the isoTWAS models, and
an LD reference panel.
Here, you need to make sure the LD reference matrix contains all SNPs that are in all models for all isoforms in the overlapping locus. If you do not have the in-sample individual-level genotypes, we recommend using 1000 Genomes. For ease in this tutorial, we have compiled the necessary LD matrix.
First, we annotate each isoform with location and select the significantly associated isoforms that are within 1 Mb of one another. We can obtain the gene annotation with the following code:
suppressPackageStartupMessages(require(biomaRt))
gene_names = unique(isoform_new$Gene)
ensembl <- useEnsembl(biomart = "ensembl",
dataset = "hsapiens_gene_ensembl",
mirror = "useast")
bm = getBM(attributes = c('ensembl_gene_id',
'chromosome_name',
'start_position',
'end_position'),
filters = 'ensembl_gene_id',
values = gene_names,
mart = ensembl)
colnames(bm) = c('Gene','Chromosome','Start','End')Next, we select the overlapping isoforms.
bm = data.frame(Gene = c('ENSG00000070831','ENSG00000162552','ENSG00000230068'),
Chromosome = 1,
Start = c(22052627,22117313,22059197),
End = c(22101360,22143969,22064199))
isoform_new = merge(bm,isoform_new,by='Gene')
isoform_new = isoform_new[order(isoform_new$Chromosome,
isoform_new$Start,
decreasing = F),]
isoform_sig = subset(isoform_new,
Screen.P.Adjusted < alpha1 &
Confirmation.P < alpha2 &
permute.P < 0.05)
keep.isoform = c()
if (nrow(isoform_sig) > 1){
for (i in 1:(nrow(isoform_sig)-1)){
if (isoform_sig$End[i] > isoform_sig$Start[i+1] - 1e6){
keep.isoform = unique(c(keep.isoform,
c(isoform_sig$Feature[c(i,i+1)])))
}
}
}
isoform_sig = subset(isoform_sig,Feature %in% keep.isoform)Significant isoforms that do not overlap within 1 Mb of another significant isoform do not need to be fine-mapped.
Now, we aggregate the models and generate a single table with the SNP-to-isoform weights of all the SNPs that predict these overlapping isoforms.
all.snps = c()
omega = c()
pos = c()
gene = c()
snp.chr = c()
### COLLECT WEIGHTS FOR SNPS IN THE MODELS
for (i in 1:nrow(isoform_sig)){
gene_in = isoform_sig$Gene[i]
model_in = readRDS(system.file("extdata",
paste0(gene_in,"_isoTWAS.RDS"),
package = "isotwas"))
model_in = subset(model_in,
Feature == isoform_sig$Feature[i])
Model = data.frame(SNP = model_in$SNP,
Chromosome = model_in$Chromosome,
Position = model_in$Position,
Effect = model_in$Weight,
A1 = model_in$ALT,
A2 = model_in$REF)
Model = subset(Model,Effect!=0)
Model = Model[!duplicated(Model$SNP),]
all.snps = c(all.snps,
as.character(Model$SNP))
omega = c(omega,
as.numeric(Model$Effect))
gene = c(gene,
rep(isoform_sig$Feature[i],nrow(Model)))
snp.chr = c(snp.chr,
as.numeric(Model$Chromosome))
pos = c(pos,as.numeric(Model$Position))
}
tot.df = data.frame(SNP = all.snps,
Gene = gene,
Effect = omega,
Chromosome = snp.chr)
model.df = as.data.frame(matrix(nrow = length(unique(all.snps)),
ncol = nrow(isoform_sig)+1))
colnames(model.df) = c('SNP',isoform_sig$Feature)
model.df$SNP = as.character(unique(all.snps))
for (q in 1:nrow(isoform_sig)){
cur.tot.df = subset(tot.df,Gene == isoform_sig$Feature[q])
cur.tot.df$SNP = as.character(cur.tot.df$SNP)
for (i in 1:nrow(model.df)){
w = which(cur.tot.df$SNP == model.df$SNP[i])
model.df[i,q+1] = ifelse(length(w) != 0,
cur.tot.df$Effect[w],
0)
}
}
model.df$Chromosome = 2
for (i in 1:nrow(model.df)){
rrr = subset(tot.df,SNP == model.df$SNP[i])
model.df$Chromosome[i] = rrr$Chromosome[1]
}
head(model.df)## SNP ENST00000400259 ENST00000431803 ENST00000290167 Chromosome
## 1 rs1002480 -3.437902e-04 8.726253e-05 0 1
## 2 rs1011380 -5.231286e-05 8.881914e-05 0 1
## 3 rs1014986 7.976357e-05 1.177811e-03 0 1
## 4 rs10157808 2.620658e-04 2.047407e-05 0 1
## 5 rs1018102 4.699228e-05 1.019356e-04 0 1
## 6 rs10218584 -6.270159e-05 3.580614e-06 0 1This step generates a data frame in model.df that
includes all SNPs in both isoTWAS models at the locus and their effect
(non-zero or otherwise) on each isoforms we are fine-mapping at the
locus.
Now, we obtain our LD matrix for the SNPs in
model.df$SNP
V = readRDS(system.file("extdata",
"test_LD_finemapping.RDS",
package = "isotwas"))Now, we run the fine-mapping process in the following steps: (1) compute the correlation between the genetically-predicted expression of the isoforms, (2) enumerate the causal configurations, (3) compute marginal likelihoods and Bayes factors for each isoform, and (4) compute posterior inclusion probabilities and credible sets.
V = V[model.df$SNP,model.df$SNP]
Omega = Matrix::Matrix(as.matrix(model.df[,-c(1,ncol(model.df))]))
zscores = isoform_sig$Z
m = length(zscores)
### COMPUTE LD BETWEEN TX ON THE GENETIC LEVEL
wcor = isotwas::estimate_cor(as.matrix(Omega),
as.matrix(V),
intercept=T)[[1]]
diag(wcor) = 1
wcor[is.na(wcor)] = 0
### COMPUTE LD INTERCEPT BETWEEN ISOFRM ON THE GENETIC LEVEL
swld = isotwas::estimate_cor(as.matrix(Omega),
as.matrix(V),
intercept=T)[[2]]
null_res = m * log(1 - 1e-3)
marginal = m * log(1 - 1e-3)
comb_list = list()
for (n in 1:min(2,length(zscores))){
comb_list = c(comb_list,
combn(1:length(zscores),n,simplify=F))
}
pips = rep(0,length(zscores))
### COMPUTE BAYES FACTORS/LIKELIHOOD AT EACH CAUSAL CONFIG
for (j in 1:length(comb_list)){
subset = comb_list[[j]]
local = isotwas::bayes_factor(zscores,
idx_set = subset,
wcor = wcor)
marginal = log(exp(local) + exp(marginal))
for (idx in subset){
if (pips[idx] == 0){
pips[idx] = local
} else {
pips[idx] = log(exp(pips[idx]) + exp(local))
}
}
}
pips = exp(pips - marginal)
null_res = exp(null_res - marginal)
isoform_sig$pip = pips
isoform_sig = isoform_sig[order(isoform_sig$pip,decreasing = T),]
npost = isoform_sig$pip/sum(isoform_sig$pip)
csum = cumsum(npost)
isoform_sig$in_cred_set = F
for (i in 1:nrow(isoform_sig)){
isoform_sig$in_cred_set[i] = T
if (i > 1){
if (csum[i] > .9 & csum[i-1] < .9){
isoform_sig$in_cred_set[i] = T
}
if (csum[i] < .9){
isoform_sig$in_cred_set[i] = T
}
if (csum[i] > .9 & csum[i-1] > .9){
isoform_sig$in_cred_set[i] = F
}
}
}
print(isoform_sig)## Gene Chromosome Start End Feature Z
## 3 ENSG00000070831 1 22052627 22101360 ENST00000400259 3.344840
## 10 ENSG00000230068 1 22059197 22064199 ENST00000431803 -3.506525
## 8 ENSG00000162552 1 22117313 22143969 ENST00000290167 -3.441877
## P permute.P topSNP topSNP.P Screen.P
## 3 0.0008232993 0.011988012 rs11582542 0.0003396706 0.0016636386
## 10 0.0004539979 0.007992008 rs11582542 0.0003396706 0.0004539979
## 8 0.0005776926 0.041958042 rs11582542 0.0003396706 0.0008302377
## Screen.P.Adjusted Confirmation.P pip in_cred_set
## 3 0.001663639 0.0049397958 1.000000e+00 TRUE
## 10 0.001245357 0.0000000000 1.000000e+00 TRUE
## 8 0.001245357 0.0005776926 1.419146e-18 FALSEAs you can see, ENST00000431803 and ENST00000290167 have
PIP = 1 and are both in the 90% credible set
(in_cred_set column).
Full scripts for association testing genome-wide are available here.