Cholesterol (quantile) - Liver
wesleycrouse
2021-09-09
Last updated: 2021-09-09
Checks: 6 1
Knit directory: ctwas_applied/
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.
The R Markdown file has unstaged changes. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.
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(20210726) 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.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
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 59e5f4d. 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:
Unstaged changes:
Modified: analysis/ukb-d-30500_irnt_Liver.Rmd
Modified: analysis/ukb-d-30600_irnt_Liver.Rmd
Modified: analysis/ukb-d-30610_irnt_Liver.Rmd
Modified: analysis/ukb-d-30620_irnt_Liver.Rmd
Modified: analysis/ukb-d-30630_irnt_Liver.Rmd
Modified: analysis/ukb-d-30640_irnt_Liver.Rmd
Modified: analysis/ukb-d-30650_irnt_Liver.Rmd
Modified: analysis/ukb-d-30660_irnt_Liver.Rmd
Modified: analysis/ukb-d-30670_irnt_Liver.Rmd
Modified: analysis/ukb-d-30680_irnt_Liver.Rmd
Modified: analysis/ukb-d-30690_irnt_Liver.Rmd
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/ukb-d-30690_irnt_Liver.Rmd) and HTML (docs/ukb-d-30690_irnt_Liver.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 | cbf7408 | wesleycrouse | 2021-09-08 | adding enrichment to reports |
| html | cbf7408 | wesleycrouse | 2021-09-08 | adding enrichment to reports |
| Rmd | 4970e3e | wesleycrouse | 2021-09-08 | updating reports |
| html | 4970e3e | wesleycrouse | 2021-09-08 | updating reports |
| Rmd | 627a4e1 | wesleycrouse | 2021-09-07 | adding heritability |
| Rmd | dfd2b5f | wesleycrouse | 2021-09-07 | regenerating reports |
| html | dfd2b5f | wesleycrouse | 2021-09-07 | regenerating reports |
| Rmd | 61b53b3 | wesleycrouse | 2021-09-06 | updated PVE calculation |
| html | 61b53b3 | wesleycrouse | 2021-09-06 | updated PVE calculation |
| Rmd | 837dd01 | wesleycrouse | 2021-09-01 | adding additional fixedsigma report |
| Rmd | d0a5417 | wesleycrouse | 2021-08-30 | adding new reports to the index |
| Rmd | 0922de7 | wesleycrouse | 2021-08-18 | updating all reports with locus plots |
| html | 0922de7 | wesleycrouse | 2021-08-18 | updating all reports with locus plots |
| html | 1c62980 | wesleycrouse | 2021-08-11 | Updating reports |
| Rmd | 5981e80 | wesleycrouse | 2021-08-11 | Adding more reports |
| html | 5981e80 | wesleycrouse | 2021-08-11 | Adding more reports |
| Rmd | 05a98b7 | wesleycrouse | 2021-08-07 | adding additional results |
| html | 05a98b7 | wesleycrouse | 2021-08-07 | adding additional results |
| html | 03e541c | wesleycrouse | 2021-07-29 | Cleaning up report generation |
| Rmd | 276893d | wesleycrouse | 2021-07-29 | Updating reports |
| html | 276893d | wesleycrouse | 2021-07-29 | Updating reports |
Overview
These are the results of a ctwas analysis of the UK Biobank trait Cholesterol (quantile) using Liver gene weights.
The GWAS was conducted by the Neale Lab, and the biomarker traits we analyzed are discussed here. Summary statistics were obtained from IEU OpenGWAS using GWAS ID: ukb-d-30690_irnt. Results were obtained from from IEU rather than Neale Lab because they are in a standardard format (GWAS VCF). Note that 3 of the 34 biomarker traits were not available from IEU and were excluded from analysis.
The weights are mashr GTEx v8 models on Liver eQTL obtained from PredictDB. We performed a full harmonization of the variants, including recovering strand ambiguous variants. This procedure is discussed in a separate document. (TO-DO: add report that describes harmonization)
LD matrices were computed from a 10% subset of Neale lab subjects. Subjects were matched using the plate and well information from genotyping. We included only biallelic variants with MAF>0.01 in the original Neale Lab GWAS. (TO-DO: add more details [number of subjects, variants, etc])
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] 10901#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
1070 768 652 417 494 611 548 408 405 434 634 629 195 365 354
16 17 18 19 20 21 22
526 663 160 859 306 114 289 #proportion of imputed weights without missing variants
mean(qclist_all$nmiss==0)[1] 0.8366205Load ctwas results
#load ctwas results
ctwas_res <- data.table::fread(paste0(results_dir, "/", analysis_id, "_ctwas.susieIrss.txt"))
#make unique identifier for regions
ctwas_res$region_tag <- paste(ctwas_res$region_tag1, ctwas_res$region_tag2, sep="_")
#compute PVE for each gene/SNP
ctwas_res$PVE = ctwas_res$susie_pip*ctwas_res$mu2/sample_size #check PVE calculation
#separate gene and SNP results
ctwas_gene_res <- ctwas_res[ctwas_res$type == "gene", ]
ctwas_gene_res <- data.frame(ctwas_gene_res)
ctwas_snp_res <- ctwas_res[ctwas_res$type == "SNP", ]
ctwas_snp_res <- data.frame(ctwas_snp_res)
#add gene information to results
sqlite <- RSQLite::dbDriver("SQLite")
db = RSQLite::dbConnect(sqlite, paste0("/project2/compbio/predictdb/mashr_models/mashr_", weight, ".db"))
query <- function(...) RSQLite::dbGetQuery(db, ...)
gene_info <- query("select gene, genename from extra")
gene_info <- query("select gene, genename, gene_type from extra")
RSQLite::dbDisconnect(db)
ctwas_gene_res <- cbind(ctwas_gene_res, gene_info[sapply(ctwas_gene_res$id, match, gene_info$gene), c("genename", "gene_type")])
#add z score to results
load(paste0(results_dir, "/", analysis_id, "_expr_z_gene.Rd"))
ctwas_gene_res$z <- z_gene[ctwas_gene_res$id,]$z
#load(paste0(results_dir, "/", analysis_id, "_expr_z_snp.Rd")) #for new version, stored after harmonization
z_snp <- readRDS(paste0(results_dir, "/", trait_id, ".RDS")) #for old version, unharmonized
z_snp <- z_snp[z_snp$id %in% ctwas_snp_res$id,] #subset snps to those included in analysis, note some are duplicated, need to match which allele was used
ctwas_snp_res$z <- z_snp$z[match(ctwas_snp_res$id, z_snp$id)] #for duplicated snps, this arbitrarily uses the first allele
ctwas_snp_res$z_flag <- as.numeric(ctwas_snp_res$id %in% z_snp$id[duplicated(z_snp$id)]) #mark the unclear z scores, flag=1
#formatting and rounding for tables
ctwas_gene_res$z <- round(ctwas_gene_res$z,2)
ctwas_snp_res$z <- round(ctwas_snp_res$z,2)
ctwas_gene_res$susie_pip <- round(ctwas_gene_res$susie_pip,3)
ctwas_snp_res$susie_pip <- round(ctwas_snp_res$susie_pip,3)
ctwas_gene_res$mu2 <- round(ctwas_gene_res$mu2,2)
ctwas_snp_res$mu2 <- round(ctwas_snp_res$mu2,2)
ctwas_gene_res$PVE <- signif(ctwas_gene_res$PVE, 2)
ctwas_snp_res$PVE <- signif(ctwas_snp_res$PVE, 2)
#merge gene and snp results with added information
ctwas_gene_res$z_flag=NA
ctwas_snp_res$genename=NA
ctwas_snp_res$gene_type=NA
ctwas_res <- rbind(ctwas_gene_res,
ctwas_snp_res[,colnames(ctwas_gene_res)])
#store columns to report
report_cols <- colnames(ctwas_gene_res)[!(colnames(ctwas_gene_res) %in% c("type", "region_tag1", "region_tag2", "cs_index", "gene_type", "z_flag", "id", "chrom", "pos"))]
first_cols <- c("genename", "region_tag")
report_cols <- c(first_cols, report_cols[!(report_cols %in% first_cols)])
report_cols_snps <- c("id", report_cols[-1])
#get number of SNPs from s1 results; adjust for thin
ctwas_res_s1 <- data.table::fread(paste0(results_dir, "/", analysis_id, "_ctwas.s1.susieIrss.txt"))
n_snps <- sum(ctwas_res_s1$type=="SNP")/thin
rm(ctwas_res_s1)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 |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
#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.0084575074 0.0001855158 #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
59.37541 11.15032 #report sample size
print(sample_size)[1] 344278#report group size
group_size <- c(nrow(ctwas_gene_res), n_snps)
print(group_size)[1] 10901 8697330#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.01590033 0.05225709 #compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))[1] 0.06108648 0.48920986Genes with highest PIPs
#distribution of PIPs
hist(ctwas_gene_res$susie_pip, xlim=c(0,1), main="Distribution of Gene PIPs")
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
#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
4435 PSRC1 1_67 1.000 1162.99 3.4e-03 -34.57
5544 CNIH4 1_114 1.000 62.85 1.8e-04 7.79
12687 RP4-781K5.7 1_121 1.000 257.07 7.5e-04 -16.29
5563 ABCG8 2_27 1.000 280.94 8.2e-04 -18.98
10657 TRIM39 6_24 1.000 78.91 2.3e-04 9.46
6391 TTC39B 9_13 1.000 96.14 2.8e-04 -9.68
7410 ABCA1 9_53 1.000 240.92 7.0e-04 15.14
5991 FADS1 11_34 1.000 306.16 8.9e-04 17.68
12008 HPR 16_38 1.000 195.86 5.7e-04 -16.51
5389 RPS11 19_34 1.000 12551.14 3.6e-02 -3.91
1999 PRKD2 19_33 0.999 36.64 1.1e-04 5.71
3721 INSIG2 2_69 0.998 176.69 5.1e-04 -8.90
8531 TNKS 8_12 0.992 87.81 2.5e-04 14.53
10625 MSH5 6_26 0.990 82.43 2.4e-04 10.14
1144 ASAP3 1_16 0.989 59.12 1.7e-04 7.73
9985 LITAF 16_12 0.980 29.57 8.4e-05 4.93
6220 PELO 5_31 0.976 47.95 1.4e-04 6.72
3212 CCND2 12_4 0.970 31.24 8.8e-05 -5.21
6996 ACP6 1_73 0.957 25.64 7.1e-05 4.62
9390 GAS6 13_62 0.949 86.25 2.4e-04 -9.89
8865 FUT2 19_33 0.947 72.34 2.0e-04 -12.26
3247 KDSR 18_35 0.945 25.34 7.0e-05 -4.57
7040 INHBB 2_70 0.943 53.31 1.5e-04 -7.10
6093 CSNK1G3 5_75 0.939 78.56 2.1e-04 8.69
4704 DDX56 7_32 0.939 43.46 1.2e-04 8.10
3562 ACVR1C 2_94 0.933 28.40 7.7e-05 -4.91
8931 CRACR2B 11_1 0.930 23.43 6.3e-05 -4.32
6778 PKN3 9_66 0.922 50.32 1.3e-04 -6.81
9062 KLHDC7A 1_13 0.895 23.71 6.2e-05 4.42
6100 ALLC 2_2 0.894 39.41 1.0e-04 6.01
2204 AKNA 9_59 0.890 27.71 7.2e-05 -4.77
3300 C10orf88 10_77 0.877 39.45 1.0e-04 -6.99
7924 PXK 3_40 0.859 35.85 8.9e-05 -4.56
10763 NYNRIN 14_3 0.851 41.98 1.0e-04 6.99
1114 SRRT 7_62 0.838 37.10 9.0e-05 5.95
9072 SPTY2D1 11_13 0.828 42.22 1.0e-04 -6.27Genes 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 |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
#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
5389 RPS11 19_34 1.000 12551.14 3.6e-02 -3.91
1227 FLT3LG 19_34 0.000 10848.93 0.0e+00 3.51
5393 RCN3 19_34 0.000 4067.51 0.0e+00 4.50
1931 FCGRT 19_34 0.000 3705.52 0.0e+00 3.50
3804 PRRG2 19_34 0.000 1833.53 0.0e+00 4.02
3803 PRMT1 19_34 0.000 1263.60 0.0e+00 3.25
3805 SCAF1 19_34 0.000 1253.48 0.0e+00 3.26
3802 IRF3 19_34 0.000 1229.05 0.0e+00 3.53
4435 PSRC1 1_67 1.000 1162.99 3.4e-03 -34.57
1940 SLC17A7 19_34 0.000 885.15 0.0e+00 2.58
5436 PSMA5 1_67 0.006 859.91 1.4e-05 -29.65
6957 USP1 1_39 0.771 484.88 1.1e-03 22.52
4317 ANGPTL3 1_39 0.238 482.73 3.3e-04 22.47
1932 PIH1D1 19_34 0.000 471.54 0.0e+00 -4.64
6476 SUPV3L1 10_46 0.000 466.70 6.6e-12 -2.21
12131 APOC4 19_32 0.000 337.87 0.0e+00 11.76
10291 CTC-301O7.4 19_34 0.000 315.81 0.0e+00 -2.19
5991 FADS1 11_34 1.000 306.16 8.9e-04 17.68
4050 TOMM40 19_32 0.000 286.06 0.0e+00 -1.44
5563 ABCG8 2_27 1.000 280.94 8.2e-04 -18.98Genes 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
5389 RPS11 19_34 1.000 12551.14 0.03600 -3.91
4435 PSRC1 1_67 1.000 1162.99 0.00340 -34.57
6957 USP1 1_39 0.771 484.88 0.00110 22.52
5991 FADS1 11_34 1.000 306.16 0.00089 17.68
5563 ABCG8 2_27 1.000 280.94 0.00082 -18.98
12687 RP4-781K5.7 1_121 1.000 257.07 0.00075 -16.29
7410 ABCA1 9_53 1.000 240.92 0.00070 15.14
12008 HPR 16_38 1.000 195.86 0.00057 -16.51
3721 INSIG2 2_69 0.998 176.69 0.00051 -8.90
4317 ANGPTL3 1_39 0.238 482.73 0.00033 22.47
6391 TTC39B 9_13 1.000 96.14 0.00028 -9.68
8531 TNKS 8_12 0.992 87.81 0.00025 14.53
10625 MSH5 6_26 0.990 82.43 0.00024 10.14
9390 GAS6 13_62 0.949 86.25 0.00024 -9.89
10657 TRIM39 6_24 1.000 78.91 0.00023 9.46
6093 CSNK1G3 5_75 0.939 78.56 0.00021 8.69
8865 FUT2 19_33 0.947 72.34 0.00020 -12.26
5544 CNIH4 1_114 1.000 62.85 0.00018 7.79
1144 ASAP3 1_16 0.989 59.12 0.00017 7.73
2092 SP4 7_19 0.592 84.35 0.00015 9.49Genes 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
4435 PSRC1 1_67 1.000 1162.99 3.4e-03 -34.57
5436 PSMA5 1_67 0.006 859.91 1.4e-05 -29.65
6957 USP1 1_39 0.771 484.88 1.1e-03 22.52
4317 ANGPTL3 1_39 0.238 482.73 3.3e-04 22.47
5563 ABCG8 2_27 1.000 280.94 8.2e-04 -18.98
5991 FADS1 11_34 1.000 306.16 8.9e-04 17.68
3441 POLK 5_45 0.003 215.49 2.1e-06 17.18
7547 LIPC 15_26 0.000 142.53 1.1e-09 -17.16
12008 HPR 16_38 1.000 195.86 5.7e-04 -16.51
6970 ATXN7L2 1_67 0.006 268.37 4.4e-06 -16.39
12687 RP4-781K5.7 1_121 1.000 257.07 7.5e-04 -16.29
4507 FADS2 11_34 0.004 260.01 2.7e-06 16.21
7955 FEN1 11_34 0.004 260.01 2.7e-06 16.21
7410 ABCA1 9_53 1.000 240.92 7.0e-04 15.14
5240 NLRC5 16_31 0.062 254.12 4.6e-05 -14.63
8531 TNKS 8_12 0.992 87.81 2.5e-04 14.53
9978 ANKDD1B 5_45 0.003 134.31 1.3e-06 14.36
11684 RP11-136O12.2 8_83 0.001 251.47 8.4e-07 13.69
2465 APOA5 11_70 0.021 160.13 1.0e-05 -13.56
8865 FUT2 19_33 0.947 72.34 2.0e-04 -12.26Comparing 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 |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
#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 |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
#proportion of significant z scores
mean(abs(ctwas_gene_res$z) > sig_thresh)[1] 0.0245849#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
4435 PSRC1 1_67 1.000 1162.99 3.4e-03 -34.57
5436 PSMA5 1_67 0.006 859.91 1.4e-05 -29.65
6957 USP1 1_39 0.771 484.88 1.1e-03 22.52
4317 ANGPTL3 1_39 0.238 482.73 3.3e-04 22.47
5563 ABCG8 2_27 1.000 280.94 8.2e-04 -18.98
5991 FADS1 11_34 1.000 306.16 8.9e-04 17.68
3441 POLK 5_45 0.003 215.49 2.1e-06 17.18
7547 LIPC 15_26 0.000 142.53 1.1e-09 -17.16
12008 HPR 16_38 1.000 195.86 5.7e-04 -16.51
6970 ATXN7L2 1_67 0.006 268.37 4.4e-06 -16.39
12687 RP4-781K5.7 1_121 1.000 257.07 7.5e-04 -16.29
4507 FADS2 11_34 0.004 260.01 2.7e-06 16.21
7955 FEN1 11_34 0.004 260.01 2.7e-06 16.21
7410 ABCA1 9_53 1.000 240.92 7.0e-04 15.14
5240 NLRC5 16_31 0.062 254.12 4.6e-05 -14.63
8531 TNKS 8_12 0.992 87.81 2.5e-04 14.53
9978 ANKDD1B 5_45 0.003 134.31 1.3e-06 14.36
11684 RP11-136O12.2 8_83 0.001 251.47 8.4e-07 13.69
2465 APOA5 11_70 0.021 160.13 1.0e-05 -13.56
8865 FUT2 19_33 0.947 72.34 2.0e-04 -12.26Locus plots for genes and SNPs
ctwas_gene_res_sortz <- ctwas_gene_res[order(-abs(ctwas_gene_res$z)),]
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])
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: 1_67"
genename region_tag susie_pip mu2 PVE z
4434 VAV3 1_67 0.068 28.25 5.6e-06 -2.36
1073 SLC25A24 1_67 0.006 5.13 8.3e-08 0.45
6966 FAM102B 1_67 0.006 5.81 1.0e-07 -0.59
3009 STXBP3 1_67 0.006 8.91 1.6e-07 1.92
3438 GPSM2 1_67 0.006 8.28 1.4e-07 -1.66
3437 CLCC1 1_67 0.006 10.66 1.7e-07 2.58
10286 TAF13 1_67 0.009 10.37 2.6e-07 -1.65
10955 TMEM167B 1_67 0.012 13.26 4.6e-07 -1.75
315 SARS 1_67 0.008 67.26 1.6e-06 7.85
5431 SYPL2 1_67 0.007 132.95 2.8e-06 -11.44
5436 PSMA5 1_67 0.006 859.91 1.4e-05 -29.65
6970 ATXN7L2 1_67 0.006 268.37 4.4e-06 -16.39
4435 PSRC1 1_67 1.000 1162.99 3.4e-03 -34.57
8615 CYB561D1 1_67 0.088 111.35 2.9e-05 9.70
9259 AMIGO1 1_67 0.008 20.79 5.0e-07 -3.57
6445 GPR61 1_67 0.015 36.64 1.6e-06 5.29
587 GNAI3 1_67 0.012 17.10 5.8e-07 -2.63
7977 GSTM4 1_67 0.006 12.36 2.1e-07 2.68
10821 GSTM2 1_67 0.006 9.24 1.5e-07 2.09
4430 GSTM1 1_67 0.006 14.64 2.7e-07 2.97
4432 GSTM5 1_67 0.007 6.42 1.2e-07 0.35
4433 GSTM3 1_67 0.006 19.63 3.3e-07 -3.84
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 1_39"
genename region_tag susie_pip mu2 PVE z
6956 TM2D1 1_39 0.076 27.56 6.1e-06 2.48
4316 KANK4 1_39 0.008 5.82 1.4e-07 0.89
6957 USP1 1_39 0.771 484.88 1.1e-03 22.52
4317 ANGPTL3 1_39 0.238 482.73 3.3e-04 22.47
3024 DOCK7 1_39 0.023 60.05 4.1e-06 7.04
3733 ATG4C 1_39 0.035 143.19 1.4e-05 -11.87
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 2_27"
genename region_tag susie_pip mu2 PVE z
12661 LINC01126 2_27 0 13.09 5.8e-10 0.52
2977 THADA 2_27 0 7.35 8.5e-11 -2.11
6208 PLEKHH2 2_27 0 16.88 1.1e-09 -2.85
11022 C1GALT1C1L 2_27 0 33.76 3.7e-09 3.46
4930 DYNC2LI1 2_27 0 7.99 6.0e-11 0.03
5563 ABCG8 2_27 1 280.94 8.2e-04 -18.98
4943 LRPPRC 2_27 0 11.95 2.4e-10 -0.96
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 11_34"
genename region_tag susie_pip mu2 PVE z
9982 FAM111B 11_34 0.003 5.89 5.6e-08 -0.39
7662 FAM111A 11_34 0.007 13.74 3.0e-07 1.26
2444 DTX4 11_34 0.003 5.61 5.1e-08 0.46
10267 MPEG1 11_34 0.003 4.93 4.2e-08 0.03
7684 PATL1 11_34 0.041 30.88 3.6e-06 3.06
7687 STX3 11_34 0.003 5.39 4.8e-08 -0.20
7688 MRPL16 11_34 0.005 9.96 1.5e-07 1.04
5997 MS4A2 11_34 0.005 9.18 1.2e-07 -0.99
2453 MS4A6A 11_34 0.003 5.55 5.1e-08 0.44
10924 MS4A4E 11_34 0.004 8.24 9.6e-08 1.00
7698 MS4A14 11_34 0.108 39.28 1.2e-05 -3.02
7697 MS4A7 11_34 0.003 5.14 4.4e-08 -0.36
2455 CCDC86 11_34 0.003 6.47 6.5e-08 -0.43
2456 PRPF19 11_34 0.004 8.29 9.6e-08 1.04
2457 TMEM109 11_34 0.003 6.62 6.6e-08 0.66
2480 SLC15A3 11_34 0.007 14.76 2.8e-07 1.91
2481 CD5 11_34 0.004 6.75 7.1e-08 0.16
7874 VPS37C 11_34 0.003 5.81 5.4e-08 0.46
7875 VWCE 11_34 0.003 5.28 4.5e-08 -0.51
5990 TMEM138 11_34 0.007 18.11 3.5e-07 -2.63
6902 CYB561A3 11_34 0.007 18.11 3.5e-07 -2.63
9789 TMEM216 11_34 0.003 5.35 4.6e-08 -0.53
11817 RP11-286N22.8 11_34 0.005 10.65 1.6e-07 -1.05
5996 CPSF7 11_34 0.003 9.33 7.9e-08 -2.11
6903 PPP1R32 11_34 0.009 13.50 3.4e-07 0.14
11812 RP11-794G24.1 11_34 0.022 21.31 1.4e-06 0.74
4508 TMEM258 11_34 0.009 78.66 2.0e-06 -8.23
4507 FADS2 11_34 0.004 260.01 2.7e-06 16.21
7955 FEN1 11_34 0.004 260.01 2.7e-06 16.21
5991 FADS1 11_34 1.000 306.16 8.9e-04 17.68
1196 GANAB 11_34 0.004 114.26 1.2e-06 -10.56
11004 FADS3 11_34 0.011 34.34 1.1e-06 4.46
7876 BEST1 11_34 0.004 36.62 4.2e-07 -5.50
3676 DKFZP434K028 11_34 0.004 12.28 1.4e-07 2.34
5994 INCENP 11_34 0.003 8.10 7.1e-08 -1.72
6904 ASRGL1 11_34 0.003 5.39 4.6e-08 -0.58
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 5_45"
genename region_tag susie_pip mu2 PVE z
8340 ENC1 5_45 0.001 4.94 1.0e-08 -0.18
7307 GFM2 5_45 0.001 5.90 1.4e-08 -0.24
7306 NSA2 5_45 0.001 11.39 4.0e-08 -1.64
10458 FAM169A 5_45 0.001 5.40 1.1e-08 -0.75
3441 POLK 5_45 0.003 215.49 2.1e-06 17.18
12287 CTC-366B18.4 5_45 0.039 102.38 1.2e-05 -10.44
9978 ANKDD1B 5_45 0.003 134.31 1.3e-06 14.36
6186 POC5 5_45 0.009 63.75 1.7e-06 -7.46
11757 AC113404.1 5_45 0.001 15.41 6.7e-08 2.46
5717 IQGAP2 5_45 0.001 9.17 2.7e-08 1.49
7281 F2RL2 5_45 0.001 5.73 1.3e-08 0.40
9219 F2R 5_45 0.001 8.48 2.6e-08 -0.88
7287 F2RL1 5_45 0.004 20.79 2.2e-07 1.84
5718 CRHBP 5_45 0.001 5.46 1.2e-08 -0.35
7288 AGGF1 5_45 0.001 8.91 2.8e-08 -0.92
4314 ZBED3 5_45 0.001 10.59 3.9e-08 -1.08
2729 PDE8B 5_45 0.001 6.43 1.6e-08 0.57
7289 WDR41 5_45 0.001 5.76 1.3e-08 -0.43
4313 AP3B1 5_45 0.018 35.83 1.8e-06 2.59
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
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
14626 rs11580527 1_34 1.000 73.69 2.1e-04 -10.13
14656 rs2495502 1_34 1.000 240.18 7.0e-04 5.54
14667 rs10888896 1_34 1.000 111.21 3.2e-04 10.83
14674 rs471705 1_34 1.000 184.84 5.4e-04 15.19
68511 rs11679386 2_12 1.000 106.04 3.1e-04 10.53
68560 rs1042034 2_13 1.000 114.09 3.3e-04 10.94
68566 rs934197 2_13 1.000 347.10 1.0e-03 29.22
68569 rs548145 2_13 1.000 557.41 1.6e-03 30.76
68646 rs1848922 2_13 1.000 191.64 5.6e-04 23.18
75620 rs72800939 2_28 1.000 46.32 1.3e-04 -7.11
158732 rs768688512 3_58 1.000 355.59 1.0e-03 2.57
225978 rs35518360 4_67 1.000 64.02 1.9e-04 -8.24
226044 rs13140033 4_68 1.000 50.70 1.5e-04 -7.51
315923 rs11376017 6_13 1.000 56.84 1.7e-04 -7.80
319652 rs115740542 6_20 1.000 143.64 4.2e-04 -11.58
320822 rs3130253 6_23 1.000 37.49 1.1e-04 6.29
344230 rs9496567 6_67 1.000 41.42 1.2e-04 -6.54
362797 rs12208357 6_103 1.000 207.53 6.0e-04 11.52
362963 rs191555775 6_104 1.000 66.64 1.9e-04 -7.78
427084 rs1495743 8_20 1.000 63.86 1.9e-04 -8.40
437357 rs140753685 8_42 1.000 58.54 1.7e-04 8.10
438753 rs4738679 8_45 1.000 116.09 3.4e-04 -12.10
458416 rs13252684 8_83 1.000 247.43 7.2e-04 12.43
491297 rs2777798 9_52 1.000 76.85 2.2e-04 8.30
498670 rs115478735 9_70 1.000 302.95 8.8e-04 18.70
523337 rs569165969 10_46 1.000 581.90 1.7e-03 3.04
523385 rs6480402 10_46 1.000 71.70 2.1e-04 -4.33
587218 rs1805735 12_9 1.000 44.45 1.3e-04 6.51
622870 rs11057830 12_76 1.000 38.59 1.1e-04 4.54
723320 rs66495554 16_31 1.000 58.57 1.7e-04 -2.33
751913 rs1801689 17_38 1.000 45.53 1.3e-04 6.91
752829 rs113408695 17_39 1.000 106.68 3.1e-04 10.72
752855 rs8070232 17_39 1.000 125.40 3.6e-04 -6.91
786250 rs73013176 19_9 1.000 198.97 5.8e-04 -14.67
786288 rs137992968 19_9 1.000 87.34 2.5e-04 -9.15
786331 rs12610628 19_10 1.000 63.91 1.9e-04 4.08
786338 rs12459555 19_10 1.000 155.18 4.5e-04 -8.87
786346 rs375323 19_10 1.000 46.16 1.3e-04 3.46
789060 rs2285626 19_15 1.000 298.01 8.7e-04 -20.00
789085 rs3794991 19_15 1.000 261.36 7.6e-04 -23.77
796501 rs62115478 19_30 1.000 123.05 3.6e-04 -11.69
796784 rs116881820 19_32 1.000 1035.31 3.0e-03 21.90
796788 rs34878901 19_32 1.000 6692.66 1.9e-02 13.17
796789 rs405509 19_32 1.000 6706.87 1.9e-02 -24.81
796793 rs814573 19_32 1.000 1943.38 5.6e-03 47.81
796798 rs12721109 19_32 1.000 479.18 1.4e-03 -37.23
807072 rs34507316 20_13 1.000 61.44 1.8e-04 -5.79
905851 rs10677712 2_69 1.000 6957.11 2.0e-02 5.51
1041306 rs964184 11_70 1.000 230.55 6.7e-04 -19.62
1094811 rs374141296 19_34 1.000 10990.56 3.2e-02 3.89
319631 rs72834643 6_20 0.999 36.76 1.1e-04 -5.12
362945 rs117733303 6_104 0.999 75.30 2.2e-04 8.54
618780 rs653178 12_67 0.999 137.16 4.0e-04 12.87
789044 rs12981966 19_15 0.999 108.45 3.1e-04 1.93
811259 rs73124945 20_24 0.999 38.14 1.1e-04 -8.30
771635 rs118043171 18_27 0.998 95.27 2.8e-04 13.33
796456 rs73036721 19_30 0.998 29.60 8.6e-05 -5.27
1094808 rs113176985 19_34 0.998 10934.53 3.2e-02 4.13
622863 rs3782287 12_76 0.997 32.27 9.3e-05 -5.97
755988 rs4969183 17_44 0.997 79.23 2.3e-04 9.38
193456 rs36205397 4_4 0.996 57.56 1.7e-04 7.28
807071 rs6075251 20_13 0.996 42.32 1.2e-04 -2.26
195681 rs2002574 4_10 0.994 28.40 8.2e-05 -5.24
630249 rs79490353 13_7 0.993 28.25 8.1e-05 -5.26
727207 rs4396539 16_37 0.993 30.82 8.9e-05 -5.09
136305 rs709149 3_9 0.992 50.58 1.5e-04 -8.29
399876 rs3197597 7_61 0.991 34.13 9.8e-05 -4.86
466469 rs7024888 9_3 0.991 28.85 8.3e-05 -5.31
143315 rs9834932 3_24 0.990 59.70 1.7e-04 -8.06
788725 rs2302209 19_14 0.990 35.86 1.0e-04 6.01
241877 rs114756490 4_100 0.989 27.51 7.9e-05 5.19
786279 rs1569372 19_9 0.989 205.33 5.9e-04 8.06
811206 rs6029132 20_24 0.989 41.60 1.2e-04 -6.82
191669 rs5855544 3_120 0.988 26.29 7.5e-05 -5.02
786271 rs147985405 19_9 0.988 1918.62 5.5e-03 -44.73
438721 rs56386732 8_45 0.987 34.07 9.8e-05 -6.98
936743 rs9884390 4_48 0.987 106.07 3.0e-04 11.39
1059648 rs10468017 15_26 0.987 207.50 5.9e-04 21.43
1078033 rs2908806 17_7 0.984 41.69 1.2e-04 -6.83
587214 rs4883198 12_9 0.983 40.46 1.2e-04 -6.24
458405 rs79658059 8_83 0.982 280.71 8.0e-04 -15.46
1041510 rs525028 11_70 0.982 134.15 3.8e-04 -17.22
999963 rs1800978 9_53 0.981 193.49 5.5e-04 -13.33
839227 rs145678077 22_17 0.979 26.38 7.5e-05 -5.39
811255 rs76981217 20_24 0.978 34.48 9.8e-05 7.65
562380 rs6591179 11_36 0.977 31.35 8.9e-05 5.83
900226 rs1260326 2_16 0.977 239.32 6.8e-04 -19.02
320118 rs6908155 6_21 0.976 25.70 7.3e-05 3.99
606560 rs148481241 12_44 0.974 24.94 7.1e-05 4.78
746170 rs55764662 17_26 0.974 25.10 7.1e-05 -4.77
277325 rs7701166 5_45 0.970 28.76 8.1e-05 -2.06
422761 rs117037226 8_11 0.970 33.98 9.6e-05 5.56
871174 rs2642438 1_112 0.968 112.00 3.1e-04 11.18
771854 rs74461650 18_28 0.966 29.40 8.2e-05 5.38
219147 rs1458038 4_54 0.965 48.58 1.4e-04 -7.15
473666 rs1556516 9_16 0.964 67.61 1.9e-04 -8.64
616873 rs1196760 12_63 0.963 26.53 7.4e-05 -4.99
490780 rs150108287 9_52 0.960 24.94 7.0e-05 4.62
68363 rs6531234 2_12 0.958 43.38 1.2e-04 -7.99
383547 rs141379002 7_33 0.957 25.19 7.0e-05 4.89
771650 rs62101781 18_27 0.954 77.11 2.1e-04 9.64
637180 rs201796 13_21 0.953 27.85 7.7e-05 -5.29
94109 rs1002015 2_66 0.944 27.11 7.4e-05 -4.75
536738 rs12244851 10_70 0.942 35.38 9.7e-05 -4.50
7471 rs79598313 1_18 0.939 23.79 6.5e-05 4.53
362988 rs374071816 6_104 0.936 123.31 3.4e-04 14.36
999627 rs4149307 9_53 0.927 172.99 4.7e-04 12.99
827139 rs2835302 21_17 0.925 26.90 7.2e-05 -4.92
905847 rs28850371 2_69 0.923 6951.55 1.9e-02 5.16
798943 rs397558 19_37 0.922 47.58 1.3e-04 7.13
786320 rs66536924 19_10 0.918 60.85 1.6e-04 4.58
167488 rs189174 3_74 0.911 59.00 1.6e-04 8.04
742371 rs117859452 17_17 0.908 23.97 6.3e-05 -4.15
581309 rs74612335 11_77 0.904 86.70 2.3e-04 10.00
319470 rs75080831 6_19 0.902 51.81 1.4e-04 -7.48
277266 rs10062361 5_45 0.897 188.75 4.9e-04 19.22
575023 rs201912654 11_59 0.897 53.30 1.4e-04 -7.45
786274 rs3745677 19_9 0.883 73.00 1.9e-04 8.19
622705 rs35480942 12_75 0.880 26.05 6.7e-05 -4.85
819896 rs62219001 21_2 0.879 23.20 5.9e-05 -4.40
12757 rs138863615 1_29 0.877 24.33 6.2e-05 4.53
811225 rs6129631 20_24 0.876 89.34 2.3e-04 -10.50
771631 rs7241918 18_27 0.871 130.39 3.3e-04 14.61
320820 rs41291774 6_23 0.870 39.02 9.9e-05 6.84
399867 rs11761624 7_61 0.868 25.97 6.5e-05 -3.22
810000 rs11167269 20_21 0.858 70.11 1.7e-04 -8.89
302172 rs12657266 5_92 0.855 201.58 5.0e-04 15.97
828276 rs149577713 21_19 0.852 33.29 8.2e-05 3.86
295077 rs546280079 5_79 0.850 27.67 6.8e-05 -4.96
634443 rs148480921 13_16 0.850 28.75 7.1e-05 -5.24
23959 rs11161548 1_52 0.845 25.61 6.3e-05 -4.87
39428 rs1795240 1_84 0.844 28.88 7.1e-05 -5.20
153306 rs6762369 3_47 0.844 31.91 7.8e-05 5.52
68563 rs78610189 2_13 0.843 53.26 1.3e-04 -8.16
225260 rs148447389 4_67 0.838 25.06 6.1e-05 4.50
231043 rs138204164 4_77 0.837 31.21 7.6e-05 -5.41
741289 rs116139938 17_16 0.831 23.26 5.6e-05 -4.22
871175 rs17850677 1_112 0.831 29.09 7.0e-05 -4.38
14657 rs1887552 1_34 0.819 276.46 6.6e-04 -8.86
730110 rs78864981 16_44 0.818 24.60 5.8e-05 4.34
171382 rs73215911 3_82 0.811 25.07 5.9e-05 -4.76
277289 rs3843482 5_45 0.810 374.00 8.8e-04 24.18SNPs 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")
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
#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
1094811 rs374141296 19_34 1.000 10990.56 3.2e-02 3.89
1094799 rs61371437 19_34 0.000 10981.84 0.0e+00 3.94
1094789 rs739349 19_34 0.000 10972.17 0.0e+00 3.93
1094790 rs756628 19_34 0.000 10971.64 0.0e+00 3.92
1094786 rs739347 19_34 0.000 10955.21 0.0e+00 3.86
1094787 rs2073614 19_34 0.000 10941.27 0.0e+00 3.84
1094782 rs4802613 19_34 0.000 10939.05 0.0e+00 3.88
1094808 rs113176985 19_34 0.998 10934.53 3.2e-02 4.13
1094801 rs35295508 19_34 0.000 10932.32 1.9e-09 4.12
1094780 rs10403394 19_34 0.000 10929.53 0.0e+00 3.91
1094781 rs17555056 19_34 0.000 10915.39 0.0e+00 3.86
1094792 rs2077300 19_34 0.000 10906.92 0.0e+00 4.04
1094815 rs2946865 19_34 0.000 10884.87 8.7e-07 4.14
1094796 rs73056059 19_34 0.000 10883.53 0.0e+00 4.01
1094806 rs73056069 19_34 0.002 10881.46 6.1e-05 4.26
1094803 rs2878354 19_34 0.000 10867.57 5.6e-09 4.22
1094816 rs60815603 19_34 0.000 10722.19 9.2e-14 4.05
1094819 rs1316885 19_34 0.000 10675.54 1.3e-14 4.10
1094821 rs60746284 19_34 0.000 10656.10 1.1e-10 4.24
1094824 rs2946863 19_34 0.000 10652.65 2.3e-13 4.14
1094817 rs35443645 19_34 0.000 10642.05 1.4e-17 3.99
1094797 rs73056062 19_34 0.000 10584.11 0.0e+00 3.52
1094827 rs553431297 19_34 0.000 10368.20 0.0e+00 3.71
1094810 rs112283514 19_34 0.000 10358.94 0.0e+00 4.15
1094812 rs11270139 19_34 0.000 10290.64 0.0e+00 3.88
1094777 rs10421294 19_34 0.000 9884.06 0.0e+00 3.91
1094776 rs8108175 19_34 0.000 9882.86 0.0e+00 3.91
1094769 rs59192944 19_34 0.000 9866.53 0.0e+00 3.95
1094775 rs1858742 19_34 0.000 9864.13 0.0e+00 3.93
1094766 rs55991145 19_34 0.000 9859.00 0.0e+00 3.93
1094761 rs3786567 19_34 0.000 9855.74 0.0e+00 3.93
1094757 rs2271952 19_34 0.000 9852.37 0.0e+00 3.93
1094760 rs4801801 19_34 0.000 9852.10 0.0e+00 3.96
1094756 rs2271953 19_34 0.000 9840.62 0.0e+00 3.92
1094758 rs2271951 19_34 0.000 9839.85 0.0e+00 3.90
1094747 rs60365978 19_34 0.000 9833.68 0.0e+00 3.93
1094753 rs4802612 19_34 0.000 9794.96 0.0e+00 4.15
1094763 rs2517977 19_34 0.000 9750.01 0.0e+00 3.77
1094750 rs55893003 19_34 0.000 9749.39 0.0e+00 4.03
1094742 rs55992104 19_34 0.000 9531.65 0.0e+00 3.29
1094736 rs60403475 19_34 0.000 9530.91 0.0e+00 3.34
1094739 rs4352151 19_34 0.000 9529.04 0.0e+00 3.33
1094733 rs11878448 19_34 0.000 9523.04 0.0e+00 3.34
1094727 rs9653100 19_34 0.000 9520.51 0.0e+00 3.31
1094723 rs4802611 19_34 0.000 9514.65 0.0e+00 3.32
1094715 rs7251338 19_34 0.000 9502.32 0.0e+00 3.33
1094714 rs59269605 19_34 0.000 9501.30 0.0e+00 3.33
1094735 rs1042120 19_34 0.000 9473.40 0.0e+00 3.51
1094731 rs113220577 19_34 0.000 9465.47 0.0e+00 3.50
1094725 rs9653118 19_34 0.000 9450.89 0.0e+00 3.48SNPs 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
1094808 rs113176985 19_34 0.998 10934.53 0.03200 4.13
1094811 rs374141296 19_34 1.000 10990.56 0.03200 3.89
905851 rs10677712 2_69 1.000 6957.11 0.02000 5.51
796788 rs34878901 19_32 1.000 6692.66 0.01900 13.17
796789 rs405509 19_32 1.000 6706.87 0.01900 -24.81
905847 rs28850371 2_69 0.923 6951.55 0.01900 5.16
905808 rs6542385 2_69 0.345 6944.41 0.00700 5.17
796793 rs814573 19_32 1.000 1943.38 0.00560 47.81
786271 rs147985405 19_9 0.988 1918.62 0.00550 -44.73
796784 rs116881820 19_32 1.000 1035.31 0.00300 21.90
905855 rs1546621 2_69 0.100 6937.30 0.00200 5.16
523337 rs569165969 10_46 1.000 581.90 0.00170 3.04
68569 rs548145 2_13 1.000 557.41 0.00160 30.76
796798 rs12721109 19_32 1.000 479.18 0.00140 -37.23
905860 rs1516006 2_69 0.069 6936.77 0.00140 5.15
523338 rs7909631 10_46 0.680 610.79 0.00120 2.48
68566 rs934197 2_13 1.000 347.10 0.00100 29.22
158732 rs768688512 3_58 1.000 355.59 0.00100 2.57
277289 rs3843482 5_45 0.810 374.00 0.00088 24.18
498670 rs115478735 9_70 1.000 302.95 0.00088 18.70
789060 rs2285626 19_15 1.000 298.01 0.00087 -20.00
458405 rs79658059 8_83 0.982 280.71 0.00080 -15.46
789085 rs3794991 19_15 1.000 261.36 0.00076 -23.77
458416 rs13252684 8_83 1.000 247.43 0.00072 12.43
14656 rs2495502 1_34 1.000 240.18 0.00070 5.54
900226 rs1260326 2_16 0.977 239.32 0.00068 -19.02
1041306 rs964184 11_70 1.000 230.55 0.00067 -19.62
14657 rs1887552 1_34 0.819 276.46 0.00066 -8.86
362797 rs12208357 6_103 1.000 207.53 0.00060 11.52
786279 rs1569372 19_9 0.989 205.33 0.00059 8.06
1059648 rs10468017 15_26 0.987 207.50 0.00059 21.43
723315 rs821840 16_31 0.673 294.55 0.00058 16.71
786250 rs73013176 19_9 1.000 198.97 0.00058 -14.67
68646 rs1848922 2_13 1.000 191.64 0.00056 23.18
999963 rs1800978 9_53 0.981 193.49 0.00055 -13.33
14674 rs471705 1_34 1.000 184.84 0.00054 15.19
523336 rs7084697 10_46 0.288 610.05 0.00051 2.43
302172 rs12657266 5_92 0.855 201.58 0.00050 15.97
277266 rs10062361 5_45 0.897 188.75 0.00049 19.22
999627 rs4149307 9_53 0.927 172.99 0.00047 12.99
786338 rs12459555 19_10 1.000 155.18 0.00045 -8.87
319652 rs115740542 6_20 1.000 143.64 0.00042 -11.58
618780 rs653178 12_67 0.999 137.16 0.00040 12.87
904638 rs6544713 2_27 0.643 214.79 0.00040 -19.78
1041510 rs525028 11_70 0.982 134.15 0.00038 -17.22
158728 rs73141241 3_58 0.376 343.70 0.00037 2.73
752855 rs8070232 17_39 1.000 125.40 0.00036 -6.91
796501 rs62115478 19_30 1.000 123.05 0.00036 -11.69
362988 rs374071816 6_104 0.936 123.31 0.00034 14.36
438753 rs4738679 8_45 1.000 116.09 0.00034 -12.10SNPs with largest z scores
#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
796793 rs814573 19_32 1.000 1943.38 5.6e-03 47.81
786271 rs147985405 19_9 0.988 1918.62 5.5e-03 -44.73
786266 rs73015020 19_9 0.009 1911.20 4.9e-05 -44.63
786264 rs138175288 19_9 0.001 1905.92 5.1e-06 -44.58
786267 rs77140532 19_9 0.001 1907.18 6.1e-06 -44.58
786265 rs138294113 19_9 0.000 1903.57 2.3e-06 -44.56
786268 rs112552009 19_9 0.000 1902.12 1.6e-06 -44.53
786269 rs10412048 19_9 0.000 1902.46 6.9e-07 -44.53
786263 rs55997232 19_9 0.000 1883.86 2.5e-09 -44.35
786272 rs17248769 19_9 0.000 1435.60 6.2e-08 -37.32
796798 rs12721109 19_32 1.000 479.18 1.4e-03 -37.23
786273 rs2228671 19_9 0.000 1425.26 4.4e-08 -37.18
856240 rs12740374 1_67 0.000 1019.85 1.1e-06 -34.67
856247 rs646776 1_67 0.000 1015.00 9.7e-07 34.60
856236 rs7528419 1_67 0.000 1015.29 1.0e-06 -34.59
856246 rs629301 1_67 0.000 1013.01 9.3e-07 34.57
796742 rs62117204 19_32 0.000 998.79 0.0e+00 -34.15
856261 rs4970836 1_67 0.000 982.29 9.4e-07 34.03
856258 rs583104 1_67 0.000 981.62 9.2e-07 34.02
856263 rs1277930 1_67 0.000 980.03 9.5e-07 33.99
856264 rs599839 1_67 0.000 979.48 9.7e-07 33.97
856244 rs3832016 1_67 0.000 951.20 7.7e-07 33.47
856241 rs660240 1_67 0.000 948.38 7.8e-07 33.42
856259 rs602633 1_67 0.000 929.19 7.8e-07 33.07
796729 rs1551891 19_32 0.000 1042.49 0.0e+00 -32.83
786262 rs9305020 19_9 0.000 1070.53 7.4e-13 -31.79
68569 rs548145 2_13 1.000 557.41 1.6e-03 30.76
68566 rs934197 2_13 1.000 347.10 1.0e-03 29.22
856227 rs4970834 1_67 0.001 710.86 1.0e-06 -28.83
68570 rs478588 2_13 0.000 511.74 2.9e-12 28.40
68597 rs532300 2_13 0.000 469.63 1.1e-12 27.62
68598 rs558130 2_13 0.000 469.61 1.1e-12 27.62
68599 rs533211 2_13 0.000 469.61 1.1e-12 27.62
68600 rs528113 2_13 0.000 469.27 1.0e-12 27.62
68620 rs574461 2_13 0.000 469.43 1.1e-12 27.62
68622 rs494465 2_13 0.000 469.33 1.1e-12 27.62
68595 rs312979 2_13 0.000 468.69 9.9e-13 27.61
68601 rs547179 2_13 0.000 469.43 1.1e-12 27.61
68605 rs1652418 2_13 0.000 468.30 9.9e-13 27.60
68607 rs563696 2_13 0.000 468.27 9.9e-13 27.60
68609 rs479545 2_13 0.000 467.98 9.7e-13 27.60
68606 rs560844 2_13 0.000 467.60 9.5e-13 27.59
68608 rs475887 2_13 0.000 467.60 9.5e-13 27.59
68603 rs1712250 2_13 0.000 466.84 9.1e-13 27.58
68610 rs548594 2_13 0.000 465.92 8.7e-13 27.56
68615 rs561850 2_13 0.000 465.80 8.7e-13 27.56
68616 rs477146 2_13 0.000 465.80 8.7e-13 27.56
68618 rs483621 2_13 0.000 465.88 8.8e-13 27.56
68611 rs488507 2_13 0.000 465.72 8.7e-13 27.55
68617 rs13411597 2_13 0.000 465.72 8.7e-13 27.55Gene 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] 36if (length(genes)>0){
GO_enrichment <- enrichr(genes, dbs)
for (db in dbs){
print(db)
df <- 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"
Term
1 peptidyl-serine phosphorylation (GO:0018105)
2 peptidyl-serine modification (GO:0018209)
3 protein phosphorylation (GO:0006468)
4 negative regulation of cholesterol storage (GO:0010887)
5 cholesterol homeostasis (GO:0042632)
6 sterol homeostasis (GO:0055092)
7 regulation of cholesterol storage (GO:0010885)
8 cellular protein modification process (GO:0006464)
9 positive regulation of cyclin-dependent protein serine/threonine kinase activity (GO:0045737)
10 activin receptor signaling pathway (GO:0032924)
11 positive regulation of cyclin-dependent protein kinase activity (GO:1904031)
12 negative regulation of lipid storage (GO:0010888)
13 phosphorylation (GO:0016310)
14 cholesterol efflux (GO:0033344)
15 regulation of DNA biosynthetic process (GO:2000278)
Overlap Adjusted.P.value Genes
1 5/156 0.002841451 CSNK1G3;TNKS;PKN3;PRKD2;GAS6
2 5/169 0.002841451 CSNK1G3;TNKS;PKN3;PRKD2;GAS6
3 7/496 0.003785591 CSNK1G3;PXK;ACVR1C;TNKS;PKN3;PRKD2;GAS6
4 2/10 0.016119155 ABCA1;TTC39B
5 3/71 0.021480586 ABCA1;ABCG8;TTC39B
6 3/72 0.021480586 ABCA1;ABCG8;TTC39B
7 2/16 0.021480586 ABCA1;TTC39B
8 8/1025 0.021480586 CSNK1G3;PXK;ACVR1C;TNKS;PKN3;PRKD2;FUT2;GAS6
9 2/17 0.021480586 CCND2;PSRC1
10 2/19 0.022430967 ACVR1C;INHBB
11 2/20 0.022430967 CCND2;PSRC1
12 2/20 0.022430967 ABCA1;TTC39B
13 5/400 0.024929070 PXK;ACVR1C;PKN3;PRKD2;GAS6
14 2/24 0.027803001 ABCA1;ABCG8
15 2/29 0.037956838 TNKS;PRKD2
[1] "GO_Cellular_Component_2021"
[1] Term Overlap Adjusted.P.value Genes
<0 rows> (or 0-length row.names)
[1] "GO_Molecular_Function_2021"
Term Overlap Adjusted.P.value
1 cholesterol transfer activity (GO:0120020) 2/18 0.01822928
2 sterol transfer activity (GO:0120015) 2/19 0.01822928
Genes
1 ABCA1;ABCG8
2 ABCA1;ABCG8PELO gene(s) from the input list not found in DisGeNET CURATEDSPTY2D1 gene(s) from the input list not found in DisGeNET CURATEDRPS11 gene(s) from the input list not found in DisGeNET CURATEDC10orf88 gene(s) from the input list not found in DisGeNET CURATEDTRIM39 gene(s) from the input list not found in DisGeNET CURATEDALLC gene(s) from the input list not found in DisGeNET CURATEDACP6 gene(s) from the input list not found in DisGeNET CURATEDDDX56 gene(s) from the input list not found in DisGeNET CURATEDCRACR2B gene(s) from the input list not found in DisGeNET CURATEDTTC39B gene(s) from the input list not found in DisGeNET CURATEDPSRC1 gene(s) from the input list not found in DisGeNET CURATEDCSNK1G3 gene(s) from the input list not found in DisGeNET CURATEDTNKS gene(s) from the input list not found in DisGeNET CURATEDNYNRIN gene(s) from the input list not found in DisGeNET CURATEDHPR gene(s) from the input list not found in DisGeNET CURATEDAKNA gene(s) from the input list not found in DisGeNET CURATEDPKN3 gene(s) from the input list not found in DisGeNET CURATEDRP4-781K5.7 gene(s) from the input list not found in DisGeNET CURATEDCNIH4 gene(s) from the input list not found in DisGeNET CURATED Description FDR Ratio
4 Blood Platelet Disorders 0.01360544 2/17
26 Hypercholesterolemia, Familial 0.01360544 2/17
42 Tangier Disease 0.01360544 1/17
56 Caliciviridae Infections 0.01360544 1/17
57 Infections, Calicivirus 0.01360544 1/17
68 Charcot-Marie-Tooth disease, Type 1C 0.01360544 1/17
76 Hypoalphalipoproteinemias 0.01360544 1/17
83 Tangier Disease Neuropathy 0.01360544 1/17
111 GALLBLADDER DISEASE 4 0.01360544 1/17
114 VITAMIN B12 PLASMA LEVEL QUANTITATIVE TRAIT LOCUS 1 0.01360544 1/17
BgRatio
4 16/9703
26 18/9703
42 1/9703
56 1/9703
57 1/9703
68 1/9703
76 1/9703
83 1/9703
111 1/9703
114 1/9703******************************************
* *
* Welcome to WebGestaltR ! *
* *
******************************************Loading the functional categories...
Loading the ID list...
Loading the reference list...
Performing the enrichment analysis...
description size overlap FDR database
1 Coronary Artery Disease 153 9 9.942959e-06 disease_GLAD4U
2 Dyslipidaemia 84 6 8.174782e-04 disease_GLAD4U
3 Coronary Disease 171 7 2.070116e-03 disease_GLAD4U
4 Arteriosclerosis 173 7 2.070116e-03 disease_GLAD4U
5 Arterial Occlusive Diseases 174 6 2.216212e-02 disease_GLAD4U
6 Myocardial Ischemia 181 6 2.302842e-02 disease_GLAD4U
userId
1 PSRC1;ABCG8;INSIG2;TTC39B;ABCA1;SPTY2D1;FADS1;NYNRIN;FUT2
2 PSRC1;ABCG8;INSIG2;TTC39B;ABCA1;FADS1
3 PSRC1;ABCG8;INSIG2;TTC39B;ABCA1;FADS1;FUT2
4 PSRC1;ABCG8;TTC39B;ABCA1;FADS1;NYNRIN;HPR
5 PSRC1;ABCG8;TTC39B;ABCA1;FADS1;NYNRIN
6 PSRC1;ABCG8;INSIG2;TTC39B;ABCA1;FADS1
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] stats graphics grDevices utils datasets methods base
other attached packages:
[1] WebGestaltR_0.4.4 disgenet2r_0.99.2 enrichR_3.0 cowplot_1.0.0
[5] ggplot2_3.3.3
loaded via a namespace (and not attached):
[1] bitops_1.0-6 matrixStats_0.57.0
[3] fs_1.3.1 bit64_4.0.5
[5] doParallel_1.0.16 progress_1.2.2
[7] httr_1.4.1 rprojroot_2.0.2
[9] GenomeInfoDb_1.20.0 doRNG_1.8.2
[11] tools_3.6.1 utf8_1.2.1
[13] R6_2.5.0 DBI_1.1.1
[15] BiocGenerics_0.30.0 colorspace_1.4-1
[17] withr_2.4.1 tidyselect_1.1.0
[19] prettyunits_1.0.2 bit_4.0.4
[21] curl_3.3 compiler_3.6.1
[23] git2r_0.26.1 Biobase_2.44.0
[25] DelayedArray_0.10.0 rtracklayer_1.44.0
[27] labeling_0.3 scales_1.1.0
[29] readr_1.4.0 apcluster_1.4.8
[31] stringr_1.4.0 digest_0.6.20
[33] Rsamtools_2.0.0 svglite_1.2.2
[35] rmarkdown_1.13 XVector_0.24.0
[37] pkgconfig_2.0.3 htmltools_0.3.6
[39] fastmap_1.1.0 BSgenome_1.52.0
[41] rlang_0.4.11 RSQLite_2.2.7
[43] generics_0.0.2 farver_2.1.0
[45] jsonlite_1.6 BiocParallel_1.18.0
[47] dplyr_1.0.7 VariantAnnotation_1.30.1
[49] RCurl_1.98-1.1 magrittr_2.0.1
[51] GenomeInfoDbData_1.2.1 Matrix_1.2-18
[53] Rcpp_1.0.6 munsell_0.5.0
[55] S4Vectors_0.22.1 fansi_0.5.0
[57] gdtools_0.1.9 lifecycle_1.0.0
[59] stringi_1.4.3 whisker_0.3-2
[61] yaml_2.2.0 SummarizedExperiment_1.14.1
[63] zlibbioc_1.30.0 plyr_1.8.4
[65] grid_3.6.1 blob_1.2.1
[67] parallel_3.6.1 promises_1.0.1
[69] crayon_1.4.1 lattice_0.20-38
[71] Biostrings_2.52.0 GenomicFeatures_1.36.3
[73] hms_1.1.0 knitr_1.23
[75] pillar_1.6.1 igraph_1.2.4.1
[77] GenomicRanges_1.36.0 rjson_0.2.20
[79] rngtools_1.5 codetools_0.2-16
[81] reshape2_1.4.3 biomaRt_2.40.1
[83] stats4_3.6.1 XML_3.98-1.20
[85] glue_1.4.2 evaluate_0.14
[87] data.table_1.14.0 foreach_1.5.1
[89] vctrs_0.3.8 httpuv_1.5.1
[91] gtable_0.3.0 purrr_0.3.4
[93] assertthat_0.2.1 cachem_1.0.5
[95] xfun_0.8 later_0.8.0
[97] tibble_3.1.2 iterators_1.0.13
[99] GenomicAlignments_1.20.1 AnnotationDbi_1.46.0
[101] memoise_2.0.0 IRanges_2.18.1
[103] workflowr_1.6.2 ellipsis_0.3.2