Atrial fibrillation - Heart Atrial Appendage
sheng Qian
2021-12-18
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.6548518Load 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)
#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.12523867Genes with highest PIPs
#distribution of PIPs
hist(ctwas_gene_res$susie_pip, xlim=c(0,1), main="Distribution of Gene PIPs")
#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 1Genes 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")
#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 2Genes 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 1Genes 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 1Comparing 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)
#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)
#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 1Locus 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
[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.8292683SNPs 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.931818SNPs 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.9647059SNPs 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.10606Gene 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] 32if (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;POPDC3Sensitivity, specificity and precision for silver standard genes
library("readxl")
known_annotations <- read_xlsx("data/summary_known_genes_annotations.xlsx", sheet="LDL")New names:
* `` -> ...4
* `` -> ...5known_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)
sensitivityctwas 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)
precisionctwas 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: stats4Loading required package: BiocGenericsLoading 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, parSapplyLBThe following objects are masked from 'package:stats':
IQR, mad, sd, var, xtabsThe 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.minLoading required package: S4Vectors
Attaching package: 'S4Vectors'The following object is masked from 'package:base':
expand.gridLoading required package: IRangesLoading required package: GenomeInfoDbensembl <- 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)
sensitivityctwas 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)
precisionctwas 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