Last updated: 2021-12-22

Checks: 6 1

Knit directory: cTWAS_analysis/

This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


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.

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

The command set.seed(20211220) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Using absolute paths to the files within your workflowr project makes it difficult for you and others to run your code on a different machine. Change the absolute path(s) below to the suggested relative path(s) to make your code more reproducible.

absolute relative
/project2/xinhe/shengqian/cTWAS/cTWAS_analysis/data/ data
/project2/xinhe/shengqian/cTWAS/cTWAS_analysis/code/ctwas_config.R code/ctwas_config.R

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 7e05bbc. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

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:


Untracked files:
    Untracked:  code/ctwas_config.R
    Untracked:  code/run_AF_analysis.sbatch
    Untracked:  code/run_AF_analysis.sh
    Untracked:  code/run_AF_ctwas_rss_LDR.R
    Untracked:  data/AF/
    Untracked:  data/mashr_Heart_Atrial_Appendage.db
    Untracked:  data/summary_known_genes_annotations.xlsx

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.


These are the previous versions of the repository in which changes were made to the R Markdown (analysis/Atrial_Fibrillation_Heart_Atrial_Appendage.Rmd) and HTML (docs/Atrial_Fibrillation_Heart_Atrial_Appendage.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File Version Author Date Message
Rmd 7e05bbc sq-96 2021-12-22 Start my new project
html 9cf0e70 sq-96 2021-12-22 Build site.
Rmd fbe71b7 sq-96 2021-12-22 Start my new project
html 77de9fb sq-96 2021-12-22 Build site.
Rmd c22d2e9 sq-96 2021-12-22 Start my new project
html 6b46de7 sq-96 2021-12-22 Build site.
Rmd 374327b sq-96 2021-12-22 Start my new project
html e22d0c9 sq-96 2021-12-21 Build site.
html a8dbae0 sq-96 2021-12-21 Build site.
Rmd 37c4593 sq-96 2021-12-21 publish files
html e080e7b sq-96 2021-12-20 Build site.
Rmd 0b4a886 sq-96 2021-12-20 Start my new project
html 9b1bb6e sq-96 2021-12-20 Build site.
html 773dcb5 sq-96 2021-12-20 Build site.
html 3052d49 sq-96 2021-12-20 Build site.
Rmd 12069cf sq-96 2021-12-20 Publish the files

Weight QC

TO-DO: add enhanced QC reporting (total number of weights, why each variant was missing for all genes)

qclist_all <- list()

qc_files <- paste0(results_dir, "/", list.files(results_dir, pattern="exprqc.Rd"))

for (i in 1:length(qc_files)){
  load(qc_files[i])
  chr <- unlist(strsplit(rev(unlist(strsplit(qc_files[i], "_")))[1], "[.]"))[1]
  qclist_all[[chr]] <- cbind(do.call(rbind, lapply(qclist,unlist)), as.numeric(substring(chr,4)))
}

qclist_all <- data.frame(do.call(rbind, qclist_all))
colnames(qclist_all)[ncol(qclist_all)] <- "chr"

rm(qclist, wgtlist, z_gene_chr)

#number of imputed weights
nrow(qclist_all)
[1] 10862
#number of imputed weights by chromosome
table(qclist_all$chr)

   1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
1103  776  625  427  501  624  518  377  399  423  648  611  203  368  356  495 
  17   18   19   20   21   22 
 676  177  840  318  126  271 
#proportion of imputed weights without missing variants
mean(qclist_all$nmiss==0)
[1] 0.6548518

Load ctwas results

Check convergence of parameters

library(ggplot2)
library(cowplot)

********************************************************
Note: As of version 1.0.0, cowplot does not change the
  default ggplot2 theme anymore. To recover the previous
  behavior, execute:
  theme_set(theme_cowplot())
********************************************************
load(paste0(results_dir, "/", analysis_id, "_ctwas.s2.susieIrssres.Rd"))
  
df <- data.frame(niter = rep(1:ncol(group_prior_rec), 2),
                 value = c(group_prior_rec[1,], group_prior_rec[2,]),
                 group = rep(c("Gene", "SNP"), each = ncol(group_prior_rec)))
df$group <- as.factor(df$group)

df$value[df$group=="SNP"] <- df$value[df$group=="SNP"]*thin #adjust parameter to account for thin argument

p_pi <- ggplot(df, aes(x=niter, y=value, group=group)) +
  geom_line(aes(color=group)) +
  geom_point(aes(color=group)) +
  xlab("Iteration") + ylab(bquote(pi)) +
  ggtitle("Prior mean") +
  theme_cowplot()

df <- data.frame(niter = rep(1:ncol(group_prior_var_rec), 2),
                 value = c(group_prior_var_rec[1,], group_prior_var_rec[2,]),
                 group = rep(c("Gene", "SNP"), each = ncol(group_prior_var_rec)))
df$group <- as.factor(df$group)
p_sigma2 <- ggplot(df, aes(x=niter, y=value, group=group)) +
  geom_line(aes(color=group)) +
  geom_point(aes(color=group)) +
  xlab("Iteration") + ylab(bquote(sigma^2)) +
  ggtitle("Prior variance") +
  theme_cowplot()

plot_grid(p_pi, p_sigma2)

Version Author Date
77de9fb sq-96 2021-12-22
6b46de7 sq-96 2021-12-22
a8dbae0 sq-96 2021-12-21
e080e7b sq-96 2021-12-20
3052d49 sq-96 2021-12-20
#estimated group prior
estimated_group_prior <- group_prior_rec[,ncol(group_prior_rec)]
names(estimated_group_prior) <- c("gene", "snp")
estimated_group_prior["snp"] <- estimated_group_prior["snp"]*thin #adjust parameter to account for thin argument
print(estimated_group_prior)
        gene          snp 
0.0208487453 0.0001927855 
#estimated group prior variance
estimated_group_prior_var <- group_prior_var_rec[,ncol(group_prior_var_rec)]
names(estimated_group_prior_var) <- c("gene", "snp")
print(estimated_group_prior_var)
    gene      snp 
9.646456 8.679052 
#report sample size
print(sample_size)
[1] 1030836
#report group size
group_size <- c(nrow(ctwas_gene_res), n_snps)
print(group_size)
[1]   10862 6839050
#estimated group PVE
estimated_group_pve <- estimated_group_prior_var*estimated_group_prior*group_size/sample_size #check PVE calculation
names(estimated_group_pve) <- c("gene", "snp")
print(estimated_group_pve)
      gene        snp 
0.00211918 0.01110076 
#compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))
[1] 0.01225055 0.12523867

Genes with highest PIPs

#distribution of PIPs
hist(ctwas_gene_res$susie_pip, xlim=c(0,1), main="Distribution of Gene PIPs")

Version Author Date
77de9fb sq-96 2021-12-22
6b46de7 sq-96 2021-12-22
a8dbae0 sq-96 2021-12-21
e080e7b sq-96 2021-12-20
3052d49 sq-96 2021-12-20
#genes with PIP>0.8 or 20 highest PIPs
head(ctwas_gene_res[order(-ctwas_gene_res$susie_pip),report_cols], max(sum(ctwas_gene_res$susie_pip>0.8), 20))
          genename region_tag susie_pip        mu2          PVE         z
2293          CAV1       7_70 1.0000000  622.84815 6.042165e-04 15.567870
3275         PRRX1       1_84 0.9998808  119.34663 1.157627e-04 14.667578
4009           DEK       6_14 0.9925129   63.47784 6.111795e-05 -9.000000
3527         CCND2       12_4 0.9906174   27.01637 2.596231e-05 -5.283784
1310           PXN      12_75 0.9826498   29.47170 2.809405e-05 -5.328302
12836 RP11-325L7.2       5_82 0.9826020 1993.73789 1.900449e-03 12.356322
9257         AGAP5      10_49 0.9779202   48.91516 4.640420e-05 11.518590
3523         KLF12      13_36 0.9770829   26.00017 2.464439e-05 -5.072464
6621         AKAP6       14_8 0.9721549   76.81846 7.244551e-05 -9.197368
6914          JAM2       21_9 0.9640543   22.19279 2.075505e-05  4.563232
9639         DLEU1      13_21 0.9586468   23.59091 2.193885e-05  4.697095
11818         DPF3      14_34 0.9574312   33.35512 3.097993e-05  6.264960
10521       FAM43A      3_120 0.9569176   29.69898 2.756935e-05 -5.487179
2444       SEC23IP      10_74 0.9517809   22.38862 2.067163e-05 -4.565228
13075    LINC01629      14_36 0.9471613   31.88580 2.929757e-05 -5.695652
2138           AES       19_4 0.9403219   20.20438 1.843030e-05  4.182804
4658        POPDC3       6_70 0.9252511   25.19514 2.261449e-05 -4.758170
10290       NKX2-5      5_103 0.9228745   61.55975 5.511247e-05 -9.391892
7515       TNFSF13       17_7 0.9135188   34.26969 3.036953e-05 -5.883117
5185          GYPC       2_74 0.9068050   38.95864 3.427111e-05 -6.380531
8248         CMTM5       14_3 0.9067300   30.09803 2.647442e-05 -5.472727
10548       SCN10A       3_28 0.8867624   77.72022 6.685774e-05 -8.814286
13967 RP5-890E16.5      17_28 0.8660551   23.11200 1.941751e-05 -4.761194
712          SP100      2_135 0.8658723   18.73540 1.573719e-05 -3.671335
9012         MTSS1       8_82 0.8594108   20.87861 1.740655e-05  4.402634
8992          MURC       9_50 0.8478869   23.34736 1.920376e-05  4.911964
5223         PSMB7       9_64 0.8430652   25.30057 2.069197e-05 -4.820896
6114       STK11IP      2_130 0.8406381   18.96825 1.546845e-05 -3.868022
10416          PGP       16_2 0.8298586   28.51218 2.295329e-05  5.943820
9691           BOK      2_144 0.8255975   19.18259 1.536335e-05  3.910125
8420          MARS      12_36 0.8180336   17.81958 1.414097e-05 -3.366197
3088          GNB4      3_110 0.8097696   30.34628 2.383842e-05 -5.583333
      num_eqtl
2293         3
3275         2
4009         1
3527         1
1310         1
12836        1
9257         2
3523         1
6621         1
6914         2
9639         1
11818        3
10521        1
2444         2
13075        1
2138         3
4658         1
10290        1
7515         1
5185         1
8248         1
10548        1
13967        1
712          2
9012         2
8992         2
5223         1
6114         2
10416        1
9691         3
8420         1
3088         1

Genes with largest effect sizes

#plot PIP vs effect size
plot(ctwas_gene_res$susie_pip, ctwas_gene_res$mu2, xlab="PIP", ylab="mu^2", main="Gene PIPs vs Effect Size")

Version Author Date
77de9fb sq-96 2021-12-22
6b46de7 sq-96 2021-12-22
a8dbae0 sq-96 2021-12-21
e080e7b sq-96 2021-12-20
3052d49 sq-96 2021-12-20
#genes with 20 largest effect sizes
head(ctwas_gene_res[order(-ctwas_gene_res$mu2),report_cols],20)
          genename region_tag    susie_pip       mu2          PVE           z
3250         WIPF1      2_105 0.000000e+00 2852.2491 0.000000e+00   8.4415584
12836 RP11-325L7.2       5_82 9.826020e-01 1993.7379 1.900449e-03  12.3563218
1487         SIRT1      10_44 3.337242e-01 1299.2406 4.206179e-04  -5.0538462
7413          IL6R       1_75 5.862422e-12 1273.1807 7.240649e-15  -4.9784946
6444         HERC4      10_44 1.013394e-01 1220.9395 1.200281e-04  -5.3595598
11106         NACA      12_35 5.127731e-03 1179.0569 5.865032e-06  -6.2408092
960          BAZ2A      12_35 2.047295e-03 1162.5562 2.308898e-06  -5.9428571
7986       SLC35A1       6_59 0.000000e+00  804.8763 0.000000e+00  -5.0724638
2507          WNT3      17_27 3.707228e-05  797.4584 2.867925e-08  -4.3894024
2292          CAV2       7_70 0.000000e+00  689.1121 0.000000e+00  14.5349429
2293          CAV1       7_70 1.000000e+00  622.8482 6.042165e-04  15.5678701
4999          ORC3       6_59 0.000000e+00  597.6005 0.000000e+00   4.1666667
10616       ARL17A      17_27 7.713881e-06  488.6853 3.656896e-09  -0.5399498
3687         KDM3B       5_82 5.275983e-06  436.2575 2.232835e-09  -6.6287277
7728        GPR155      2_105 0.000000e+00  371.9331 0.000000e+00  -1.4546586
7276      ARHGAP27      17_27 5.875984e-06  189.4420 1.079860e-09   0.1576482
10739         MAPT      17_27 2.791224e-04  166.5127 4.508711e-08   3.8724621
7729          PMVK       1_76 9.693335e-05  161.1464 1.515319e-08 -12.1029412
7730        PBXIP1       1_76 9.335873e-05  156.1414 1.414111e-08 -11.8676471
341         FAM13B       5_82 1.826415e-01  153.9825 2.728233e-05   7.6009242
      num_eqtl
3250         1
12836        1
1487         1
7413         1
6444         2
11106        2
960          1
7986         1
2507         2
2292         2
2293         3
4999         1
10616        2
3687         2
7728         2
7276         2
10739        2
7729         1
7730         1
341          2

Genes with highest PVE

#genes with 20 highest pve
head(ctwas_gene_res[order(-ctwas_gene_res$PVE),report_cols],20)
          genename region_tag susie_pip        mu2          PVE          z
12836 RP11-325L7.2       5_82 0.9826020 1993.73789 1.900449e-03  12.356322
2293          CAV1       7_70 1.0000000  622.84815 6.042165e-04  15.567870
1487         SIRT1      10_44 0.3337242 1299.24064 4.206179e-04  -5.053846
6444         HERC4      10_44 0.1013394 1220.93950 1.200281e-04  -5.359560
3275         PRRX1       1_84 0.9998808  119.34663 1.157627e-04  14.667578
6621         AKAP6       14_8 0.9721549   76.81846 7.244551e-05  -9.197368
10548       SCN10A       3_28 0.8867624   77.72022 6.685774e-05  -8.814286
12013       ZSWIM8      10_49 0.7522384   87.25963 6.367652e-05 -11.216495
4009           DEK       6_14 0.9925129   63.47784 6.111795e-05  -9.000000
8298       SYNPO2L      10_49 0.6265676   99.72467 6.061512e-05 -11.945652
10290       NKX2-5      5_103 0.9228745   61.55975 5.511247e-05  -9.391892
9257         AGAP5      10_49 0.9779202   48.91516 4.640420e-05  11.518590
5185          GYPC       2_74 0.9068050   38.95864 3.427111e-05  -6.380531
3665         KCNJ5      11_80 0.5037532   67.73277 3.309993e-05  -8.748092
11818         DPF3      14_34 0.9574312   33.35512 3.097993e-05   6.264960
10515      C5orf47      5_104 0.7661297   41.61108 3.092586e-05   6.680556
7515       TNFSF13       17_7 0.9135188   34.26969 3.036953e-05  -5.883117
13075    LINC01629      14_36 0.9471613   31.88580 2.929757e-05  -5.695652
1310           PXN      12_75 0.9826498   29.47170 2.809405e-05  -5.328302
10521       FAM43A      3_120 0.9569176   29.69898 2.756935e-05  -5.487179
      num_eqtl
12836        1
2293         3
1487         1
6444         2
3275         2
6621         1
10548        1
12013        1
4009         1
8298         1
10290        1
9257         2
5185         1
3665         1
11818        3
10515        1
7515         1
13075        1
1310         1
10521        1

Genes with largest z scores

#genes with 20 largest z scores
head(ctwas_gene_res[order(-abs(ctwas_gene_res$z)),report_cols],20)
          genename region_tag    susie_pip        mu2          PVE          z
2293          CAV1       7_70 1.000000e+00  622.84815 6.042165e-04  15.567870
3275         PRRX1       1_84 9.998808e-01  119.34663 1.157627e-04  14.667578
2292          CAV2       7_70 0.000000e+00  689.11215 0.000000e+00  14.534943
12836 RP11-325L7.2       5_82 9.826020e-01 1993.73789 1.900449e-03  12.356322
7729          PMVK       1_76 9.693335e-05  161.14638 1.515319e-08 -12.102941
8298       SYNPO2L      10_49 6.265676e-01   99.72467 6.061512e-05 -11.945652
7730        PBXIP1       1_76 9.335873e-05  156.14141 1.414111e-08 -11.867647
9805         MYOZ1      10_49 3.497772e-03   88.98387 3.019349e-07  11.759348
9257         AGAP5      10_49 9.779202e-01   48.91516 4.640420e-05  11.518590
12013       ZSWIM8      10_49 7.522384e-01   87.25963 6.367652e-05 -11.216495
6395        C9orf3       9_48 4.919856e-02   94.74228 4.521751e-06  10.671642
7407        ZBTB7B       1_76 9.636563e-05  122.09475 1.141378e-08  10.638889
13556   RP1-79C4.4       1_84 8.254209e-03   58.95140 4.720414e-07   9.586404
10290       NKX2-5      5_103 9.228745e-01   61.55975 5.511247e-05  -9.391892
9719        SEC24C      10_49 9.597200e-04   46.26092 4.306945e-08   9.271945
6621         AKAP6       14_8 9.721549e-01   76.81846 7.244551e-05  -9.197368
4009           DEK       6_14 9.925129e-01   63.47784 6.111795e-05  -9.000000
10548       SCN10A       3_28 8.867624e-01   77.72022 6.685774e-05  -8.814286
3665         KCNJ5      11_80 5.037532e-01   67.73277 3.309993e-05  -8.748092
7733         DCST2       1_76 2.613168e-04   89.76014 2.275418e-08  -8.718310
      num_eqtl
2293         3
3275         2
2292         2
12836        1
7729         1
8298         1
7730         1
9805         2
9257         2
12013        1
6395         1
7407         1
13556        2
10290        1
9719         2
6621         1
4009         1
10548        1
3665         1
7733         1

Comparing z scores and PIPs

#set nominal signifiance threshold for z scores
alpha <- 0.05

#bonferroni adjusted threshold for z scores
sig_thresh <- qnorm(1-(alpha/nrow(ctwas_gene_res)/2), lower=T)

#Q-Q plot for z scores
obs_z <- ctwas_gene_res$z[order(ctwas_gene_res$z)]
exp_z <- qnorm((1:nrow(ctwas_gene_res))/nrow(ctwas_gene_res))

plot(exp_z, obs_z, xlab="Expected z", ylab="Observed z", main="Gene z score Q-Q plot")
abline(a=0,b=1)

Version Author Date
77de9fb sq-96 2021-12-22
6b46de7 sq-96 2021-12-22
a8dbae0 sq-96 2021-12-21
e080e7b sq-96 2021-12-20
3052d49 sq-96 2021-12-20
#plot z score vs PIP
plot(abs(ctwas_gene_res$z), ctwas_gene_res$susie_pip, xlab="abs(z)", ylab="PIP")
abline(v=sig_thresh, col="red", lty=2)

Version Author Date
77de9fb sq-96 2021-12-22
6b46de7 sq-96 2021-12-22
a8dbae0 sq-96 2021-12-21
e080e7b sq-96 2021-12-20
3052d49 sq-96 2021-12-20
#proportion of significant z scores
mean(abs(ctwas_gene_res$z) > sig_thresh)
[1] 0.008654023
#genes with most significant z scores
head(ctwas_gene_res[order(-abs(ctwas_gene_res$z)),report_cols],20)
          genename region_tag    susie_pip        mu2          PVE          z
2293          CAV1       7_70 1.000000e+00  622.84815 6.042165e-04  15.567870
3275         PRRX1       1_84 9.998808e-01  119.34663 1.157627e-04  14.667578
2292          CAV2       7_70 0.000000e+00  689.11215 0.000000e+00  14.534943
12836 RP11-325L7.2       5_82 9.826020e-01 1993.73789 1.900449e-03  12.356322
7729          PMVK       1_76 9.693335e-05  161.14638 1.515319e-08 -12.102941
8298       SYNPO2L      10_49 6.265676e-01   99.72467 6.061512e-05 -11.945652
7730        PBXIP1       1_76 9.335873e-05  156.14141 1.414111e-08 -11.867647
9805         MYOZ1      10_49 3.497772e-03   88.98387 3.019349e-07  11.759348
9257         AGAP5      10_49 9.779202e-01   48.91516 4.640420e-05  11.518590
12013       ZSWIM8      10_49 7.522384e-01   87.25963 6.367652e-05 -11.216495
6395        C9orf3       9_48 4.919856e-02   94.74228 4.521751e-06  10.671642
7407        ZBTB7B       1_76 9.636563e-05  122.09475 1.141378e-08  10.638889
13556   RP1-79C4.4       1_84 8.254209e-03   58.95140 4.720414e-07   9.586404
10290       NKX2-5      5_103 9.228745e-01   61.55975 5.511247e-05  -9.391892
9719        SEC24C      10_49 9.597200e-04   46.26092 4.306945e-08   9.271945
6621         AKAP6       14_8 9.721549e-01   76.81846 7.244551e-05  -9.197368
4009           DEK       6_14 9.925129e-01   63.47784 6.111795e-05  -9.000000
10548       SCN10A       3_28 8.867624e-01   77.72022 6.685774e-05  -8.814286
3665         KCNJ5      11_80 5.037532e-01   67.73277 3.309993e-05  -8.748092
7733         DCST2       1_76 2.613168e-04   89.76014 2.275418e-08  -8.718310
      num_eqtl
2293         3
3275         2
2292         2
12836        1
7729         1
8298         1
7730         1
9805         2
9257         2
12013        1
6395         1
7407         1
13556        2
10290        1
9719         2
6621         1
4009         1
10548        1
3665         1
7733         1

Locus plots for genes and SNPs

ctwas_gene_res_sortz <- ctwas_gene_res[order(-abs(ctwas_gene_res$z)),]
report_cols_region <- report_cols[!(report_cols %in% c("num_eqtl"))]

n_plots <- 5
for (region_tag_plot in head(unique(ctwas_gene_res_sortz$region_tag), n_plots)){
  ctwas_res_region <-  ctwas_res[ctwas_res$region_tag==region_tag_plot,]
  start <- min(ctwas_res_region$pos)
  end <- max(ctwas_res_region$pos)
  
  ctwas_res_region <- ctwas_res_region[order(ctwas_res_region$pos),]
  ctwas_res_region_gene <- ctwas_res_region[ctwas_res_region$type=="gene",]
  ctwas_res_region_snp <- ctwas_res_region[ctwas_res_region$type=="SNP",]
  
  #region name
  print(paste0("Region: ", region_tag_plot))
  
  #table of genes in region
  print(ctwas_res_region_gene[,report_cols_region])
  
  par(mfrow=c(4,1))
  
  #gene z scores
  plot(ctwas_res_region_gene$pos, abs(ctwas_res_region_gene$z), xlab="Position", ylab="abs(gene_z)", xlim=c(start,end),
   ylim=c(0,max(sig_thresh, abs(ctwas_res_region_gene$z))),
   main=paste0("Region: ", region_tag_plot))
  abline(h=sig_thresh,col="red",lty=2)
  
  #significance threshold for SNPs
  alpha_snp <- 5*10^(-8)
  sig_thresh_snp <- qnorm(1-alpha_snp/2, lower=T)
  
  #snp z scores
  plot(ctwas_res_region_snp$pos, abs(ctwas_res_region_snp$z), xlab="Position", ylab="abs(snp_z)",xlim=c(start,end),
   ylim=c(0,max(sig_thresh_snp, max(abs(ctwas_res_region_snp$z)))))
  abline(h=sig_thresh_snp,col="purple",lty=2)
  
  #gene pips
  plot(ctwas_res_region_gene$pos, ctwas_res_region_gene$susie_pip, xlab="Position", ylab="Gene PIP", xlim=c(start,end), ylim=c(0,1))
  abline(h=0.8,col="blue",lty=2)
  
  #snp pips
  plot(ctwas_res_region_snp$pos, ctwas_res_region_snp$susie_pip, xlab="Position", ylab="SNP PIP", xlim=c(start,end), ylim=c(0,1))
  abline(h=0.8,col="blue",lty=2)
}
[1] "Region: 7_70"
       genename region_tag susie_pip        mu2          PVE          z
9316      LRRN3       7_70         0   4.798139 0.0000000000 -0.2554113
10486    IMMP2L       7_70         0   9.726840 0.0000000000 -1.1251408
10096    LSMEM1       7_70         0   4.653263 0.0000000000  0.1730349
6325    TMEM168       7_70         0  13.129092 0.0000000000  1.4473363
8015      GPR85       7_70         0  59.638319 0.0000000000  3.6640335
8014       BMT2       7_70         0  14.144825 0.0000000000  1.5304348
12439    HRAT17       7_70         0  10.070641 0.0000000000 -1.1617647
11985 LINC00998       7_70         0   9.138407 0.0000000000 -1.0595238
6892    PPP1R3A       7_70         0  11.619381 0.0000000000 -1.3141517
4991      MDFIC       7_70         0   9.711018 0.0000000000  1.1940299
4990        TES       7_70         0   5.365475 0.0000000000  0.2279955
2292       CAV2       7_70         0 689.112146 0.0000000000 14.5349429
2293       CAV1       7_70         1 622.848152 0.0006042165 15.5678701
11541    CAPZA2       7_70         0   4.817249 0.0000000000 -0.5522388

Version Author Date
77de9fb sq-96 2021-12-22
6b46de7 sq-96 2021-12-22
a8dbae0 sq-96 2021-12-21
e080e7b sq-96 2021-12-20
3052d49 sq-96 2021-12-20
[1] "Region: 1_84"
        genename region_tag   susie_pip        mu2          PVE          z
3275       PRRX1       1_84 0.999880773 119.346630 1.157627e-04 14.6675780
13556 RP1-79C4.4       1_84 0.008254209  58.951404 4.720414e-07  9.5864041
132         FMO3       1_84 0.007235866   6.979886 4.899472e-08  0.7058824
1450        FMO2       1_84 0.005819057   4.561901 2.575187e-08 -0.1257130
964         FMO4       1_84 0.006067492   5.455974 3.211382e-08  0.8619618
364         MYOC       1_84 0.005790338   4.650929 2.612486e-08  0.2537313
3427       VAMP4       1_84 0.024556213  14.423816 3.435991e-07 -1.8006622
11344       DNM3       1_84 0.028756639  18.635294 5.198581e-07  1.9131944
5075        PIGC       1_84 0.009407288   9.222146 8.416022e-08 -1.1709099
10093   C1orf105       1_84 0.006032874   4.956744 2.900889e-08 -0.5632184
1451        SUCO       1_84 0.010489256   9.463884 9.629961e-08  0.9829060

Version Author Date
77de9fb sq-96 2021-12-22
[1] "Region: 5_82"
          genename region_tag    susie_pip         mu2          PVE           z
1029        PKD2L2       5_82 1.074133e-06   44.303563 4.616439e-11 -3.01646091
12836 RP11-325L7.2       5_82 9.826020e-01 1993.737895 1.900449e-03 12.35632184
341         FAM13B       5_82 1.826415e-01  153.982522 2.728233e-05  7.60092424
2972          NME5       5_82 8.035979e-08    7.216916 5.626015e-13 -0.78899083
3687         KDM3B       5_82 5.275983e-06  436.257478 2.232835e-09 -6.62872770
4678         REEP2       5_82 3.273136e-07   34.258854 1.087795e-11 -2.33668342
3688          EGR1       5_82 1.582507e-07   28.468689 4.370424e-12 -0.46666667
2975         HSPA9       5_82 1.191309e-07   15.195222 1.756070e-12 -1.47659574
425         CTNNA1       5_82 8.323910e-08   55.376623 4.471614e-12 -0.68493151
6256        LRRTM2       5_82 8.644342e-08   53.238157 4.464423e-12 -0.57534247
8964       SLC23A1       5_82 7.435045e-08    9.459617 6.822878e-13 -0.45945946
12284        PROB1       5_82 7.443824e-08    7.340816 5.300915e-13 -0.37662338
8961       SPATA24       5_82 7.971111e-08   12.023543 9.297404e-13  0.74647887
8959       DNAJC18       5_82 7.430306e-08    6.843115 4.932544e-13  0.07692308
10457      TMEM173       5_82 1.052563e-07    9.743512 9.948880e-13  0.52710454

Version Author Date
77de9fb sq-96 2021-12-22
[1] "Region: 1_76"
      genename region_tag    susie_pip        mu2          PVE           z
7729      PMVK       1_76 9.693335e-05 161.146376 1.515319e-08 -12.1029412
7730    PBXIP1       1_76 9.335873e-05 156.141409 1.414111e-08 -11.8676471
7409      SHC1       1_76 1.029606e-04  12.706941 1.269178e-09  -0.5945946
7407    ZBTB7B       1_76 9.636563e-05 122.094754 1.141378e-08  10.6388889
6001    ADAM15       1_76 1.054076e-04  19.807533 2.025408e-09   3.0000000
7733     DCST2       1_76 2.613168e-04  89.760135 2.275418e-08  -8.7183099
6011     EFNA3       1_76 1.089596e-04  11.240179 1.188089e-09  -2.2724138
8787     EFNA1       1_76 1.606785e-04  11.329113 1.765892e-09  -1.8067354
8786   SLC50A1       1_76 1.003646e-04   6.652241 6.476776e-10   0.8410502
10571     MUC1       1_76 2.518113e-04  20.984333 5.126026e-09  -0.9716981
9323      MTX1       1_76 1.791435e-04  16.586200 2.882427e-09   0.5555556
9782       GBA       1_76 1.026910e-04   7.311893 7.284044e-10  -0.5299145
3310    SCAMP3       1_76 2.912632e-04  23.792281 6.722520e-09  -3.3289474
7735    YY1AP1       1_76 2.425271e-04  16.618716 3.909923e-09  -1.9970399
4695      DAP3       1_76 2.664936e-04  29.575646 7.645948e-09  -3.5890411
3313     GON4L       1_76 2.496430e-03  41.586291 1.007117e-07  -3.9452055
4703     SYT11       1_76 2.496430e-03  41.586291 1.007117e-07   3.9452055
6015      RIT1       1_76 2.304738e-04  20.358641 4.551775e-09   2.5915493
4696  KIAA0907       1_76 1.475526e-04   7.222145 1.033769e-09   0.5050776
3314   ARHGEF2       1_76 2.545536e-03  38.509404 9.509474e-08   3.1323529
7755      SSR2       1_76 1.009922e-04  10.918250 1.069673e-09   1.8503401
3315   LAMTOR2       1_76 2.403669e-04  11.662362 2.719390e-09   0.7234637
11033   SEMA4A       1_76 3.065448e-04  12.239053 3.639588e-09  -0.4358974
7419  SLC25A44       1_76 1.133154e-04   6.396866 7.031801e-10   0.9130609
12648    BGLAP       1_76 1.086835e-04  12.509241 1.318879e-09   2.3348238
7418     PAQR6       1_76 1.717449e-04  10.265336 1.710281e-09   0.2561399
7754    TMEM79       1_76 9.629572e-05   9.278082 8.667136e-10   2.2625000
11561     SMG5       1_76 9.629572e-05   9.278082 8.667136e-10  -2.2625000

Version Author Date
77de9fb sq-96 2021-12-22
[1] "Region: 10_49"
            genename region_tag    susie_pip       mu2          PVE           z
2448           PALD1      10_49 6.930061e-05 10.464471 7.035011e-10   1.2023810
5427        ADAMTS14      10_49 4.041175e-05  5.577877 2.186689e-10  -0.5051546
8279           SGPL1      10_49 4.383695e-05  6.314746 2.685386e-10   0.6594096
8282           PCBD1      10_49 3.637634e-05  4.624703 1.631974e-10  -0.1515152
2449           UNC5B      10_49 4.079184e-05  5.662476 2.240733e-10   0.5252525
2450           CDH23      10_49 5.372325e-05  8.157613 4.251437e-10  -0.9400479
2451            VSIR      10_49 4.235290e-05  5.884211 2.417585e-10  -0.4533295
12014      C10orf105      10_49 3.693800e-05  4.763447 1.706889e-10   0.2383178
11301           PSAP      10_49 4.143656e-05  5.328503 2.141901e-10   0.6926829
3840           CHST3      10_49 3.920994e-05  6.023333 2.291097e-10  -1.2205882
2452          SPOCK2      10_49 3.803373e-05  5.050485 1.863427e-10   0.5000000
8296         ANAPC16      10_49 7.630441e-05 11.415357 8.449861e-10   1.5096021
5424           ASCC1      10_49 3.691548e-05  5.519486 1.976595e-10  -0.3622163
8632           DDIT4      10_49 4.114673e-05  6.155464 2.457008e-10   0.4788732
6453         DNAJB12      10_49 3.871607e-05  7.240038 2.719209e-10  -1.2362964
7004             MCU      10_49 3.659951e-05  5.173114 1.836698e-10   0.6527778
5426            OIT3      10_49 4.946253e-05  8.617426 4.134894e-10  -0.6470588
13645  RP11-344N10.5      10_49 3.675546e-05  4.856516 1.731638e-10   0.8848495
8299          NUDT13      10_49 4.246247e-05 16.333612 6.728184e-10  -4.1574570
10190         MRPS16      10_49 8.192817e-05 15.875112 1.261713e-09   2.4854195
11923         DNAJC9      10_49 1.447809e-04 21.802775 3.062200e-09  -2.9593023
7006          CFAP70      10_49 4.129670e-05 16.906306 6.772897e-10  -4.3161765
5422           ANXA7      10_49 4.129111e-05 16.129161 6.460689e-10  -4.2137405
13668  RP11-464F9.22      10_49 5.378631e-05 19.460835 1.015415e-09   4.8622754
9805           MYOZ1      10_49 3.497772e-03 88.983867 3.019349e-07  11.7593478
9257           AGAP5      10_49 9.779202e-01 48.915157 4.640420e-05  11.5185898
8298         SYNPO2L      10_49 6.265676e-01 99.724672 6.061512e-05 -11.9456522
13589 RP11-574K11.29      10_49 5.552979e-05 15.390089 8.290441e-10   1.2828283
9719          SEC24C      10_49 9.597200e-04 46.260923 4.306945e-08   9.2719446
11176          FUT11      10_49 5.375832e-04 21.866281 1.140331e-08  -6.5522388
12013         ZSWIM8      10_49 7.522384e-01 87.259631 6.367652e-05 -11.2164948
3838            PLAU      10_49 2.066160e-04 15.107410 3.028060e-09   3.5310286
10503          AP3M1      10_49 3.334117e-04 23.844396 7.712187e-09  -2.0202020
7013             ADK      10_49 3.465243e-04 24.183065 8.129342e-09   2.0408163
1049          DUSP13      10_49 5.692456e-05  8.684736 4.795862e-10  -1.7916667
8167          COMTD1      10_49 1.542145e-04 10.638192 1.591489e-09   2.1111111
13481 RP11-399K21.11      10_49 2.039823e-04 18.379908 3.637025e-09   0.6792247
6445        C10orf11      10_49 4.685159e-05  6.813717 3.096841e-10  -1.0581482
13654 RP11-399K21.14      10_49 3.636678e-05  4.614057 1.627790e-10  -0.1324675

Version Author Date
77de9fb sq-96 2021-12-22

SNPs with highest PIPs

#snps with PIP>0.8 or 20 highest PIPs
head(ctwas_snp_res[order(-ctwas_snp_res$susie_pip),report_cols_snps],
max(sum(ctwas_snp_res$susie_pip>0.8), 20))
                id region_tag susie_pip        mu2          PVE           z
28439   rs12404927       1_75 1.0000000 1437.87698 1.394865e-03  -1.5517241
28440    rs7536152       1_75 1.0000000 1473.01903 1.428956e-03  -6.4626866
28533   rs34515871       1_76 1.0000000  299.28422 2.903316e-04  17.0000000
32410   rs12142529       1_83 1.0000000  145.61337 1.412576e-04  12.0396040
32435  rs112797273       1_83 1.0000000  101.67077 9.862944e-05   9.4903846
183533   rs1906615       4_72 1.0000000 1756.73093 1.704181e-03  45.1604938
183538   rs7440714       4_72 1.0000000  992.10110 9.624238e-04  -9.2183406
272612 rs111990232       6_59 1.0000000 1710.90898 1.659730e-03   0.5392670
414968  rs76106073      10_44 1.0000000 1273.72173 1.235620e-03   0.6572238
550561 rs140798119      15_35 1.0000000  285.26104 2.767279e-04   3.9950739
550564  rs74022964      15_35 1.0000000  342.92230 3.326643e-04  12.5777778
681435   rs4972703      2_105 1.0000000 4299.78452 4.171163e-03  -0.7842670
739005   rs9770220       7_70 1.0000000  835.76918 8.107683e-04  -5.2814815
788780  rs61937778      12_35 1.0000000 1668.54978 1.618637e-03  -1.3056995
28537    rs2878412       1_76 1.0000000  105.57502 1.024169e-04   5.4250000
586354  rs62056842      17_27 1.0000000 1130.21900 1.096410e-03   1.1269036
28535   rs11576820       1_76 1.0000000  133.71597 1.297160e-04   8.0243902
422177  rs60469668      10_66 1.0000000  172.10824 1.669599e-04  16.9647059
110779   rs9852222        3_9 1.0000000   83.84308 8.133503e-05   8.3823529
339255    rs373983       8_21 0.9999998   67.64329 6.561982e-05   8.5588235
298378  rs12112152       7_15 0.9999995   45.13406 4.378392e-05   6.9358974
422175   rs7094488      10_66 0.9999992  110.01432 1.067233e-04  14.1060606
28540   rs10908445       1_76 0.9999991  193.56990 1.877794e-04 -13.1470588
570907   rs6499606      16_39 0.9999990  103.94144 1.008321e-04  12.6376812
826958   rs8005417       14_9 0.9999978  427.56249 4.147716e-04   0.9863946
852634 rs140185678       16_2 0.9999910   45.86064 4.448839e-05   7.6100917
97513    rs1975584      2_118 0.9999847   71.76048 6.961280e-05   1.2311828
580602  rs72811292      17_11 0.9999789   40.04445 3.884576e-05  -6.5233645
279525  rs11756438       6_79 0.9999723   47.69440 4.626641e-05  -8.2089552
117531 rs116202356       3_27 0.9999305   36.10604 3.502355e-05   6.0813559
176218   rs1458038       4_54 0.9998791   34.99562 3.394467e-05   6.0277778
121102  rs12330500       3_40 0.9997290  519.11865 5.034535e-04   0.5641026
396844   rs1886296       9_73 0.9995279   34.54923 3.349992e-05   5.9863014
39581    rs6427989      1_102 0.9993189   36.39724 3.528442e-05  -4.4142857
653873    rs464901       22_4 0.9992684   43.84562 4.250292e-05  -7.0555556
708945 rs574293775       5_82 0.9987918 1853.69252 1.796069e-03   6.2411765
39582   rs12353975      1_102 0.9982647   35.87479 3.474126e-05  -4.3000000
184120   rs7700110       4_73 0.9979235   45.19125 4.374839e-05   7.0666667
232199 rs199992924       5_68 0.9969497 1340.40151 1.296339e-03  -0.6309524
655632    rs133902       22_7 0.9963501   31.05883 3.001978e-05   6.1617647
232019    rs338623       5_68 0.9953592   67.12738 6.481716e-05  -7.9855072
414975  rs12360521      10_44 0.9953399 1283.59990 1.239400e-03   5.3740458
696421   rs7374540       3_28 0.9939021   45.86079 4.421764e-05   4.7794118
570915    rs876727      16_39 0.9926586   48.82838 4.702000e-05 -10.0000000
275680   rs9496567       6_67 0.9918710   26.51232 2.551017e-05   5.1125000
586570  rs75230966      17_27 0.9917637   49.43748 4.756362e-05   6.0582524
183553   rs4631108       4_72 0.9915393  306.03072 2.943645e-04 -17.6376812
40258    rs4951023      1_104 0.9906227   31.13453 2.991996e-05   5.5970149
238444   rs6894302       5_84 0.9877503   34.35180 3.291601e-05  -7.2933333
666039   rs2885697       1_25 0.9874950   44.28194 4.242013e-05  -6.2714286
570920  rs60602157      16_39 0.9854801   38.72262 3.701885e-05  10.4415584
339226    rs403894       8_19 0.9793925   56.27052 5.346236e-05   6.8529412
650809   rs7282237      21_16 0.9782270   25.00237 2.372636e-05  -4.7605634
529120   rs2738413      14_29 0.9782089  112.52732 1.067825e-04 -11.6119403
570912   rs4788691      16_39 0.9747904   40.13308 3.795108e-05  -1.6900000
329588  rs35760656       7_94 0.9746629   35.32804 3.340291e-05  -6.3536585
470281  rs17380837      12_18 0.9737999   44.42786 4.196967e-05  -6.9583333
366399   rs7460121       8_88 0.9725824   24.19453 2.282727e-05   4.7500000
2727      rs284278        1_7 0.9713320   35.00958 3.298869e-05  -6.0857143
826964   rs8011559       14_9 0.9693907  435.54854 4.095867e-04   5.6913580
283467    rs958747       6_89 0.9688123   24.37662 2.290992e-05  -4.9710145
556130  rs12898337      15_48 0.9643924   35.44339 3.315885e-05  -6.0579710
279482  rs77435894       6_78 0.9642414   28.80796 2.694690e-05   6.1111111
549931    rs745636      15_33 0.9571831   25.61925 2.378876e-05  -5.0740741
25783    rs4839174       1_69 0.9569309   33.97079 3.153528e-05   5.9459459
183505   rs1823291       4_72 0.9568146  286.06204 2.655207e-04 -19.4520548
604563  rs17794590      18_24 0.9550117   24.59348 2.278448e-05  -4.8242424
100283  rs35880620      2_125 0.9544571   39.77143 3.682460e-05  -6.2168675
58743    rs7578482       2_16 0.9542037   32.14647 2.975670e-05   5.7435897
738958   rs1633714       7_70 0.9499686  778.22440 7.171740e-04  11.7536232
97488   rs11889306      2_118 0.9487414   45.59873 4.196730e-05   6.7391304
826887  rs73241997       14_9 0.9472521  276.10117 2.537139e-04   7.8817204
681445   rs1367220      2_105 0.9467808 4303.78606 3.952852e-03   7.2207792
784225   rs2291437      12_17 0.9462388   66.11326 6.068757e-05   9.1826923
32457    rs7522387       1_83 0.9460866   59.86066 5.493926e-05   6.6375000
696586   rs7373065       3_28 0.9355464   65.79596 5.971384e-05  -8.0637450
28495     rs906280       1_75 0.9353945   59.57170 5.405616e-05   7.9111111
490521  rs12425471      12_69 0.9303264   41.49244 3.744680e-05  -5.2500000
614935   rs7256735       19_3 0.9241700   23.45925 2.103180e-05   4.4952381
279432     rs89107       6_78 0.9189049   47.36791 4.222457e-05  -9.5606061
372071   rs1594768        9_9 0.9139999   24.21822 2.147330e-05  -4.6911765
5627    rs10917072       1_15 0.9043759   33.53922 2.942472e-05   5.8235294
217602 rs114414434       5_31 0.8994692   26.16276 2.282866e-05  -4.9718310
232192   rs4073838       5_68 0.8923424 1368.92743 1.185011e-03  -5.9545455
272609    rs371814       6_59 0.8922509 1722.76723 1.491159e-03   6.1911765
252160  rs73724866       6_13 0.8855608   94.00847 8.075990e-05 -10.4141414
479761    rs776211      12_43 0.8831880   28.12821 2.409937e-05  -6.3424658
550553    rs519946      15_34 0.8819650   34.47055 2.949239e-05   5.9242424
494561   rs6560886      12_82 0.8817026   28.94415 2.475673e-05   5.6666667
28531   rs12128882       1_76 0.8779576   32.56646 2.773668e-05   1.0731707
68039     rs243080       2_40 0.8749115   24.16683 2.051135e-05  -4.5522388
467254  rs12821447      12_12 0.8716947   24.08112 2.036346e-05   4.5128205
615414  rs73919353       19_5 0.8614285   28.38347 2.371894e-05   5.4076433
453904    rs565449      11_55 0.8575434   28.38859 2.361622e-05  -5.2535211
123427   rs1091584       3_45 0.8532454   30.27502 2.505929e-05   5.5882353
244628  rs62377226      5_100 0.8422335   33.87890 2.768039e-05   5.8493151
113948  rs73032373       3_18 0.8375231   33.29814 2.705373e-05  -5.9178082
369034  rs72693377       8_94 0.8312985   24.72583 1.993968e-05   4.2762431
659021  rs11705586      22_15 0.8290134   26.32683 2.117242e-05   4.8988764
121105      rs6924       3_40 0.8277296  522.73041 4.197364e-04   4.2133333
304278 rs145593380       7_26 0.8170209   29.74581 2.357596e-05   5.3333333
495691   rs7321083       13_3 0.8149736   24.65878 1.949510e-05   4.8292683

SNPs with largest effect sizes

#plot PIP vs effect size
#plot(ctwas_snp_res$susie_pip, ctwas_snp_res$mu2, xlab="PIP", ylab="mu^2", main="SNP PIPs vs Effect Size")

#SNPs with 50 largest effect sizes
head(ctwas_snp_res[order(-ctwas_snp_res$mu2),report_cols_snps],50)
                id region_tag    susie_pip      mu2          PVE         z
681445   rs1367220      2_105 9.467808e-01 4303.786 3.952852e-03  7.220779
681435   rs4972703      2_105 1.000000e+00 4299.785 4.171163e-03 -0.784267
681446   rs1367219      2_105 2.069799e-01 4294.507 8.622874e-04  7.103896
681441  rs10168156      2_105 3.312313e-01 4292.416 1.379252e-03  7.233766
681439   rs6713018      2_105 9.687858e-02 4291.538 4.033213e-04  7.168831
681430   rs1864453      2_105 1.876497e-01 4289.961 7.809293e-04  7.207792
681422   rs2033315      2_105 1.145049e-01 4287.502 4.762543e-04  7.064935
681433   rs6707162      2_105 3.496928e-02 4279.920 1.451887e-04  7.090909
681434   rs6735680      2_105 2.455791e-02 4279.461 1.019509e-04  7.077922
681423   rs2033314      2_105 2.432678e-02 4278.387 1.009660e-04  7.155844
681425  rs28485554      2_105 8.747105e-03 4277.770 3.629879e-05  7.000000
681456  rs10803884      2_105 7.083440e-03 4217.525 2.898093e-05  7.285714
681464   rs6738901      2_105 7.743316e-03 4201.159 3.155778e-05 -7.272727
681449  rs12466643      2_105 3.135672e-04 4196.140 1.276412e-06  7.272727
681480  rs35368253      2_105 1.998443e-04 4178.444 8.100592e-07 -7.256410
681466  rs10197521      2_105 1.763290e-11 4171.648 7.135785e-14 -6.547619
681468  rs34661753      2_105 3.194389e-08 4119.395 1.276532e-10 -7.076923
681419   rs7590328      2_105 2.272449e-11 4103.428 9.045890e-14  6.784810
681417   rs6433497      2_105 0.000000e+00 3465.041 0.000000e+00  6.345679
681418   rs2115874      2_105 0.000000e+00 3462.201 0.000000e+00  6.419753
681410   rs1430185      2_105 0.000000e+00 3387.075 0.000000e+00  6.185185
681397   rs1991601      2_105 0.000000e+00 3363.058 0.000000e+00  5.950617
681395   rs6759870      2_105 0.000000e+00 2978.209 0.000000e+00  5.987342
681486  rs35215597      2_105 0.000000e+00 2793.497 0.000000e+00 -8.441558
681506  rs56181519      2_105 0.000000e+00 2703.083 0.000000e+00 -8.597403
681308  rs13032076      2_105 0.000000e+00 2402.122 0.000000e+00 -5.250000
681408  rs13024657      2_105 0.000000e+00 2366.705 0.000000e+00 -5.258065
681509   rs2358891      2_105 0.000000e+00 2265.703 0.000000e+00 -7.500000
681526  rs35808589      2_105 0.000000e+00 2030.164 0.000000e+00 -5.229885
681558   rs1376875      2_105 0.000000e+00 2019.770 0.000000e+00 -5.000000
681541  rs35444726      2_105 0.000000e+00 2016.551 0.000000e+00 -5.172414
681550  rs13420492      2_105 0.000000e+00 1989.902 0.000000e+00 -4.908046
681552  rs71417497      2_105 0.000000e+00 1980.143 0.000000e+00 -4.772727
708976   rs2040862       5_82 1.932759e-01 1961.536 3.677770e-04 12.459770
708860  rs17171711       5_82 5.061870e-02 1958.969 9.619419e-05 12.482759
708884  rs13355516       5_82 4.398123e-03 1956.728 8.348495e-06 12.379310
708835   rs9327807       5_82 2.294241e-03 1953.970 4.348779e-06 12.367816
708869  rs77915370       5_82 3.425950e-04 1952.684 6.489682e-07 12.356322
708740  rs78081438       5_82 6.691465e-04 1951.053 1.266487e-06 12.390805
708691  rs73300168       5_82 3.565820e-04 1950.786 6.748068e-07 12.298851
708817  rs11959181       5_82 9.541911e-05 1949.634 1.804674e-07 12.241379
708671  rs10076361       5_82 2.676249e-04 1949.238 5.060598e-07 12.298851
708836 rs113331339       5_82 8.891833e-05 1949.202 1.681352e-07 12.241379
708789  rs73299268       5_82 8.169862e-05 1949.085 1.544742e-07 12.241379
708892 rs148378888       5_82 1.992837e-04 1948.655 3.767186e-07 12.344828
708714  rs73299210       5_82 5.888915e-05 1946.246 1.111843e-07 12.287356
708717  rs73299219       5_82 6.204555e-05 1945.825 1.171183e-07 12.287356
708737  rs13362264       5_82 2.545517e-05 1945.354 4.803802e-08 12.264368
681335   rs2303891      2_105 0.000000e+00 1921.675 0.000000e+00 -4.185185
708899  rs10479176       5_82 6.731608e-01 1866.560 1.218908e-03  6.931818

SNPs with highest PVE

#SNPs with 50 highest pve
head(ctwas_snp_res[order(-ctwas_snp_res$PVE),report_cols_snps],50)
                id region_tag  susie_pip       mu2          PVE           z
681435   rs4972703      2_105 1.00000000 4299.7845 0.0041711625  -0.7842670
681445   rs1367220      2_105 0.94678083 4303.7861 0.0039528520   7.2207792
708945 rs574293775       5_82 0.99879179 1853.6925 0.0017960693   6.2411765
183533   rs1906615       4_72 1.00000000 1756.7309 0.0017041808  45.1604938
272612 rs111990232       6_59 1.00000000 1710.9090 0.0016597296   0.5392670
788780  rs61937778      12_35 1.00000000 1668.5498 0.0016186375  -1.3056995
272609    rs371814       6_59 0.89225093 1722.7672 0.0014911593   6.1911765
28440    rs7536152       1_75 1.00000000 1473.0190 0.0014289557  -6.4626866
28439   rs12404927       1_75 1.00000000 1437.8770 0.0013948649  -1.5517241
681441  rs10168156      2_105 0.33123129 4292.4155 0.0013792517   7.2337662
232199 rs199992924       5_68 0.99694972 1340.4015 0.0012963390  -0.6309524
414975  rs12360521      10_44 0.99533986 1283.5999 0.0012394000   5.3740458
414968  rs76106073      10_44 1.00000000 1273.7217 0.0012356201   0.6572238
708899  rs10479176       5_82 0.67316075 1866.5597 0.0012189085   6.9318182
232192   rs4073838       5_68 0.89234244 1368.9274 0.0011850110  -5.9545455
586354  rs62056842      17_27 1.00000000 1130.2190 0.0010964101   1.1269036
788720   rs2860482      12_35 0.64493512 1697.2111 0.0010618479  -7.1052632
788775   rs7313074      12_35 0.58220672 1705.1065 0.0009630285  -6.9868421
183538   rs7440714       4_72 1.00000000  992.1011 0.0009624238  -9.2183406
681446   rs1367219      2_105 0.20697994 4294.5073 0.0008622874   7.1038961
739005   rs9770220       7_70 1.00000000  835.7692 0.0008107683  -5.2814815
681430   rs1864453      2_105 0.18764974 4289.9609 0.0007809293   7.2077922
738958   rs1633714       7_70 0.94996864  778.2244 0.0007171740  11.7536232
788778   rs4759256      12_35 0.41979796 1704.6254 0.0006941921  -6.9736842
272613    rs384318       6_59 0.37976257 1721.2083 0.0006340975   6.1029412
738995  rs10255816       7_70 0.78914390  709.0714 0.0005428210  -9.7123288
121102  rs12330500       3_40 0.99972902  519.1187 0.0005034535   0.5641026
681422   rs2033315      2_105 0.11450492 4287.5023 0.0004762543   7.0649351
272614   rs1145714       6_59 0.26478699 1717.6032 0.0004411943   6.1029412
121105      rs6924       3_40 0.82772961  522.7304 0.0004197364   4.2133333
826958   rs8005417       14_9 0.99999777  427.5625 0.0004147716   0.9863946
826964   rs8011559       14_9 0.96939072  435.5485 0.0004095867   5.6913580
681439   rs6713018      2_105 0.09687858 4291.5379 0.0004033213   7.1688312
232191   rs4235764       5_68 0.29136909 1367.2078 0.0003864457  -5.8787879
708976   rs2040862       5_82 0.19327595 1961.5362 0.0003677770  12.4597701
788715   rs7978685      12_35 0.21731446 1694.8625 0.0003573004  -7.0657895
550564  rs74022964      15_35 1.00000000  342.9223 0.0003326643  12.5777778
232194   rs4235768       5_68 0.24121476 1367.3620 0.0003199616  -5.8030303
272571   rs9444476       6_59 0.19008883 1698.6999 0.0003132447  -6.3088235
183553   rs4631108       4_72 0.99153932  306.0307 0.0002943645 -17.6376812
28533   rs34515871       1_76 1.00000000  299.2842 0.0002903316  17.0000000
272596   rs9444488       6_59 0.16974071 1703.6881 0.0002805347  -6.2647059
550561 rs140798119      15_35 1.00000000  285.2610 0.0002767279   3.9950739
183505   rs1823291       4_72 0.95681461  286.0620 0.0002655207 -19.4520548
826887  rs73241997       14_9 0.94725212  276.1012 0.0002537139   7.8817204
708944 rs191830974       5_82 0.12101272 1865.7774 0.0002190288   6.9206349
121104   rs1554125       3_40 0.38957721  520.5125 0.0001967139   4.1333333
28540   rs10908445       1_76 0.99999905  193.5699 0.0001877794 -13.1470588
272583   rs9444481       6_59 0.10849613 1701.3375 0.0001790668  -6.2500000
422177  rs60469668      10_66 1.00000000  172.1082 0.0001669599  16.9647059

SNPs with largest z scores

#histogram of (abs) SNP z scores
hist(abs(ctwas_snp_res$z))

Version Author Date
77de9fb sq-96 2021-12-22
#SNPs with 50 largest z scores
head(ctwas_snp_res[order(-abs(ctwas_snp_res$z)),report_cols_snps],50)
                id region_tag    susie_pip       mu2          PVE         z
183533   rs1906615       4_72 1.000000e+00 1756.7309 1.704181e-03  45.16049
183534  rs75725917       4_72 0.000000e+00  588.7458 0.000000e+00  42.33663
183532   rs1906611       4_72 0.000000e+00  526.9126 0.000000e+00  40.97143
183531  rs28521134       4_72 0.000000e+00  512.6248 0.000000e+00  40.22772
183527  rs10019689       4_72 0.000000e+00  508.5039 0.000000e+00  40.11000
183528  rs76013973       4_72 0.000000e+00  502.3963 0.000000e+00  40.04902
183530  rs12639820       4_72 0.000000e+00  499.7947 0.000000e+00  40.03922
183526  rs12647393       4_72 0.000000e+00  500.8094 0.000000e+00  39.98020
183529  rs74496596       4_72 0.000000e+00  484.7331 0.000000e+00  39.75728
183525  rs12647316       4_72 0.000000e+00  476.7021 0.000000e+00  39.55882
183524   rs4529121       4_72 0.000000e+00  480.5177 0.000000e+00  39.55340
183518  rs12650829       4_72 0.000000e+00  450.0533 0.000000e+00  31.54444
183541   rs3866831       4_72 0.000000e+00  919.4634 0.000000e+00 -30.10000
183540   rs6533530       4_72 0.000000e+00  904.8879 0.000000e+00 -29.91429
183523  rs12644107       4_72 0.000000e+00  153.0712 0.000000e+00  25.98851
183521   rs2723318       4_72 0.000000e+00  365.6985 0.000000e+00 -25.45833
183516   rs2218698       4_72 0.000000e+00  337.5178 0.000000e+00 -25.09589
183519   rs2197814       4_72 0.000000e+00  339.8970 0.000000e+00  24.59722
183517   rs1448799       4_72 0.000000e+00  339.2379 0.000000e+00  24.56944
183515 rs112927894       4_72 0.000000e+00  320.7713 0.000000e+00  24.50685
183505   rs1823291       4_72 9.568146e-01  286.0620 2.655207e-04 -19.45205
183509   rs2723296       4_72 4.318269e-02  278.0630 1.164832e-05 -19.02740
183514  rs11724067       4_72 0.000000e+00  321.1887 0.000000e+00 -18.94059
183508   rs2044674       4_72 2.641186e-06  253.4848 6.494734e-10  18.58108
183553   rs4631108       4_72 9.915393e-01  306.0307 2.943645e-04 -17.63768
183552   rs1906613       4_72 8.460684e-03  295.1509 2.422479e-06  17.34783
28533   rs34515871       1_76 1.000000e+00  299.2842 2.903316e-04  17.00000
422177  rs60469668      10_66 1.000000e+00  172.1082 1.669599e-04  16.96471
183513  rs13111704       4_72 0.000000e+00  393.8288 0.000000e+00 -16.91743
28536   rs12058931       1_76 1.643627e-04  221.4339 3.530676e-08  16.54054
183537   rs4124159       4_72 0.000000e+00 1159.1777 0.000000e+00 -16.14414
183535   rs1906606       4_72 0.000000e+00 1154.4323 0.000000e+00 -16.02703
183542   rs4032974       4_72 0.000000e+00 1146.0129 0.000000e+00 -15.86726
739346  rs11773845       7_70 0.000000e+00  555.7324 0.000000e+00  15.73134
183562  rs17513625       4_72 0.000000e+00  176.0291 0.000000e+00  15.71130
183539  rs10006881       4_72 0.000000e+00 1136.3170 0.000000e+00 -15.68750
739366   rs1997571       7_70 0.000000e+00  553.9040 0.000000e+00  15.68657
739329   rs3807989       7_70 0.000000e+00  551.5677 0.000000e+00  15.62687
422173  rs12572965      10_66 2.273800e-04  132.7631 2.928466e-08  15.27059
422172  rs56965730      10_66 2.300323e-04  131.7918 2.940950e-08  15.21176
739367   rs1997572       7_70 0.000000e+00  525.1335 0.000000e+00  14.50685
739223   rs2270188       7_70 0.000000e+00  711.6787 0.000000e+00  14.43939
739224   rs2270189       7_70 0.000000e+00  711.6066 0.000000e+00  14.43939
739272   rs2109514       7_70 0.000000e+00  687.7806 0.000000e+00  14.42424
739276  rs55883210       7_70 0.000000e+00  687.5295 0.000000e+00  14.42424
739239  rs10271007       7_70 0.000000e+00  707.4203 0.000000e+00  14.39394
739256   rs6466579       7_70 0.000000e+00  703.1231 0.000000e+00  14.33333
739267   rs7795510       7_70 0.000000e+00  696.7170 0.000000e+00  14.22727
183507  rs12642151       4_72 2.282285e-12  237.6299 5.261159e-16  14.14286
422175   rs7094488      10_66 9.999992e-01  110.0143 1.067233e-04  14.10606

Gene set enrichment for genes with PIP>0.8

#GO enrichment analysis
library(enrichR)
Welcome to enrichR
Checking connection ... 
Enrichr ... Connection is Live!
FlyEnrichr ... Connection is available!
WormEnrichr ... Connection is available!
YeastEnrichr ... Connection is available!
FishEnrichr ... Connection is available!
dbs <- c("GO_Biological_Process_2021", "GO_Cellular_Component_2021", "GO_Molecular_Function_2021")
genes <- ctwas_gene_res$genename[ctwas_gene_res$susie_pip>0.8]

#number of genes for gene set enrichment
length(genes)
[1] 32
if (length(genes)>0){
  GO_enrichment <- enrichr(genes, dbs)

  for (db in dbs){
    print(db)
    df <- GO_enrichment[[db]]
    print(plotEnrich(GO_enrichment[[db]]))
    df <- df[df$Adjusted.P.value<0.05,c("Term", "Overlap", "Adjusted.P.value", "Genes")]
    print(df)
  }
  
  #DisGeNET enrichment
  
  # devtools::install_bitbucket("ibi_group/disgenet2r")
  library(disgenet2r)
  
  disgenet_api_key <- get_disgenet_api_key(
                    email = "wesleycrouse@gmail.com", 
                    password = "uchicago1" )
  
  Sys.setenv(DISGENET_API_KEY= disgenet_api_key)
  
  res_enrich <-disease_enrichment(entities=genes, vocabulary = "HGNC",
                               database = "CURATED" )
  
  df <- res_enrich@qresult[1:10, c("Description", "FDR", "Ratio",  "BgRatio")]
  print(df)
  
  #WebGestalt enrichment
  library(WebGestaltR)
  
  background <- ctwas_gene_res$genename
  
  #listGeneSet()
  databases <- c("pathway_KEGG", "disease_GLAD4U", "disease_OMIM")
  
  enrichResult <- WebGestaltR(enrichMethod="ORA", organism="hsapiens",
                              interestGene=genes, referenceGene=background,
                              enrichDatabase=databases, interestGeneType="genesymbol",
                              referenceGeneType="genesymbol", isOutput=F)
  print(enrichResult[,c("description", "size", "overlap", "FDR", "database", "userId")])
}
Uploading data to Enrichr... Done.
  Querying GO_Biological_Process_2021... Done.
  Querying GO_Cellular_Component_2021... Done.
  Querying GO_Molecular_Function_2021... Done.
Parsing results... Done.
[1] "GO_Biological_Process_2021"

Version Author Date
77de9fb sq-96 2021-12-22
[1] Term             Overlap          Adjusted.P.value Genes           
<0 rows> (or 0-length row.names)
[1] "GO_Cellular_Component_2021"

Version Author Date
77de9fb sq-96 2021-12-22
[1] Term             Overlap          Adjusted.P.value Genes           
<0 rows> (or 0-length row.names)
[1] "GO_Molecular_Function_2021"
[1] Term             Overlap          Adjusted.P.value Genes           
<0 rows> (or 0-length row.names)
MURC gene(s) from the input list not found in DisGeNET CURATEDJAM2 gene(s) from the input list not found in DisGeNET CURATEDAGAP5 gene(s) from the input list not found in DisGeNET CURATEDCMTM5 gene(s) from the input list not found in DisGeNET CURATEDDLEU1 gene(s) from the input list not found in DisGeNET CURATEDSEC23IP gene(s) from the input list not found in DisGeNET CURATEDMARS gene(s) from the input list not found in DisGeNET CURATEDPOPDC3 gene(s) from the input list not found in DisGeNET CURATEDAES gene(s) from the input list not found in DisGeNET CURATEDRP5-890E16.5 gene(s) from the input list not found in DisGeNET CURATEDSTK11IP gene(s) from the input list not found in DisGeNET CURATEDSP100 gene(s) from the input list not found in DisGeNET CURATEDRP11-325L7.2 gene(s) from the input list not found in DisGeNET CURATEDLINC01629 gene(s) from the input list not found in DisGeNET CURATEDBOK gene(s) from the input list not found in DisGeNET CURATED
                                                                   Description
5                                                          Atrial Fibrillation
85                                              Paroxysmal atrial fibrillation
156                                             Persistent atrial fibrillation
171                                               familial atrial fibrillation
37                                                                Cardiomegaly
130                                                        Cardiac Hypertrophy
58                                                    Congenital retrognathism
157                      HYPOTHYROIDISM, CONGENITAL, NONGOITROUS, 5 (disorder)
158                              Lipodystrophy, Congenital Generalized, Type 3
167 ATRIAL SEPTAL DEFECT 7 WITH OR WITHOUT ATRIOVENTRICULAR CONDUCTION DEFECTS
             FDR Ratio  BgRatio
5   7.915974e-11  9/17 160/9703
85  7.915974e-11  9/17 156/9703
156 7.915974e-11  9/17 156/9703
171 7.915974e-11  9/17 156/9703
37  1.217805e-02  3/17  82/9703
130 1.217805e-02  3/17  82/9703
58  1.467481e-02  1/17   1/9703
157 1.467481e-02  1/17   1/9703
158 1.467481e-02  1/17   1/9703
167 1.467481e-02  1/17   1/9703
******************************************

*                                        *

*          Welcome to WebGestaltR !      *

*                                        *

******************************************

Version Author Date
77de9fb sq-96 2021-12-22
Loading the functional categories...
Loading the ID list...
Loading the reference list...
Performing the enrichment analysis...
                 description size overlap         FDR       database
1            Isaacs Syndrome   55       5 0.001525439 disease_GLAD4U
2        Atrial Fibrillation   52       4 0.022003754 disease_GLAD4U
3 Atrioventricular block NOS   22       3 0.031360061 disease_GLAD4U
                          userId
1 SCN10A;GNB4;NKX2-5;CCND2;KLF12
2       PRRX1;SCN10A;NKX2-5;CAV1
3           SCN10A;NKX2-5;POPDC3

Sensitivity, specificity and precision for silver standard genes

library("readxl")

known_annotations <- read_xlsx("data/summary_known_genes_annotations.xlsx", sheet="LDL")
New names:
* `` -> ...4
* `` -> ...5
known_annotations <- unique(known_annotations$`Gene Symbol`)

unrelated_genes <- ctwas_gene_res$genename[!(ctwas_gene_res$genename %in% known_annotations)]

#number of genes in known annotations
print(length(known_annotations))
[1] 69
#number of genes in known annotations with imputed expression
print(sum(known_annotations %in% ctwas_gene_res$genename))
[1] 44
#assign ctwas, TWAS, and bystander genes
ctwas_genes <- ctwas_gene_res$genename[ctwas_gene_res$susie_pip>0.8]
twas_genes <- ctwas_gene_res$genename[abs(ctwas_gene_res$z)>sig_thresh]
novel_genes <- ctwas_genes[!(ctwas_genes %in% twas_genes)]

#significance threshold for TWAS
print(sig_thresh)
[1] 4.582104
#number of ctwas genes
length(ctwas_genes)
[1] 32
#number of TWAS genes
length(twas_genes)
[1] 94
#show novel genes (ctwas genes with not in TWAS genes)
ctwas_gene_res[ctwas_gene_res$genename %in% novel_genes,report_cols]
     genename region_tag susie_pip      mu2          PVE         z num_eqtl
6114  STK11IP      2_130 0.8406381 18.96825 1.546845e-05 -3.868022        2
712     SP100      2_135 0.8658723 18.73540 1.573719e-05 -3.671335        2
9691      BOK      2_144 0.8255975 19.18259 1.536335e-05  3.910125        3
9012    MTSS1       8_82 0.8594108 20.87861 1.740655e-05  4.402634        2
2444  SEC23IP      10_74 0.9517809 22.38862 2.067163e-05 -4.565228        2
8420     MARS      12_36 0.8180336 17.81958 1.414097e-05 -3.366197        1
2138      AES       19_4 0.9403219 20.20438 1.843030e-05  4.182804        3
6914     JAM2       21_9 0.9640543 22.19279 2.075505e-05  4.563232        2
#sensitivity / recall
sensitivity <- rep(NA,2)
names(sensitivity) <- c("ctwas", "TWAS")
sensitivity["ctwas"] <- sum(ctwas_genes %in% known_annotations)/length(known_annotations)
sensitivity["TWAS"] <- sum(twas_genes %in% known_annotations)/length(known_annotations)
sensitivity
ctwas  TWAS 
    0     0 
#specificity
specificity <- rep(NA,2)
names(specificity) <- c("ctwas", "TWAS")
specificity["ctwas"] <- sum(!(unrelated_genes %in% ctwas_genes))/length(unrelated_genes)
specificity["TWAS"] <- sum(!(unrelated_genes %in% twas_genes))/length(unrelated_genes)
specificity
    ctwas      TWAS 
0.9970420 0.9913108 
#precision / PPV
precision <- rep(NA,2)
names(precision) <- c("ctwas", "TWAS")
precision["ctwas"] <- sum(ctwas_genes %in% known_annotations)/length(ctwas_genes)
precision["TWAS"] <- sum(twas_genes %in% known_annotations)/length(twas_genes)
precision
ctwas  TWAS 
    0     0 
#ROC curves

pip_range <- (0:1000)/1000
sensitivity <- rep(NA, length(pip_range))
specificity <- rep(NA, length(pip_range))

for (index in 1:length(pip_range)){
  pip <- pip_range[index]
  ctwas_genes <- ctwas_gene_res$genename[ctwas_gene_res$susie_pip>=pip]
  sensitivity[index] <- sum(ctwas_genes %in% known_annotations)/length(known_annotations)
  specificity[index] <- sum(!(unrelated_genes %in% ctwas_genes))/length(unrelated_genes)
}

plot(1-specificity, sensitivity, type="l", xlim=c(0,1), ylim=c(0,1))

sig_thresh_range <- seq(from=0, to=max(abs(ctwas_gene_res$z)), length.out=length(pip_range))

for (index in 1:length(sig_thresh_range)){
  sig_thresh_plot <- sig_thresh_range[index]
  twas_genes <- ctwas_gene_res$genename[abs(ctwas_gene_res$z)>=sig_thresh_plot]
  sensitivity[index] <- sum(twas_genes %in% known_annotations)/length(known_annotations)
  specificity[index] <- sum(!(unrelated_genes %in% twas_genes))/length(unrelated_genes)
}

lines(1-specificity, sensitivity, xlim=c(0,1), ylim=c(0,1), col="red", lty=2)

Version Author Date
9cf0e70 sq-96 2021-12-22

Sensitivity, specificity and precision for silver standard genes - bystanders only

This section first uses imputed silver standard genes to identify bystander genes within 1Mb. The bystander gene list is then subset to only genes with imputed expression in this analysis. Then, the ctwas and TWAS gene lists from this analysis are subset to only genes that are in the (subset) silver standard and bystander genes. These gene lists are then used to compute sensitivity, specificity and precision for ctwas and TWAS.

library(biomaRt)
library(GenomicRanges)
Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'
The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB
The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs
The following objects are masked from 'package:base':

    anyDuplicated, append, as.data.frame, basename, cbind, colnames,
    dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
    grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
    rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min
Loading required package: S4Vectors

Attaching package: 'S4Vectors'
The following object is masked from 'package:base':

    expand.grid
Loading required package: IRanges
Loading required package: GenomeInfoDb
ensembl <- useEnsembl(biomart="ENSEMBL_MART_ENSEMBL", dataset="hsapiens_gene_ensembl")
G_list <- getBM(filters= "chromosome_name", attributes= c("hgnc_symbol","chromosome_name","start_position","end_position","gene_biotype"), values=1:22, mart=ensembl)
G_list <- G_list[G_list$hgnc_symbol!="",]
G_list <- G_list[G_list$gene_biotype %in% c("protein_coding","lncRNA"),]
G_list$start <- G_list$start_position
G_list$end <- G_list$end_position
G_list_granges <- makeGRangesFromDataFrame(G_list, keep.extra.columns=T)

#remove genes without imputed expression from gene lists
known_annotations <- known_annotations[known_annotations %in% ctwas_gene_res$genename]

known_annotations_positions <- G_list[G_list$hgnc_symbol %in% known_annotations,]
half_window <- 1000000
known_annotations_positions$start <- known_annotations_positions$start_position - half_window
known_annotations_positions$end <- known_annotations_positions$end_position + half_window
known_annotations_positions$start[known_annotations_positions$start<1] <- 1
known_annotations_granges <- makeGRangesFromDataFrame(known_annotations_positions, keep.extra.columns=T)

bystanders <- findOverlaps(known_annotations_granges,G_list_granges)
bystanders <- unique(subjectHits(bystanders))
bystanders <- G_list$hgnc_symbol[bystanders]
bystanders <- unique(bystanders[!(bystanders %in% known_annotations)])
unrelated_genes <- bystanders

#save gene lists
save(known_annotations, file=paste0(results_dir, "/known_annotations.Rd"))
save(unrelated_genes, file=paste0(results_dir, "/bystanders.Rd"))


load(paste0(results_dir, "/known_annotations.Rd"))
load(paste0(results_dir, "/bystanders.Rd"))


#remove genes without imputed expression from bystander list
unrelated_genes <- unrelated_genes[unrelated_genes %in% ctwas_gene_res$genename]

#number of genes in known annotations (with imputed expression)
print(length(known_annotations))
[1] 44
#number of bystander genes (with imputed expression)
print(length(unrelated_genes))
[1] 605
#subset results to genes in known annotations or bystanders
ctwas_gene_res_subset <- ctwas_gene_res[ctwas_gene_res$genename %in% c(known_annotations, unrelated_genes),]

#assign ctwas and TWAS genes
ctwas_genes <- ctwas_gene_res_subset$genename[ctwas_gene_res_subset$susie_pip>0.8]
twas_genes <- ctwas_gene_res_subset$genename[abs(ctwas_gene_res_subset$z)>sig_thresh]

#significance threshold for TWAS
print(sig_thresh)
[1] 4.582104
#number of ctwas genes (in known annotations or bystanders)
length(ctwas_genes)
[1] 1
#number of TWAS genes (in known annotations or bystanders)
length(twas_genes)
[1] 2
#sensitivity / recall
sensitivity <- rep(NA,2)
names(sensitivity) <- c("ctwas", "TWAS")
sensitivity["ctwas"] <- sum(ctwas_genes %in% known_annotations)/length(known_annotations)
sensitivity["TWAS"] <- sum(twas_genes %in% known_annotations)/length(known_annotations)
sensitivity
ctwas  TWAS 
    0     0 
#specificity / (1 - False Positive Rate)
specificity <- rep(NA,2)
names(specificity) <- c("ctwas", "TWAS")
specificity["ctwas"] <- sum(!(unrelated_genes %in% ctwas_genes))/length(unrelated_genes)
specificity["TWAS"] <- sum(!(unrelated_genes %in% twas_genes))/length(unrelated_genes)
specificity
    ctwas      TWAS 
0.9983471 0.9966942 
#precision / PPV / (1 - False Discovery Rate)
precision <- rep(NA,2)
names(precision) <- c("ctwas", "TWAS")
precision["ctwas"] <- sum(ctwas_genes %in% known_annotations)/length(ctwas_genes)
precision["TWAS"] <- sum(twas_genes %in% known_annotations)/length(twas_genes)
precision
ctwas  TWAS 
    0     0 
#store sensitivity and specificity calculations for plots
sensitivity_plot <- sensitivity
specificity_plot <- specificity

#precision / PPV by PIP bin
pip_range <- c(0.2, 0.4, 0.6, 0.8, 1)
precision_range <- rep(NA, length(pip_range))

for (i in 1:length(pip_range)){
  pip_upper <- pip_range[i]

  if (i==1){
    pip_lower <- 0
  } else {
    pip_lower <- pip_range[i-1]
  }
  
  #assign ctwas genes in PIP bin
  ctwas_genes <- ctwas_gene_res_subset$genename[ctwas_gene_res_subset$susie_pip>=pip_lower & ctwas_gene_res_subset$susie_pip<pip_upper]
  
  precision_range[i] <- sum(ctwas_genes %in% known_annotations)/length(ctwas_genes)
}

names(precision_range) <- paste(c(0, pip_range[-length(pip_range)]), pip_range,sep=" - ")

barplot(precision_range, ylim=c(0,1), main="Precision by PIP Range", xlab="PIP Range", ylab="Precision")
abline(h=0.2, lty=2)
abline(h=0.4, lty=2)
abline(h=0.6, lty=2)
abline(h=0.8, lty=2)
barplot(precision_range, add=T, col="darkgrey")

#precision / PPV by PIP threshold
#pip_range <- c(0.2, 0.4, 0.6, 0.8, 1)
pip_range <- c(0.5, 0.8, 1)
precision_range <- rep(NA, length(pip_range))
number_detected <- rep(NA, length(pip_range))

for (i in 1:length(pip_range)){
  pip_upper <- pip_range[i]

  if (i==1){
    pip_lower <- 0
  } else {
    pip_lower <- pip_range[i-1]
  }
  
  #assign ctwas genes using PIP threshold
  ctwas_genes <- ctwas_gene_res_subset$genename[ctwas_gene_res_subset$susie_pip>=pip_lower]
  
  number_detected[i] <- length(ctwas_genes)
  precision_range[i] <- sum(ctwas_genes %in% known_annotations)/length(ctwas_genes)
}

names(precision_range) <- paste0(">= ", c(0, pip_range[-length(pip_range)]))

precision_range <- precision_range*100

precision_range <- c(precision_range, precision["TWAS"]*100)
names(precision_range)[4] <- "TWAS Bonferroni"
number_detected <- c(number_detected, length(twas_genes))

barplot(precision_range, ylim=c(0,100), main="Precision for Distinguishing Silver Standard and Bystander Genes", xlab="PIP Threshold for Detection", ylab="% of Detected Genes in Silver Standard")
abline(h=20, lty=2)
abline(h=40, lty=2)
abline(h=60, lty=2)
abline(h=80, lty=2)
xx <- barplot(precision_range, add=T, col=c(rep("darkgrey",3), "white"))
text(x = xx, y = rep(0, length(number_detected)), label = paste0(number_detected, " detected"), pos = 3, cex=0.8)

#text(x = xx, y = precision_range, label = paste0(round(precision_range,1), "%"), pos = 3, cex=0.8, offset = 1.5)

#false discovery rate by PIP threshold

barplot(100-precision_range, ylim=c(0,100), main="False Discovery Rate for Distinguishing Silver Standard and Bystander Genes", xlab="PIP Threshold for Detection", ylab="% Bystanders in Detected Genes")
abline(h=20, lty=2)
abline(h=40, lty=2)
abline(h=60, lty=2)
abline(h=80, lty=2)
xx <- barplot(100-precision_range, add=T, col=c(rep("darkgrey",3), "white"))
text(x = xx, y = rep(0, length(number_detected)), label = paste0(number_detected, " detected"), pos = 3, cex=0.8)

#text(x = xx, y = precision_range, label = paste0(round(precision_range,1), "%"), pos = 3, cex=0.8, offset = 1.5)

#ROC curves

pip_range <- (0:1000)/1000
sensitivity <- rep(NA, length(pip_range))
specificity <- rep(NA, length(pip_range))

for (index in 1:length(pip_range)){
  pip <- pip_range[index]
  ctwas_genes <- ctwas_gene_res_subset$genename[ctwas_gene_res_subset$susie_pip>=pip]
  sensitivity[index] <- sum(ctwas_genes %in% known_annotations)/length(known_annotations)
  specificity[index] <- sum(!(unrelated_genes %in% ctwas_genes))/length(unrelated_genes)
}

plot(1-specificity, sensitivity, type="l", xlim=c(0,1), ylim=c(0,1), main="", xlab="1 - Specificity", ylab="Sensitivity")
title(expression("ROC Curve for cTWAS (black) and TWAS (" * phantom("red") * ")"))
title(expression(phantom("ROC Curve for cTWAS (black) and TWAS (") * "red" * phantom(")")), col.main="red")

sig_thresh_range <- seq(from=0, to=max(abs(ctwas_gene_res_subset$z)), length.out=length(pip_range))

for (index in 1:length(sig_thresh_range)){
  sig_thresh_plot <- sig_thresh_range[index]
  twas_genes <- ctwas_gene_res_subset$genename[abs(ctwas_gene_res_subset$z)>=sig_thresh_plot]
  sensitivity[index] <- sum(twas_genes %in% known_annotations)/length(known_annotations)
  specificity[index] <- sum(!(unrelated_genes %in% twas_genes))/length(unrelated_genes)
}

lines(1-specificity, sensitivity, xlim=c(0,1), ylim=c(0,1), col="red", lty=1)

abline(a=0,b=1,lty=3)

#add previously computed points from the analysis
ctwas_genes <- ctwas_gene_res_subset$genename[ctwas_gene_res_subset$susie_pip>0.8]
twas_genes <- ctwas_gene_res_subset$genename[abs(ctwas_gene_res_subset$z)>sig_thresh]

points(1-specificity_plot["ctwas"], sensitivity_plot["ctwas"], pch=21, bg="black")
points(1-specificity_plot["TWAS"], sensitivity_plot["TWAS"], pch=21, bg="red")


sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] GenomicRanges_1.36.1 GenomeInfoDb_1.20.0  IRanges_2.18.1      
 [4] S4Vectors_0.22.1     BiocGenerics_0.30.0  biomaRt_2.40.1      
 [7] readxl_1.3.1         WebGestaltR_0.4.4    disgenet2r_0.99.2   
[10] enrichR_3.0          cowplot_1.0.0        ggplot2_3.3.5       
[13] workflowr_1.6.2     

loaded via a namespace (and not attached):
 [1] bitops_1.0-7           fs_1.5.2               bit64_4.0.5           
 [4] progress_1.2.2         doParallel_1.0.16      httr_1.4.2            
 [7] rprojroot_2.0.2        tools_3.6.1            doRNG_1.8.2           
[10] utf8_1.2.2             R6_2.5.1               DBI_1.1.1             
[13] colorspace_2.0-2       withr_2.4.3            prettyunits_1.1.1     
[16] tidyselect_1.1.1       bit_4.0.4              curl_4.3.2            
[19] compiler_3.6.1         git2r_0.26.1           Biobase_2.44.0        
[22] labeling_0.4.2         scales_1.1.1           readr_2.1.1           
[25] stringr_1.4.0          apcluster_1.4.8        digest_0.6.29         
[28] rmarkdown_2.11         svglite_1.2.2          XVector_0.24.0        
[31] pkgconfig_2.0.3        htmltools_0.5.2        fastmap_1.1.0         
[34] highr_0.9              rlang_0.4.12           RSQLite_2.2.8         
[37] jquerylib_0.1.4        farver_2.1.0           generics_0.1.1        
[40] jsonlite_1.7.2         vroom_1.5.7            dplyr_1.0.7           
[43] RCurl_1.98-1.5         magrittr_2.0.1         GenomeInfoDbData_1.2.1
[46] Matrix_1.2-18          Rcpp_1.0.7             munsell_0.5.0         
[49] fansi_0.5.0            gdtools_0.1.9          lifecycle_1.0.1       
[52] stringi_1.7.6          whisker_0.3-2          yaml_2.2.1            
[55] zlibbioc_1.30.0        plyr_1.8.6             grid_3.6.1            
[58] blob_1.2.2             promises_1.0.1         crayon_1.4.2          
[61] lattice_0.20-38        hms_1.1.1              knitr_1.36            
[64] pillar_1.6.4           igraph_1.2.10          rjson_0.2.20          
[67] rngtools_1.5.2         reshape2_1.4.4         codetools_0.2-16      
[70] XML_3.99-0.3           glue_1.5.1             evaluate_0.14         
[73] data.table_1.14.2      vctrs_0.3.8            tzdb_0.2.0            
[76] httpuv_1.5.1           foreach_1.5.1          cellranger_1.1.0      
[79] gtable_0.3.0           purrr_0.3.4            assertthat_0.2.1      
[82] cachem_1.0.6           xfun_0.29              later_0.8.0           
[85] tibble_3.1.6           iterators_1.0.13       AnnotationDbi_1.46.0  
[88] memoise_2.0.1          ellipsis_0.3.2