Apoliprotein A (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
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-30630_irnt_Liver.Rmd) and HTML (docs/ukb-d-30630_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 Apoliprotein A (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-30630_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.0186832223 0.0001879761 #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
23.40897 19.71911 #report sample size
print(sample_size)[1] 313387#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.01521316 0.10287148 #compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))[1] 0.04124006 2.48957722Genes 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
1144 ASAP3 1_16 1.000 52.40 1.7e-04 8.08
12613 GPIHBP1 8_94 1.000 164.31 5.2e-04 -13.58
6391 TTC39B 9_13 1.000 270.90 8.6e-04 -16.60
8531 TNKS 8_12 0.999 117.19 3.7e-04 14.94
7410 ABCA1 9_53 0.994 289.70 9.2e-04 25.98
11684 RP11-136O12.2 8_83 0.987 49.01 1.5e-04 -7.14
5210 PCSK6 15_50 0.985 33.09 1.0e-04 5.12
12467 RP11-219B17.3 15_27 0.982 33.74 1.1e-04 -5.55
7040 INHBB 2_70 0.978 29.32 9.1e-05 4.71
3210 LDAH 2_12 0.976 35.03 1.1e-04 -5.47
10104 SULF2 20_29 0.975 101.39 3.2e-04 -9.53
6509 NTAN1 16_15 0.967 37.16 1.1e-04 -5.83
11699 RP11-10A14.4 8_11 0.957 26.69 8.2e-05 1.99
906 UBE2K 4_32 0.953 122.48 3.7e-04 6.31
9006 BEND3 6_71 0.953 22.05 6.7e-05 -4.36
2731 PCDHB15 5_83 0.951 20.68 6.3e-05 -4.05
2204 AKNA 9_59 0.950 61.27 1.9e-04 -7.00
9478 KMT5A 12_75 0.945 933.88 2.8e-03 -12.14
6100 ALLC 2_2 0.943 48.98 1.5e-04 7.07
9404 PTTG1IP 21_23 0.937 84.80 2.5e-04 9.12
4435 PSRC1 1_67 0.933 227.59 6.8e-04 16.49
3774 ZNF436 1_16 0.925 27.20 8.0e-05 -6.19
10870 IRF9 14_3 0.914 23.26 6.8e-05 4.68
7203 GNL3 3_36 0.913 33.06 9.6e-05 7.22
12229 RP11-346C20.3 16_39 0.912 21.76 6.3e-05 -4.20
4239 TRIM5 11_4 0.906 26.80 7.7e-05 4.08
3855 IDUA 4_2 0.894 121.73 3.5e-04 5.77
12687 RP4-781K5.7 1_121 0.844 81.41 2.2e-04 -8.90
10848 CLIC1 6_26 0.835 115.16 3.1e-04 9.69
1960 SNAPC2 19_7 0.833 23.40 6.2e-05 -4.39
2274 PITRM1 10_4 0.824 19.00 5.0e-05 -3.70
1603 PYGB 20_18 0.817 20.36 5.3e-05 3.94
6792 ADAR 1_75 0.802 53.28 1.4e-04 -7.59Genes 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
7241 MST1R 3_35 0.000 38631.61 0.0e+00 -8.12
35 RBM6 3_35 0.000 37685.47 0.0e+00 9.37
5389 RPS11 19_34 0.000 33558.54 0.0e+00 -4.94
9460 TRAIP 3_35 0.000 30424.48 0.0e+00 -4.69
1227 FLT3LG 19_34 0.000 28943.78 0.0e+00 3.83
11255 GPX1 3_35 0.000 26498.51 0.0e+00 -0.12
656 RHOA 3_35 0.000 26486.69 0.0e+00 -0.11
7235 APEH 3_35 0.000 26033.81 0.0e+00 0.37
9902 SLC38A3 3_35 0.000 16474.17 0.0e+00 5.03
11109 BSN-AS2 3_35 0.000 15062.59 0.0e+00 1.60
7843 CDK2AP2 11_37 0.000 14614.38 1.3e-07 4.29
5666 NICN1 3_35 0.000 14046.25 0.0e+00 -1.59
8555 GMPPB 3_35 0.000 13438.04 0.0e+00 -2.66
11381 LRRC37A2 17_27 0.000 13361.48 0.0e+00 -3.39
8762 RPS6KB2 11_37 0.013 12630.84 5.4e-04 -5.71
5393 RCN3 19_34 0.000 10972.81 0.0e+00 5.07
1931 FCGRT 19_34 0.000 9968.24 0.0e+00 3.34
12683 HCP5B 6_24 0.000 9945.53 2.0e-11 -7.35
7842 NDUFV1 11_37 0.000 9525.28 1.3e-12 1.53
7844 NUDT8 11_37 0.000 9136.69 1.8e-12 1.72Genes 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
9478 KMT5A 12_75 0.945 933.88 0.00280 -12.14
7410 ABCA1 9_53 0.994 289.70 0.00092 25.98
6391 TTC39B 9_13 1.000 270.90 0.00086 -16.60
4435 PSRC1 1_67 0.933 227.59 0.00068 16.49
4317 ANGPTL3 1_39 0.687 254.63 0.00056 16.90
8762 RPS6KB2 11_37 0.013 12630.84 0.00054 -5.71
12613 GPIHBP1 8_94 1.000 164.31 0.00052 -13.58
5991 FADS1 11_34 0.661 201.22 0.00042 14.77
906 UBE2K 4_32 0.953 122.48 0.00037 6.31
8531 TNKS 8_12 0.999 117.19 0.00037 14.94
3855 IDUA 4_2 0.894 121.73 0.00035 5.77
10104 SULF2 20_29 0.975 101.39 0.00032 -9.53
10848 CLIC1 6_26 0.835 115.16 0.00031 9.69
9404 PTTG1IP 21_23 0.937 84.80 0.00025 9.12
12687 RP4-781K5.7 1_121 0.844 81.41 0.00022 -8.90
2204 AKNA 9_59 0.950 61.27 0.00019 -7.00
7656 CATSPER2 15_16 0.545 109.34 0.00019 -10.07
7329 DAGLB 7_9 0.535 101.96 0.00017 10.23
537 DGAT2 11_42 0.501 109.16 0.00017 13.56
1144 ASAP3 1_16 1.000 52.40 0.00017 8.08Genes 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
5240 NLRC5 16_31 0.000 2848.55 5.4e-17 -66.39
1120 CETP 16_31 0.000 1318.73 0.0e+00 -54.09
7547 LIPC 15_26 0.000 1244.96 4.9e-18 -42.31
7410 ABCA1 9_53 0.994 289.70 9.2e-04 25.98
7898 PAFAH1B2 11_70 0.000 334.67 0.0e+00 21.90
3154 APOA1 11_70 0.000 246.24 0.0e+00 -17.98
8739 LPL 8_21 0.000 416.96 0.0e+00 17.48
4317 ANGPTL3 1_39 0.687 254.63 5.6e-04 16.90
6957 USP1 1_39 0.102 249.98 8.1e-05 16.75
6391 TTC39B 9_13 1.000 270.90 8.6e-04 -16.60
4435 PSRC1 1_67 0.933 227.59 6.8e-04 16.49
11738 RP11-115J16.2 8_12 0.016 191.89 9.6e-06 15.73
9360 DDX28 16_36 0.014 215.22 9.6e-06 15.40
6005 SIDT2 11_70 0.000 147.00 0.0e+00 -15.38
8531 TNKS 8_12 0.999 117.19 3.7e-04 14.94
5991 FADS1 11_34 0.661 201.22 4.2e-04 14.77
1894 TRPS1 8_78 0.093 246.49 7.3e-05 -14.52
1231 PABPC4 1_24 0.125 183.15 7.3e-05 14.45
6713 PSKH1 16_36 0.014 169.12 7.3e-06 -14.36
12613 GPIHBP1 8_94 1.000 164.31 5.2e-04 -13.58Comparing 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.02807082#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
5240 NLRC5 16_31 0.000 2848.55 5.4e-17 -66.39
1120 CETP 16_31 0.000 1318.73 0.0e+00 -54.09
7547 LIPC 15_26 0.000 1244.96 4.9e-18 -42.31
7410 ABCA1 9_53 0.994 289.70 9.2e-04 25.98
7898 PAFAH1B2 11_70 0.000 334.67 0.0e+00 21.90
3154 APOA1 11_70 0.000 246.24 0.0e+00 -17.98
8739 LPL 8_21 0.000 416.96 0.0e+00 17.48
4317 ANGPTL3 1_39 0.687 254.63 5.6e-04 16.90
6957 USP1 1_39 0.102 249.98 8.1e-05 16.75
6391 TTC39B 9_13 1.000 270.90 8.6e-04 -16.60
4435 PSRC1 1_67 0.933 227.59 6.8e-04 16.49
11738 RP11-115J16.2 8_12 0.016 191.89 9.6e-06 15.73
9360 DDX28 16_36 0.014 215.22 9.6e-06 15.40
6005 SIDT2 11_70 0.000 147.00 0.0e+00 -15.38
8531 TNKS 8_12 0.999 117.19 3.7e-04 14.94
5991 FADS1 11_34 0.661 201.22 4.2e-04 14.77
1894 TRPS1 8_78 0.093 246.49 7.3e-05 -14.52
1231 PABPC4 1_24 0.125 183.15 7.3e-05 14.45
6713 PSKH1 16_36 0.014 169.12 7.3e-06 -14.36
12613 GPIHBP1 8_94 1.000 164.31 5.2e-04 -13.58Locus 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: 16_31"
genename region_tag susie_pip mu2 PVE z
6688 CES5A 16_31 0 36.06 0.0e+00 4.62
1124 GNAO1 16_31 0 6.27 0.0e+00 -1.31
11561 RP11-461O7.1 16_31 0 9.45 0.0e+00 0.86
6695 AMFR 16_31 0 15.64 0.0e+00 -1.50
7710 NUDT21 16_31 0 9.79 0.0e+00 2.36
3681 BBS2 16_31 0 12.09 0.0e+00 1.27
1122 MT3 16_31 0 12.72 0.0e+00 -4.83
8094 MT1E 16_31 0 17.37 0.0e+00 -0.88
10727 MT1M 16_31 0 81.04 0.0e+00 -7.97
10725 MT1A 16_31 0 11.78 0.0e+00 -4.12
10386 MT1F 16_31 0 38.11 0.0e+00 5.39
9805 MT1X 16_31 0 19.57 0.0e+00 1.82
1740 NUP93 16_31 0 16.15 0.0e+00 -5.10
438 HERPUD1 16_31 0 162.71 0.0e+00 -10.04
1120 CETP 16_31 0 1318.73 0.0e+00 -54.09
5240 NLRC5 16_31 0 2848.55 5.4e-17 -66.39
5239 CPNE2 16_31 0 26.57 0.0e+00 2.21
8472 FAM192A 16_31 0 11.84 0.0e+00 2.03
6698 RSPRY1 16_31 0 17.37 0.0e+00 3.33
1745 PLLP 16_31 0 8.23 0.0e+00 4.95
81 CX3CL1 16_31 0 22.38 0.0e+00 2.63
1747 CCL17 16_31 0 9.91 0.0e+00 2.04
52 CIAPIN1 16_31 0 7.01 0.0e+00 0.46
1154 COQ9 16_31 0 5.07 0.0e+00 -0.60
3685 DOK4 16_31 0 77.59 0.0e+00 -3.36
4628 CCDC102A 16_31 0 6.59 0.0e+00 0.00
10722 ADGRG1 16_31 0 24.57 0.0e+00 2.22
9366 ADGRG3 16_31 0 8.66 0.0e+00 -0.92
5241 KATNB1 16_31 0 35.18 0.0e+00 2.57
5242 KIFC3 16_31 0 5.15 0.0e+00 0.28
1754 USB1 16_31 0 14.42 0.0e+00 -1.45
7571 ZNF319 16_31 0 15.83 0.0e+00 -1.55
1753 MMP15 16_31 0 7.97 0.0e+00 0.83
729 CFAP20 16_31 0 14.99 0.0e+00 -1.49
730 CSNK2A2 16_31 0 5.28 0.0e+00 -0.32
9278 GINS3 16_31 0 5.13 0.0e+00 -0.27
1757 NDRG4 16_31 0 8.00 0.0e+00 -0.83
3680 CNOT1 16_31 0 56.96 0.0e+00 -3.37
1759 SLC38A7 16_31 0 7.74 0.0e+00 0.80
3684 GOT2 16_31 0 50.28 0.0e+00 -3.14
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 15_26"
genename region_tag susie_pip mu2 PVE z
7547 LIPC 15_26 0 1244.96 4.9e-18 -42.31
4905 ADAM10 15_26 0 7.10 2.3e-20 0.20
6536 RNF111 15_26 0 19.27 3.1e-19 0.90
4889 SLTM 15_26 0 62.47 1.1e-17 -5.04
8386 LDHAL6B 15_26 0 7.17 2.5e-20 0.04
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 9_53"
genename region_tag susie_pip mu2 PVE z
7410 ABCA1 9_53 0.994 289.70 9.2e-04 25.98
2193 FKTN 9_53 0.000 56.89 2.0e-20 -3.24
1314 TMEM38B 9_53 0.000 30.00 0.0e+00 -3.46
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 11_70"
genename region_tag susie_pip mu2 PVE z
4868 BUD13 11_70 0 152.74 0.0e+00 -0.84
2465 APOA5 11_70 0 210.32 0.0e+00 6.04
3154 APOA1 11_70 0 246.24 0.0e+00 -17.98
7898 PAFAH1B2 11_70 0 334.67 0.0e+00 21.90
6005 SIDT2 11_70 0 147.00 0.0e+00 -15.38
6006 TAGLN 11_70 0 98.49 0.0e+00 3.91
6785 PCSK7 11_70 0 156.37 3.4e-17 11.85
7745 RNF214 11_70 0 63.46 0.0e+00 -5.11
2466 CEP164 11_70 0 22.37 0.0e+00 1.40
9720 BACE1 11_70 0 64.24 0.0e+00 0.87
4881 FXYD2 11_70 0 16.65 0.0e+00 1.29
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 8_21"
genename region_tag susie_pip mu2 PVE z
5836 CSGALNACT1 8_21 0 45.28 0 8.37
1906 INTS10 8_21 0 15.47 0 -2.61
8739 LPL 8_21 0 416.96 0 17.48
| 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
53977 rs2642420 1_112 1.000 44.74 1.4e-04 -7.64
56631 rs2103827 1_117 1.000 179.19 5.7e-04 19.40
56632 rs11122453 1_117 1.000 360.35 1.1e-03 22.94
57113 rs766167074 1_118 1.000 6750.47 2.2e-02 3.00
69357 rs1042034 2_13 1.000 447.45 1.4e-03 -22.04
69358 rs1801699 2_13 1.000 40.25 1.3e-04 -7.66
71093 rs780093 2_16 1.000 93.81 3.0e-04 -11.04
179767 rs9817452 3_97 1.000 53.69 1.7e-04 7.45
216334 rs7696472 4_48 1.000 73.78 2.4e-04 8.45
227093 rs35518360 4_67 1.000 221.65 7.1e-04 -15.55
227159 rs13140033 4_68 1.000 128.87 4.1e-04 -11.78
271710 rs62369502 5_28 1.000 48.58 1.5e-04 -6.90
364838 rs191555775 6_104 1.000 61.65 2.0e-04 -9.19
410803 rs6961342 7_80 1.000 152.20 4.9e-04 -11.19
425100 rs1402522 8_13 1.000 53.87 1.7e-04 7.83
425861 rs17810889 8_15 1.000 55.67 1.8e-04 7.91
430062 rs75835816 8_21 1.000 438.39 1.4e-03 -21.65
430098 rs11986461 8_21 1.000 295.35 9.4e-04 19.74
491202 rs2777798 9_52 1.000 389.14 1.2e-03 14.36
491210 rs2777804 9_52 1.000 157.08 5.0e-04 4.64
514152 rs71007692 10_28 1.000 6429.39 2.1e-02 -2.78
536281 rs17875416 10_71 1.000 60.30 1.9e-04 6.96
558784 rs7123635 11_28 1.000 50.67 1.6e-04 -7.04
564942 rs695110 11_42 1.000 108.52 3.5e-04 -14.31
601536 rs6581124 12_35 1.000 53.58 1.7e-04 8.14
601555 rs7397189 12_36 1.000 85.02 2.7e-04 11.12
605596 rs2137537 12_44 1.000 41.11 1.3e-04 -5.98
621085 rs61941677 12_76 1.000 141.65 4.5e-04 -14.41
621091 rs12230146 12_76 1.000 61.86 2.0e-04 7.50
673554 rs13379043 14_34 1.000 73.64 2.3e-04 8.46
694679 rs72737411 15_25 1.000 43.97 1.4e-04 -6.21
722747 rs12925793 16_36 1.000 63.78 2.0e-04 8.37
723116 rs2276329 16_37 1.000 50.74 1.6e-04 -6.87
741834 rs55764662 17_26 1.000 188.64 6.0e-04 -13.87
753750 rs146424771 18_3 1.000 67.99 2.2e-04 2.47
766432 rs11082766 18_27 1.000 233.47 7.5e-04 13.66
766452 rs6507938 18_27 1.000 656.73 2.1e-03 32.32
766453 rs118043171 18_27 1.000 691.07 2.2e-03 27.52
766672 rs74461650 18_28 1.000 112.50 3.6e-04 10.92
779926 rs111500536 19_8 1.000 45.49 1.5e-04 6.62
779929 rs116843064 19_8 1.000 227.86 7.3e-04 16.71
780915 rs1865063 19_10 1.000 218.08 7.0e-04 -14.13
780917 rs3745683 19_10 1.000 90.90 2.9e-04 -13.33
787908 rs56287732 19_23 1.000 72.31 2.3e-04 -7.25
787946 rs889140 19_23 1.000 71.65 2.3e-04 8.63
791323 rs62117204 19_32 1.000 94.72 3.0e-04 12.88
791341 rs111794050 19_32 1.000 79.05 2.5e-04 9.35
791374 rs814573 19_32 1.000 729.74 2.3e-03 -29.76
791380 rs4803775 19_32 1.000 185.47 5.9e-04 14.22
807664 rs147591082 20_28 1.000 58.83 1.9e-04 -7.73
861393 rs140584594 1_67 1.000 159.89 5.1e-04 14.42
899551 rs534602560 3_9 1.000 21485.52 6.9e-02 3.60
899628 rs56233336 3_9 1.000 22046.51 7.0e-02 5.05
907149 rs142955295 3_35 1.000 62503.79 2.0e-01 -6.08
923637 rs35945848 4_2 1.000 2312.94 7.4e-03 3.13
927284 rs761645978 4_32 1.000 2036.26 6.5e-03 2.66
927288 rs77240109 4_32 1.000 2030.71 6.5e-03 2.78
931686 rs1611236 6_24 1.000 20314.09 6.5e-02 -4.07
996882 rs4149307 9_53 1.000 493.21 1.6e-03 22.75
997116 rs11789603 9_53 1.000 331.76 1.1e-03 19.85
997195 rs2740488 9_53 1.000 663.11 2.1e-03 -31.28
1022895 rs11601507 11_4 1.000 52.70 1.7e-04 -6.79
1037372 rs57808037 11_37 1.000 20196.22 6.4e-02 4.04
1037377 rs146923372 11_37 1.000 20197.81 6.4e-02 4.02
1041137 rs964184 11_70 1.000 595.31 1.9e-03 6.29
1041363 rs613808 11_70 1.000 1001.00 3.2e-03 -32.58
1051358 rs10507274 12_72 1.000 47.74 1.5e-04 6.73
1058934 rs547584892 12_75 1.000 934.25 3.0e-03 1.19
1074677 rs261290 15_26 1.000 1991.71 6.4e-03 -52.03
1074775 rs12708454 15_26 1.000 749.58 2.4e-03 34.54
1092865 rs193695 16_31 1.000 1101.75 3.5e-03 50.33
1092959 rs5883 16_31 1.000 859.07 2.7e-03 19.54
1092977 rs117427818 16_31 1.000 739.41 2.4e-03 -44.88
1118535 rs11556624 17_23 1.000 46.51 1.5e-04 3.65
1124583 rs67479069 17_27 1.000 15128.84 4.8e-02 -4.07
1139021 rs61371437 19_34 1.000 33893.35 1.1e-01 4.89
1139030 rs113176985 19_34 1.000 33923.30 1.1e-01 4.85
1139033 rs374141296 19_34 1.000 33973.65 1.1e-01 5.19
1144766 rs12975366 19_37 1.000 81.54 2.6e-04 -10.99
1165848 rs780018294 22_10 1.000 1299.38 4.1e-03 -0.71
324689 rs181268076 6_27 0.999 62.98 2.0e-04 -7.07
416795 rs6977416 7_93 0.999 39.45 1.3e-04 -6.32
430048 rs11986942 8_21 0.999 415.62 1.3e-03 32.15
621069 rs3782287 12_76 0.999 62.55 2.0e-04 -11.72
621099 rs11834751 12_76 0.999 50.17 1.6e-04 -4.27
726805 rs12443634 16_45 0.999 70.10 2.2e-04 10.14
766472 rs8093206 18_27 0.999 82.23 2.6e-04 -8.45
780946 rs35753011 19_10 0.999 90.60 2.9e-04 -0.92
6540 rs603412 1_14 0.998 39.44 1.3e-04 -6.26
32964 rs185073199 1_73 0.998 31.45 1.0e-04 5.49
35488 rs4657041 1_79 0.998 33.92 1.1e-04 -5.77
499073 rs115478735 9_70 0.998 92.56 2.9e-04 9.87
885272 rs79800183 2_12 0.998 56.38 1.8e-04 7.16
1074113 rs28690720 15_26 0.998 388.37 1.2e-03 -24.99
722927 rs200561116 16_36 0.997 230.44 7.3e-04 16.19
790273 rs11879413 19_28 0.996 31.46 1.0e-04 5.64
581005 rs4937122 11_77 0.995 32.96 1.0e-04 -5.50
539749 rs10901802 10_78 0.994 34.46 1.1e-04 5.87
682400 rs35007880 14_52 0.994 50.30 1.6e-04 -7.17
652539 rs185932947 13_52 0.993 28.01 8.9e-05 5.10
741577 rs151322266 17_25 0.993 30.50 9.7e-05 -4.95
747647 rs113408695 17_39 0.993 35.59 1.1e-04 -5.86
1069294 rs2494732 14_55 0.993 106.73 3.4e-04 1.84
1074653 rs11071376 15_26 0.993 253.89 8.0e-04 7.19
228914 rs111349657 4_72 0.992 33.21 1.1e-04 -5.71
695700 rs11071771 15_29 0.990 36.55 1.2e-04 -5.86
1042791 rs555766713 11_70 0.990 174.37 5.5e-04 19.20
562564 rs6591179 11_36 0.989 39.39 1.2e-04 6.69
779931 rs62117512 19_8 0.989 57.14 1.8e-04 10.40
624255 rs76734539 12_82 0.985 30.13 9.5e-05 5.68
968560 rs118027010 8_80 0.985 33.09 1.0e-04 4.89
142885 rs4681065 3_19 0.984 29.09 9.1e-05 5.17
579487 rs1219430 11_74 0.983 30.99 9.7e-05 -5.72
766468 rs62101781 18_27 0.982 285.64 8.9e-04 19.53
1041340 rs12721030 11_70 0.978 656.64 2.0e-03 31.91
1092885 rs183130 16_31 0.978 3552.25 1.1e-02 74.18
700253 rs16972386 15_38 0.977 27.73 8.6e-05 -4.97
752053 rs72854483 17_46 0.977 26.60 8.3e-05 -4.92
168989 rs334563 3_74 0.973 46.62 1.4e-04 6.67
323549 rs3095311 6_25 0.973 89.45 2.8e-04 -10.67
369119 rs11971790 7_2 0.973 59.63 1.9e-04 -6.93
738875 rs2227735 17_17 0.971 84.49 2.6e-04 -9.43
383910 rs9490 7_28 0.970 27.19 8.4e-05 4.97
574069 rs72980276 11_59 0.969 26.96 8.3e-05 -4.99
558964 rs7111419 11_29 0.968 30.65 9.5e-05 5.80
379665 rs75348547 7_22 0.965 28.06 8.6e-05 -4.86
430044 rs6999813 8_21 0.965 447.06 1.4e-03 31.97
474552 rs145804707 9_18 0.963 25.20 7.7e-05 -4.69
90786 rs9248 2_54 0.961 39.34 1.2e-04 6.29
1125445 rs113667149 17_27 0.961 15123.01 4.6e-02 -4.35
323170 rs1233480 6_23 0.960 98.00 3.0e-04 -9.65
553109 rs12288512 11_19 0.960 54.21 1.7e-04 -7.27
179525 rs816547 3_97 0.959 25.59 7.8e-05 -4.76
324393 rs77424484 6_25 0.958 57.74 1.8e-04 -7.74
562360 rs12294913 11_36 0.957 30.95 9.4e-05 -5.80
273625 rs113088001 5_31 0.956 33.18 1.0e-04 -5.21
726812 rs11641142 16_45 0.955 48.13 1.5e-04 9.01
325605 rs58737173 6_28 0.947 38.01 1.1e-04 5.14
327948 rs142449754 6_32 0.947 32.84 9.9e-05 -6.44
194053 rs13116176 4_4 0.943 29.91 9.0e-05 -4.46
813570 rs759884584 20_38 0.940 31.65 9.5e-05 2.74
322926 rs3132390 6_22 0.938 90.29 2.7e-04 -9.77
1074923 rs2070895 15_26 0.934 2450.05 7.3e-03 50.25
454017 rs2570952 8_69 0.931 30.04 8.9e-05 5.25
274561 rs173964 5_33 0.930 82.21 2.4e-04 -7.90
379842 rs2699814 7_23 0.929 27.06 8.0e-05 4.98
769328 rs41292412 18_31 0.928 30.95 9.2e-05 -5.49
370812 rs10262715 7_6 0.927 25.52 7.5e-05 4.82
303664 rs4958365 5_90 0.925 31.50 9.3e-05 4.70
70800 rs577641882 2_15 0.924 33.29 9.8e-05 -5.71
429836 rs113231830 8_20 0.924 26.24 7.7e-05 -4.89
293520 rs11064 5_72 0.921 25.97 7.6e-05 4.83
55141 rs12132342 1_115 0.919 24.67 7.2e-05 -4.53
829382 rs12321 22_9 0.919 25.39 7.4e-05 4.48
601491 rs1874888 12_35 0.918 34.16 1.0e-04 -6.19
53975 rs883847 1_112 0.914 33.40 9.7e-05 6.84
72448 rs11124265 2_20 0.912 25.56 7.4e-05 4.48
328700 rs581465 6_34 0.908 25.98 7.5e-05 4.83
732166 rs3760230 17_3 0.908 33.44 9.7e-05 -5.67
832665 rs531420711 22_17 0.899 25.77 7.4e-05 -4.81
361516 rs6905582 6_99 0.897 24.41 7.0e-05 -4.49
521138 rs10761739 10_42 0.897 41.77 1.2e-04 6.40
590471 rs964974 12_15 0.896 37.37 1.1e-04 -5.99
1035543 rs4930352 11_37 0.893 621.13 1.8e-03 7.80
791037 rs73036721 19_30 0.892 23.98 6.8e-05 4.35
621062 rs838876 12_76 0.889 127.63 3.6e-04 -13.65
354053 rs112120898 6_84 0.881 24.36 6.8e-05 4.42
81544 rs753446737 2_36 0.880 24.72 6.9e-05 4.43
546948 rs2923096 11_8 0.879 33.78 9.5e-05 5.90
808141 rs11699481 20_28 0.876 31.75 8.9e-05 -6.33
601578 rs140734681 12_36 0.875 23.58 6.6e-05 -1.53
780891 rs4423524 19_9 0.870 37.81 1.0e-04 7.81
56624 rs6678475 1_117 0.863 26.49 7.3e-05 -1.21
672015 rs194749 14_32 0.860 42.17 1.2e-04 6.09
438087 rs62515418 8_40 0.859 24.57 6.7e-05 4.31
904195 rs115787626 3_33 0.859 55.72 1.5e-04 -7.73
739210 rs2066713 17_18 0.856 35.87 9.8e-05 -5.77
285860 rs538659275 5_57 0.855 27.28 7.4e-05 4.87
457764 rs10095930 8_78 0.855 77.37 2.1e-04 5.62
794288 rs34637812 19_38 0.855 25.85 7.0e-05 -4.61
883409 rs144977853 2_12 0.853 36.30 9.9e-05 3.35
296608 rs10057561 5_77 0.852 26.97 7.3e-05 -4.91
1041384 rs141469619 11_70 0.852 111.38 3.0e-04 -12.44
605717 rs1707498 12_44 0.850 29.84 8.1e-05 5.11
685364 rs11161243 15_4 0.850 32.74 8.9e-05 5.48
324560 rs9276189 6_27 0.842 38.49 1.0e-04 3.68
737309 rs78792395 17_15 0.841 26.24 7.0e-05 4.54
795594 rs2316866 20_1 0.841 25.44 6.8e-05 -4.54
498011 rs111472765 9_67 0.839 24.28 6.5e-05 4.38
1048445 rs7305678 12_7 0.835 43.34 1.2e-04 -6.39
323810 rs2524096 6_25 0.828 50.71 1.3e-04 9.03
512611 rs145407036 10_24 0.823 29.14 7.6e-05 -5.44
111561 rs1862069 2_102 0.819 30.19 7.9e-05 5.12
148231 rs11709009 3_30 0.813 23.79 6.2e-05 -4.09
367574 rs34207042 6_110 0.813 24.85 6.5e-05 4.28
512834 rs12765967 10_25 0.811 29.93 7.7e-05 -5.17
211154 rs723585 4_40 0.809 36.95 9.5e-05 5.82
202749 rs56147366 4_22 0.807 44.94 1.2e-04 -6.56
280876 rs3733890 5_46 0.802 25.50 6.5e-05 -4.52SNPs 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
907149 rs142955295 3_35 1.000 62503.79 2.0e-01 -6.08
907115 rs9853458 3_35 0.350 62398.81 7.0e-02 6.07
907113 rs9876508 3_35 0.251 62398.61 5.0e-02 6.07
907114 rs9815766 3_35 0.135 62395.93 2.7e-02 6.06
907086 rs1049256 3_35 0.087 62395.02 1.7e-02 6.06
907083 rs7634902 3_35 0.071 62394.92 1.4e-02 6.06
907080 rs3811696 3_35 0.042 62392.27 8.4e-03 6.06
907081 rs3811695 3_35 0.017 62391.83 3.4e-03 6.05
907079 rs4855850 3_35 0.031 62390.77 6.2e-03 6.06
907116 rs7374277 3_35 0.153 62389.89 3.0e-02 6.07
907174 rs34451146 3_35 0.056 62388.71 1.1e-02 -6.07
907187 rs9814765 3_35 0.044 62388.59 8.7e-03 -6.07
907188 rs11130221 3_35 0.044 62388.59 8.7e-03 -6.07
907194 rs13063621 3_35 0.040 62388.55 8.0e-03 -6.07
907203 rs9871654 3_35 0.032 62388.43 6.3e-03 -6.06
907117 rs7374183 3_35 0.084 62387.21 1.7e-02 6.07
907142 rs7634886 3_35 0.048 62386.44 9.6e-03 -6.07
907175 rs57648519 3_35 0.029 62386.01 5.8e-03 -6.07
907165 rs6446295 3_35 0.003 62384.80 5.0e-04 -6.04
907065 rs3749240 3_35 0.108 62384.42 2.2e-02 6.09
907072 rs34614773 3_35 0.049 62382.94 9.7e-03 6.08
907160 rs7431106 3_35 0.001 62382.15 2.8e-04 -6.04
907146 rs9865480 3_35 0.010 62381.97 2.0e-03 -6.05
907039 rs1491986 3_35 0.080 62379.28 1.6e-02 6.09
907151 rs6809431 3_35 0.003 62379.04 6.6e-04 -6.05
907147 rs60205400 3_35 0.003 62379.00 6.3e-04 -6.05
907071 rs11130219 3_35 0.011 62378.90 2.3e-03 6.07
907169 rs6766836 3_35 0.003 62378.67 6.5e-04 -6.05
907154 rs9859153 3_35 0.003 62378.65 6.3e-04 -6.05
907144 rs9882639 3_35 0.003 62378.32 6.5e-04 -6.05
907070 rs11130218 3_35 0.006 62375.60 1.1e-03 6.07
907056 rs12381242 3_35 0.167 62374.25 3.3e-02 6.10
907089 rs10632976 3_35 0.017 62374.01 3.4e-03 6.06
907128 rs7372730 3_35 0.124 62373.14 2.5e-02 6.09
907132 rs9855505 3_35 0.119 62373.12 2.4e-02 6.09
907123 rs7429353 3_35 0.027 62370.00 5.3e-03 6.08
907127 rs7372725 3_35 0.025 62369.98 5.1e-03 6.08
907048 rs11709680 3_35 0.017 62367.31 3.3e-03 6.09
907047 rs11716575 3_35 0.017 62367.28 3.3e-03 6.09
907059 rs4855862 3_35 0.007 62366.46 1.4e-03 6.08
907034 rs6785549 3_35 0.187 62364.86 3.7e-02 6.12
907046 rs3749241 3_35 0.014 62362.22 2.8e-03 6.08
907130 rs9872864 3_35 0.010 62362.20 1.9e-03 6.08
907051 rs4855841 3_35 0.001 62356.06 2.2e-04 6.07
907140 rs12490656 3_35 0.294 62347.55 5.8e-02 -6.14
907135 rs35365539 3_35 0.000 62336.66 7.1e-05 6.05
907121 rs7372966 3_35 0.000 62331.01 3.6e-05 6.05
907030 rs9873183 3_35 0.001 62329.27 2.2e-04 6.09
907124 rs7426497 3_35 0.000 62328.94 5.5e-05 6.06
907053 rs4855867 3_35 0.000 62318.58 5.3e-08 6.01SNPs 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
907149 rs142955295 3_35 1.000 62503.79 0.2000 -6.08
1139021 rs61371437 19_34 1.000 33893.35 0.1100 4.89
1139030 rs113176985 19_34 1.000 33923.30 0.1100 4.85
1139033 rs374141296 19_34 1.000 33973.65 0.1100 5.19
899628 rs56233336 3_9 1.000 22046.51 0.0700 5.05
907115 rs9853458 3_35 0.350 62398.81 0.0700 6.07
899551 rs534602560 3_9 1.000 21485.52 0.0690 3.60
931686 rs1611236 6_24 1.000 20314.09 0.0650 -4.07
1037372 rs57808037 11_37 1.000 20196.22 0.0640 4.04
1037377 rs146923372 11_37 1.000 20197.81 0.0640 4.02
907140 rs12490656 3_35 0.294 62347.55 0.0580 -6.14
907113 rs9876508 3_35 0.251 62398.61 0.0500 6.07
1124583 rs67479069 17_27 1.000 15128.84 0.0480 -4.07
1125445 rs113667149 17_27 0.961 15123.01 0.0460 -4.35
899613 rs17035943 3_9 0.560 21886.54 0.0390 4.04
907034 rs6785549 3_35 0.187 62364.86 0.0370 6.12
907056 rs12381242 3_35 0.167 62374.25 0.0330 6.10
899625 rs28897670 3_9 0.432 21858.09 0.0300 4.04
907116 rs7374277 3_35 0.153 62389.89 0.0300 6.07
907114 rs9815766 3_35 0.135 62395.93 0.0270 6.06
899620 rs73813138 3_9 0.370 21858.26 0.0260 4.04
907128 rs7372730 3_35 0.124 62373.14 0.0250 6.09
907132 rs9855505 3_35 0.119 62373.12 0.0240 6.09
57113 rs766167074 1_118 1.000 6750.47 0.0220 3.00
907065 rs3749240 3_35 0.108 62384.42 0.0220 6.09
514152 rs71007692 10_28 1.000 6429.39 0.0210 -2.78
1124549 rs111164082 17_27 0.390 15112.66 0.0190 -4.28
899623 rs28897669 3_9 0.261 21857.09 0.0180 4.04
907086 rs1049256 3_35 0.087 62395.02 0.0170 6.06
907117 rs7374183 3_35 0.084 62387.21 0.0170 6.07
907039 rs1491986 3_35 0.080 62379.28 0.0160 6.09
931655 rs1633020 6_24 0.242 20303.97 0.0160 -4.09
899619 rs35839040 3_9 0.211 21858.27 0.0150 4.03
907083 rs7634902 3_35 0.071 62394.92 0.0140 6.06
931617 rs2844838 6_24 0.182 20305.63 0.0120 -4.08
907174 rs34451146 3_35 0.056 62388.71 0.0110 -6.07
931672 rs1611228 6_24 0.172 20305.67 0.0110 -4.08
1092885 rs183130 16_31 0.978 3552.25 0.0110 74.18
514158 rs2472183 10_28 0.489 6464.38 0.0100 -2.84
514161 rs11011452 10_28 0.478 6464.51 0.0099 -2.83
514151 rs2474565 10_28 0.473 6464.33 0.0098 -2.84
907072 rs34614773 3_35 0.049 62382.94 0.0097 6.08
931604 rs1633033 6_24 0.149 20305.68 0.0097 -4.08
907142 rs7634886 3_35 0.048 62386.44 0.0096 -6.07
907187 rs9814765 3_35 0.044 62388.59 0.0087 -6.07
907188 rs11130221 3_35 0.044 62388.59 0.0087 -6.07
907080 rs3811696 3_35 0.042 62392.27 0.0084 6.06
931659 rs1633018 6_24 0.125 20303.52 0.0081 -4.08
907194 rs13063621 3_35 0.040 62388.55 0.0080 -6.07
57110 rs10489611 1_118 0.358 6707.80 0.0077 3.25SNPs 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
1092885 rs183130 16_31 0.978 3552.25 1.1e-02 74.18
1092897 rs3764261 16_31 0.007 3542.86 7.5e-05 74.09
1092899 rs821840 16_31 0.016 3544.17 1.8e-04 74.05
1092882 rs247617 16_31 0.000 3528.10 4.0e-07 74.01
1092907 rs17231506 16_31 0.000 3533.67 1.2e-06 73.99
1092906 rs36229491 16_31 0.000 3524.21 3.3e-08 73.95
1092879 rs247616 16_31 0.000 3530.92 3.0e-07 73.94
1092872 rs12446515 16_31 0.000 3484.88 1.1e-15 73.66
1092873 rs56156922 16_31 0.000 3492.97 1.1e-14 73.56
1092895 rs12149545 16_31 0.000 3443.09 0.0e+00 72.85
1092874 rs56228609 16_31 0.000 3423.52 0.0e+00 72.66
1092875 rs173539 16_31 0.000 3450.29 0.0e+00 72.51
1092913 rs1800775 16_31 0.000 2559.22 0.0e+00 68.24
1092915 rs3816117 16_31 0.000 2522.94 0.0e+00 68.09
1092950 rs1532625 16_31 0.000 2874.85 0.0e+00 67.51
1092949 rs7205804 16_31 0.000 2859.05 0.0e+00 67.43
1092916 rs711752 16_31 0.000 2928.27 1.2e-15 67.42
1092917 rs708272 16_31 0.000 2921.25 8.3e-16 67.35
1092951 rs1532624 16_31 0.000 2853.85 0.0e+00 67.28
1092918 rs34620476 16_31 0.000 2884.29 4.9e-17 67.13
1092927 rs11508026 16_31 0.000 2828.29 1.0e-18 66.59
1092925 rs12720926 16_31 0.000 2825.42 1.0e-18 66.58
1092933 rs4784741 16_31 0.000 2805.82 0.0e+00 66.40
1092936 rs12444012 16_31 0.000 2805.41 0.0e+00 66.39
1092866 rs72786786 16_31 0.000 2738.46 0.0e+00 65.91
1092930 rs8045855 16_31 0.359 1217.80 1.4e-03 -64.00
1092931 rs12720922 16_31 0.048 1208.25 1.8e-04 -63.89
1092955 rs11076175 16_31 0.374 1202.03 1.4e-03 -63.88
1092956 rs7499892 16_31 0.199 1202.28 7.6e-04 -63.81
1092919 rs1864163 16_31 0.000 1166.30 0.0e+00 -63.79
1092934 rs12720908 16_31 0.020 1199.32 7.8e-05 -63.63
1092947 rs9939224 16_31 0.000 1257.69 2.1e-16 62.90
1092920 rs5817082 16_31 0.000 1150.56 0.0e+00 -62.71
1092926 rs7203984 16_31 0.000 1153.59 1.0e-06 -60.89
1092957 rs289713 16_31 0.000 1118.87 3.5e-09 60.66
1092932 rs118146573 16_31 0.000 880.05 0.0e+00 -55.94
1092881 rs12923459 16_31 0.000 1298.78 0.0e+00 -54.09
1092902 rs711751 16_31 0.000 1757.12 0.0e+00 53.65
1092960 rs11076176 16_31 0.000 1156.22 0.0e+00 -53.07
1092946 rs9926440 16_31 0.000 1437.37 0.0e+00 52.47
1092870 rs7203286 16_31 0.000 1200.93 0.0e+00 -52.12
1074677 rs261290 15_26 1.000 1991.71 6.4e-03 -52.03
1092864 rs9989419 16_31 0.000 1176.35 2.1e-07 51.57
1074684 rs7350789 15_26 0.000 1957.87 1.4e-13 51.52
1074687 rs7177289 15_26 0.000 1955.50 3.7e-14 51.50
1074686 rs261291 15_26 0.000 1953.61 5.3e-15 51.48
1092867 rs12448528 16_31 0.000 1282.14 5.6e-10 51.18
1074690 rs35853021 15_26 0.000 1890.18 0.0e+00 51.03
1074693 rs2043085 15_26 0.000 1776.25 0.0e+00 -50.93
1074698 rs753768034 15_26 0.000 1766.65 0.0e+00 -50.87Gene 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] 33if (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 negative regulation of cholesterol storage (GO:0010887)
2 cholesterol homeostasis (GO:0042632)
3 sterol homeostasis (GO:0055092)
4 regulation of lipase activity (GO:0060191)
5 regulation of cholesterol storage (GO:0010885)
6 negative regulation of lipid storage (GO:0010888)
7 regulation of lipoprotein lipase activity (GO:0051004)
8 regulation of DNA damage response, signal transduction by p53 class mediator (GO:0043516)
9 positive regulation of telomere maintenance (GO:0032206)
Overlap Adjusted.P.value Genes
1 2/10 0.02099954 ABCA1;TTC39B
2 3/71 0.02099954 ABCA1;GPIHBP1;TTC39B
3 3/72 0.02099954 ABCA1;GPIHBP1;TTC39B
4 2/14 0.02099954 PCSK6;GPIHBP1
5 2/16 0.02210770 ABCA1;TTC39B
6 2/20 0.02749245 ABCA1;TTC39B
7 2/21 0.02749245 PCSK6;GPIHBP1
8 2/28 0.04298914 PTTG1IP;KMT5A
9 2/32 0.04993501 TNKS;GNL3
[1] "GO_Cellular_Component_2021"
[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)RP11-346C20.3 gene(s) from the input list not found in DisGeNET CURATEDRP11-219B17.3 gene(s) from the input list not found in DisGeNET CURATEDIRF9 gene(s) from the input list not found in DisGeNET CURATEDALLC gene(s) from the input list not found in DisGeNET CURATEDRP11-10A14.4 gene(s) from the input list not found in DisGeNET CURATEDGNL3 gene(s) from the input list not found in DisGeNET CURATEDNTAN1 gene(s) from the input list not found in DisGeNET CURATEDUBE2K gene(s) from the input list not found in DisGeNET CURATEDTTC39B gene(s) from the input list not found in DisGeNET CURATEDPCSK6 gene(s) from the input list not found in DisGeNET CURATEDPSRC1 gene(s) from the input list not found in DisGeNET CURATEDBEND3 gene(s) from the input list not found in DisGeNET CURATEDKMT5A gene(s) from the input list not found in DisGeNET CURATEDTNKS gene(s) from the input list not found in DisGeNET CURATEDRP11-136O12.2 gene(s) from the input list not found in DisGeNET CURATEDAKNA gene(s) from the input list not found in DisGeNET CURATEDSNAPC2 gene(s) from the input list not found in DisGeNET CURATEDRP4-781K5.7 gene(s) from the input list not found in DisGeNET CURATEDPTTG1IP gene(s) from the input list not found in DisGeNET CURATED Description FDR Ratio BgRatio
10 Hypercholesterolemia 0.005875077 2/14 39/9703
14 Mucopolysaccharidosis I 0.005875077 1/14 1/9703
16 Mucopolysaccharidosis V 0.005875077 1/14 1/9703
22 Tangier Disease 0.005875077 1/14 1/9703
23 Hurler-Scheie Syndrome 0.005875077 1/14 1/9703
25 Pfaundler-Hurler Syndrome 0.005875077 1/14 1/9703
34 Symmetrical dyschromatosis of extremities 0.005875077 1/14 1/9703
35 Hypoalphalipoproteinemias 0.005875077 1/14 1/9703
38 Tangier Disease Neuropathy 0.005875077 1/14 1/9703
46 alpha-L-Iduronidase Deficiency 0.005875077 1/14 1/9703******************************************
* *
* Welcome to WebGestaltR ! *
* *
******************************************Loading the functional categories...
Loading the ID list...
Loading the reference list...
Performing the enrichment analysis...Warning in oraEnrichment(interestGeneList, referenceGeneList, geneSet,
minNum = minNum, : No significant gene set is identified based on FDR 0.05!NULL
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