Phosphate (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
Modified: analysis/ukb-d-30700_irnt_Liver.Rmd
Modified: analysis/ukb-d-30710_irnt_Liver.Rmd
Modified: analysis/ukb-d-30720_irnt_Liver.Rmd
Modified: analysis/ukb-d-30730_irnt_Liver.Rmd
Modified: analysis/ukb-d-30740_irnt_Liver.Rmd
Modified: analysis/ukb-d-30750_irnt_Liver.Rmd
Modified: analysis/ukb-d-30760_irnt_Liver.Rmd
Modified: analysis/ukb-d-30770_irnt_Liver.Rmd
Modified: analysis/ukb-d-30780_irnt_Liver.Rmd
Modified: analysis/ukb-d-30790_irnt_Liver.Rmd
Modified: analysis/ukb-d-30800_irnt_Liver.Rmd
Modified: analysis/ukb-d-30810_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-30810_irnt_Liver.Rmd) and HTML (docs/ukb-d-30810_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 Phosphate (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-30810_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.0194001836 0.0001642049 #estimated group prior variance
estimated_group_prior_var <- group_prior_var_rec[,ncol(group_prior_var_rec)]
names(estimated_group_prior_var) <- c("gene", "snp")
print(estimated_group_prior_var) gene snp
9.275474 15.211986 #report sample size
print(sample_size)[1] 314658#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.006234039 0.069042946 #compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))[1] 0.03140854 2.20274272Genes 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
5788 RSPO3 6_84 0.997 30.08 9.5e-05 6.17
721 WIPI1 17_39 0.995 36.90 1.2e-04 7.09
18 TMEM176A 7_93 0.993 31.50 9.9e-05 -5.81
3212 CCND2 12_4 0.988 27.07 8.5e-05 -5.31
12467 RP11-219B17.3 15_27 0.984 33.97 1.1e-04 -6.69
10050 SUPT3H 6_34 0.983 30.17 9.4e-05 7.53
9855 PALM3 19_11 0.979 39.85 1.2e-04 6.53
4454 GRHL1 2_6 0.972 36.46 1.1e-04 -6.46
10338 PRIM1 12_35 0.963 21.14 6.5e-05 -4.34
394 WDR37 10_2 0.944 19.82 5.9e-05 -4.13
2474 CBL 11_71 0.923 28.23 8.3e-05 -5.31
5373 AKT1 14_55 0.921 19.40 5.7e-05 -4.03
1945 TRMT1 19_10 0.920 22.51 6.6e-05 -4.55
5415 SYTL1 1_19 0.895 21.95 6.2e-05 -4.58
7353 CHMP4C 8_58 0.888 45.37 1.3e-04 -7.00
11538 MIATNB 22_8 0.882 18.48 5.2e-05 -3.90
2887 NRBP1 2_16 0.873 20.57 5.7e-05 4.57
3965 ICE2 15_27 0.861 22.70 6.2e-05 5.51
12242 RP11-1109F11.3 12_54 0.857 20.26 5.5e-05 4.03
8173 ELP5 17_6 0.843 28.86 7.7e-05 -5.28
2436 C11orf63 11_75 0.820 19.31 5.0e-05 3.97
5822 GIGYF1 7_62 0.809 20.43 5.3e-05 -4.34
4701 ZDHHC4 7_9 0.804 26.03 6.7e-05 -5.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
11381 LRRC37A2 17_27 0.00 44009.46 0.00000 -9.14
9663 ARL17A 17_27 0.00 29599.47 0.00000 -7.09
11199 LINC00271 6_89 0.00 25642.10 0.00000 -5.85
4604 AHI1 6_89 0.00 8890.73 0.00000 -4.30
3310 KANSL1 17_27 0.00 7065.31 0.00000 3.14
12113 RP11-798G7.6 17_27 0.00 7065.31 0.00000 3.14
8846 LRRC37A 17_27 0.00 4692.35 0.00000 2.50
802 NSF 17_27 0.00 4547.72 0.00000 3.26
9773 MAPT 17_27 0.00 2144.27 0.00000 1.98
12583 AC142472.6 17_27 0.00 1188.58 0.00000 -1.04
6873 IP6K3 6_28 0.07 1065.01 0.00024 36.39
6678 ARHGAP27 17_27 0.00 964.23 0.00000 -1.89
2310 GOSR2 17_27 0.00 950.07 0.00000 -1.97
8499 DCAKD 17_27 0.00 949.31 0.00000 -2.51
5418 NBPF3 1_15 0.00 521.33 0.00000 -22.48
2661 HBS1L 6_89 0.00 461.96 0.00000 -2.51
3183 ALDH8A1 6_89 0.00 443.09 0.00000 -0.46
2301 WNT3 17_27 0.00 310.83 0.00000 -1.03
10511 TBKBP1 17_27 0.00 271.39 0.00000 -7.06
2309 KPNB1 17_27 0.00 222.35 0.00000 -6.91Genes 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
6873 IP6K3 6_28 0.070 1065.01 2.4e-04 36.39
2338 SYNGR2 17_44 0.768 53.37 1.3e-04 7.88
7353 CHMP4C 8_58 0.888 45.37 1.3e-04 -7.00
721 WIPI1 17_39 0.995 36.90 1.2e-04 7.09
9855 PALM3 19_11 0.979 39.85 1.2e-04 6.53
3411 CSTA 3_76 0.277 130.33 1.1e-04 17.53
8545 EHBP1L1 11_36 0.668 53.24 1.1e-04 -8.48
4454 GRHL1 2_6 0.972 36.46 1.1e-04 -6.46
12467 RP11-219B17.3 15_27 0.984 33.97 1.1e-04 -6.69
18 TMEM176A 7_93 0.993 31.50 9.9e-05 -5.81
5788 RSPO3 6_84 0.997 30.08 9.5e-05 6.17
10050 SUPT3H 6_34 0.983 30.17 9.4e-05 7.53
3212 CCND2 12_4 0.988 27.07 8.5e-05 -5.31
2474 CBL 11_71 0.923 28.23 8.3e-05 -5.31
12035 GTF2I 7_48 0.766 32.72 8.0e-05 -5.76
8173 ELP5 17_6 0.843 28.86 7.7e-05 -5.28
7355 BRI3 7_60 0.746 31.60 7.5e-05 -5.98
6133 MFSD6 2_113 0.478 45.06 6.8e-05 -6.91
4701 ZDHHC4 7_9 0.804 26.03 6.7e-05 -5.59
1945 TRMT1 19_10 0.920 22.51 6.6e-05 -4.55Genes 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
6873 IP6K3 6_28 0.070 1065.01 2.4e-04 36.39
5418 NBPF3 1_15 0.000 521.33 0.0e+00 -22.48
3411 CSTA 3_76 0.277 130.33 1.1e-04 17.53
8037 LMAN2 5_106 0.003 86.38 7.0e-07 -14.31
4999 PARP9 3_76 0.000 57.47 7.6e-09 -11.81
11994 MAFTRR 16_44 0.037 117.40 1.4e-05 11.59
3731 MED1 17_23 0.020 91.87 5.7e-06 -10.80
2297 FBXL20 17_23 0.050 91.89 1.4e-05 10.58
2794 KPNA1 3_76 0.000 48.65 2.8e-10 10.24
11381 LRRC37A2 17_27 0.000 44009.46 0.0e+00 -9.14
8545 EHBP1L1 11_36 0.668 53.24 1.1e-04 -8.48
9030 FAM219B 15_35 0.131 51.46 2.1e-05 8.10
8040 THBS3 1_76 0.177 53.22 3.0e-05 -7.95
2338 SYNGR2 17_44 0.768 53.37 1.3e-04 7.88
5212 SCAMP2 15_35 0.060 42.78 8.1e-06 7.74
6849 PGAP3 17_23 0.019 45.02 2.7e-06 -7.72
9025 RPP25 15_35 0.113 46.33 1.7e-05 -7.72
10157 PDLIM7 5_106 0.290 38.82 3.6e-05 7.64
9036 MPI 15_35 0.072 43.04 9.8e-06 7.62
10050 SUPT3H 6_34 0.983 30.17 9.4e-05 7.53Comparing 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.01119163#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
6873 IP6K3 6_28 0.070 1065.01 2.4e-04 36.39
5418 NBPF3 1_15 0.000 521.33 0.0e+00 -22.48
3411 CSTA 3_76 0.277 130.33 1.1e-04 17.53
8037 LMAN2 5_106 0.003 86.38 7.0e-07 -14.31
4999 PARP9 3_76 0.000 57.47 7.6e-09 -11.81
11994 MAFTRR 16_44 0.037 117.40 1.4e-05 11.59
3731 MED1 17_23 0.020 91.87 5.7e-06 -10.80
2297 FBXL20 17_23 0.050 91.89 1.4e-05 10.58
2794 KPNA1 3_76 0.000 48.65 2.8e-10 10.24
11381 LRRC37A2 17_27 0.000 44009.46 0.0e+00 -9.14
8545 EHBP1L1 11_36 0.668 53.24 1.1e-04 -8.48
9030 FAM219B 15_35 0.131 51.46 2.1e-05 8.10
8040 THBS3 1_76 0.177 53.22 3.0e-05 -7.95
2338 SYNGR2 17_44 0.768 53.37 1.3e-04 7.88
5212 SCAMP2 15_35 0.060 42.78 8.1e-06 7.74
6849 PGAP3 17_23 0.019 45.02 2.7e-06 -7.72
9025 RPP25 15_35 0.113 46.33 1.7e-05 -7.72
10157 PDLIM7 5_106 0.290 38.82 3.6e-05 7.64
9036 MPI 15_35 0.072 43.04 9.8e-06 7.62
10050 SUPT3H 6_34 0.983 30.17 9.4e-05 7.53Locus 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: 6_28"
genename region_tag susie_pip mu2 PVE z
11210 RPS18 6_28 0.000 11.37 1.8e-09 0.11
11018 VPS52 6_28 0.000 5.59 3.9e-10 -0.62
11111 WDR46 6_28 0.000 5.77 5.0e-10 0.77
11218 TAPBP 6_28 0.001 38.35 1.5e-07 -6.08
11330 ZBTB22 6_28 0.000 30.89 3.1e-08 5.35
10581 DAXX 6_28 0.001 35.40 1.2e-07 -5.78
2675 CUTA 6_28 0.000 14.34 2.0e-09 -0.45
1340 ITPR3 6_28 0.000 19.19 1.2e-08 -0.77
4831 UQCC2 6_28 0.001 27.34 6.7e-08 -1.69
6873 IP6K3 6_28 0.070 1065.01 2.4e-04 36.39
6874 LEMD2 6_28 0.000 45.28 5.4e-09 -6.89
12314 NUDT3 6_28 0.000 13.77 1.1e-09 2.08
3651 RPS10 6_28 0.000 10.32 1.4e-09 -0.99
3656 SPDEF 6_28 0.000 5.48 3.8e-10 0.73
10147 C6orf106 6_28 0.000 13.24 6.3e-09 2.01
3639 SNRPC 6_28 0.000 10.25 1.6e-09 0.67
586 UHRF1BP1 6_28 0.000 8.84 1.2e-09 -2.35
580 TAF11 6_28 0.000 9.78 8.2e-10 -1.72
581 ANKS1A 6_28 0.000 4.58 3.0e-10 0.32
583 ZNF76 6_28 0.000 16.61 3.6e-09 0.98
282 DEF6 6_28 0.000 8.08 8.2e-10 0.79
2621 PPARD 6_28 0.000 15.70 2.3e-09 -1.81
2622 FANCE 6_28 0.000 10.62 1.0e-09 -1.33
10454 RPL10A 6_28 0.000 6.79 5.4e-10 0.88
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 1_15"
genename region_tag susie_pip mu2 PVE z
5418 NBPF3 1_15 0 521.33 0 -22.48
1235 USP48 1_15 0 7.38 0 -0.82
9856 LDLRAD2 1_15 0 69.50 0 5.06
5419 HSPG2 1_15 0 11.92 0 3.00
5417 CELA3A 1_15 0 10.18 0 -1.68
10971 LINC00339 1_15 0 27.27 0 4.03
735 CDC42 1_15 0 5.92 0 0.84
6947 WNT4 1_15 0 34.11 0 -1.46
9541 ZBTB40 1_15 0 4.57 0 0.19
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 3_76"
genename region_tag susie_pip mu2 PVE z
3411 CSTA 3_76 0.277 130.33 1.1e-04 17.53
2792 FAM162A 3_76 0.000 19.47 5.4e-10 -1.47
6734 CCDC58 3_76 0.000 27.01 8.9e-10 -6.39
10167 WDR5B 3_76 0.000 22.12 2.0e-10 3.89
2794 KPNA1 3_76 0.000 48.65 2.8e-10 10.24
4999 PARP9 3_76 0.000 57.47 7.6e-09 -11.81
8518 PARP15 3_76 0.000 20.26 1.5e-10 2.02
8517 PARP14 3_76 0.000 9.57 5.2e-11 0.67
8023 HSPBAP1 3_76 0.000 6.34 2.6e-11 -1.57
4996 DIRC2 3_76 0.000 6.23 2.7e-11 -1.78
12407 LINC02035 3_76 0.000 8.10 3.0e-11 -0.66
1008 SEMA5B 3_76 0.000 35.34 1.8e-09 4.31
3410 SEC22A 3_76 0.000 4.84 1.8e-11 0.33
10775 HACD2 3_76 0.000 6.02 2.8e-11 -0.75
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 5_106"
genename region_tag susie_pip mu2 PVE z
8146 SIMC1 5_106 0.004 9.47 1.3e-07 -0.89
3450 KIAA1191 5_106 0.004 10.08 1.4e-07 -1.02
8737 ARL10 5_106 0.002 4.57 3.3e-08 -0.05
5758 HIGD2A 5_106 0.002 4.55 3.3e-08 0.07
8738 CLTB 5_106 0.003 5.33 4.6e-08 -0.09
8046 GPRIN1 5_106 0.003 5.93 5.3e-08 -0.16
403 TSPAN17 5_106 0.004 7.99 8.9e-08 0.68
2780 UNC5A 5_106 0.003 9.09 9.3e-08 2.22
6811 HK3 5_106 0.002 4.87 3.8e-08 -0.01
1119 UIMC1 5_106 0.003 6.25 5.3e-08 -1.03
2779 ZNF346 5_106 0.003 6.24 5.7e-08 -0.80
6807 FGFR4 5_106 0.115 20.20 7.4e-06 2.86
7484 NSD1 5_106 0.003 7.66 7.8e-08 -0.79
8039 PRELID1 5_106 0.003 13.27 1.4e-07 -3.36
10820 MXD3 5_106 0.002 5.30 3.9e-08 1.29
8037 LMAN2 5_106 0.003 86.38 7.0e-07 -14.31
10107 PFN3 5_106 0.003 20.21 2.1e-07 0.96
4159 F12 5_106 0.003 29.70 2.5e-07 -5.68
4160 PRR7 5_106 0.005 12.55 2.1e-07 0.45
2778 DBN1 5_106 0.006 13.36 2.4e-07 0.14
10157 PDLIM7 5_106 0.290 38.82 3.6e-05 7.64
12333 RP11-1277A3.3 5_106 0.146 23.60 1.1e-05 -4.54
301 B4GALT7 5_106 0.015 12.79 6.1e-07 -2.39
8144 FAM153A 5_106 0.004 6.60 7.8e-08 0.73
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 16_44"
genename region_tag susie_pip mu2 PVE z
9016 MAF 16_44 0.013 5.21 2.1e-07 0.86
11994 MAFTRR 16_44 0.037 117.40 1.4e-05 11.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
6758 rs72654902 1_14 1.000 155.97 5.0e-04 11.58
6797 rs2873296 1_14 1.000 97.04 3.1e-04 -2.62
6832 rs1780323 1_15 1.000 1267.71 4.0e-03 34.33
6834 rs192212396 1_15 1.000 206.66 6.6e-04 -11.43
6851 rs4654961 1_15 1.000 171.12 5.4e-04 -15.27
6876 rs34745509 1_15 1.000 195.56 6.2e-04 6.73
6881 rs4654973 1_15 1.000 417.00 1.3e-03 -15.10
13771 rs3176461 1_32 1.000 43.61 1.4e-04 -6.81
37286 rs28383573 1_78 1.000 96.74 3.1e-04 -9.56
62906 rs12044944 1_126 1.000 34.83 1.1e-04 5.79
95036 rs6730773 2_57 1.000 48.37 1.5e-04 7.83
112390 rs4664862 2_94 1.000 104.31 3.3e-04 8.98
175786 rs7638322 3_75 1.000 52.32 1.7e-04 -9.09
175829 rs76947531 3_76 1.000 90.11 2.9e-04 9.65
196410 rs56328339 3_115 1.000 79.61 2.5e-04 9.80
305769 rs4958244 5_80 1.000 51.07 1.6e-04 7.25
364935 rs199804242 6_89 1.000 126187.85 4.0e-01 11.21
451853 rs6994124 8_54 1.000 52.19 1.7e-04 7.31
480505 rs10115804 9_12 1.000 36.77 1.2e-04 -5.46
599078 rs61909253 12_5 1.000 430.40 1.4e-03 17.10
629712 rs7315286 12_65 1.000 61.64 2.0e-04 -8.67
748219 rs11078597 17_2 1.000 117.46 3.7e-04 10.97
827695 rs209955 20_32 1.000 53.31 1.7e-04 7.71
883936 rs114024083 6_27 1.000 68.21 2.2e-04 -8.23
977067 rs67479069 17_27 1.000 55245.58 1.8e-01 -9.47
328935 rs115740542 6_20 0.999 33.69 1.1e-04 -5.14
559883 rs4243928 11_9 0.999 47.53 1.5e-04 7.03
566358 rs369062552 11_21 0.999 34.18 1.1e-04 5.30
566368 rs34830202 11_21 0.999 41.02 1.3e-04 -6.45
594668 rs4937122 11_77 0.999 32.22 1.0e-04 5.47
976374 rs79406732 17_27 0.999 54818.62 1.7e-01 -9.37
977907 rs140412994 17_27 0.999 54794.05 1.7e-01 -9.32
491179 rs11144105 9_35 0.998 31.13 9.9e-05 5.40
739144 rs12595893 16_39 0.998 44.57 1.4e-04 -5.95
804914 rs35496032 19_26 0.998 31.12 9.9e-05 -5.37
13797 rs57401684 1_32 0.997 55.86 1.8e-04 5.79
965197 rs34933034 15_35 0.997 86.96 2.8e-04 -10.60
747919 rs2663349 17_1 0.996 32.77 1.0e-04 5.64
334013 rs1187117 6_28 0.995 90.71 2.9e-04 8.77
269464 rs835146 5_12 0.994 29.33 9.3e-05 4.80
374128 rs4709741 6_106 0.993 28.53 9.0e-05 5.18
134829 rs7584554 2_137 0.992 98.47 3.1e-04 11.53
171298 rs16853539 3_67 0.992 28.11 8.9e-05 -5.08
215117 rs16852326 4_33 0.992 29.59 9.3e-05 5.29
390948 rs6974574 7_28 0.991 41.92 1.3e-04 -6.05
599104 rs145615184 12_5 0.990 50.03 1.6e-04 3.19
51469 rs1105822 1_104 0.988 27.70 8.7e-05 -5.15
337095 rs75111243 6_35 0.988 37.47 1.2e-04 -6.31
307699 rs34433004 5_84 0.987 34.54 1.1e-04 -5.78
364951 rs6923513 6_89 0.987 125824.41 3.9e-01 11.61
6724 rs3026894 1_14 0.986 85.98 2.7e-04 -1.50
637210 rs2701627 12_80 0.986 39.00 1.2e-04 -7.88
848284 rs5754136 22_12 0.986 27.06 8.5e-05 -4.97
356878 rs2354558 6_71 0.983 26.22 8.2e-05 4.88
480565 rs33917322 9_12 0.982 28.74 9.0e-05 -5.11
823528 rs4812458 20_24 0.981 36.98 1.2e-04 5.70
6765 rs77025042 1_14 0.979 132.88 4.1e-04 8.34
369648 rs3020333 6_99 0.978 76.75 2.4e-04 -10.38
312981 rs13167291 5_93 0.975 34.24 1.1e-04 -5.51
630511 rs653178 12_67 0.975 107.90 3.3e-04 -10.17
758370 rs2931274 17_26 0.974 32.22 1.0e-04 5.11
175836 rs60992117 3_76 0.970 49.02 1.5e-04 7.54
364934 rs2327654 6_89 0.970 125820.92 3.9e-01 11.60
503317 rs1326895 9_56 0.970 26.30 8.1e-05 -4.84
543837 rs78723267 10_64 0.968 51.48 1.6e-04 -7.23
629614 rs11608460 12_65 0.968 32.76 1.0e-04 -5.73
276947 rs12518871 5_25 0.967 35.07 1.1e-04 -6.28
357236 rs9285397 6_73 0.967 91.69 2.8e-04 -10.11
408187 rs17308346 7_58 0.967 26.14 8.0e-05 4.87
878342 rs116402366 5_106 0.966 67.59 2.1e-04 -9.24
476442 rs447124 9_5 0.964 28.01 8.6e-05 5.11
599118 rs2970810 12_5 0.961 720.82 2.2e-03 32.03
761940 rs201245566 17_35 0.958 53.32 1.6e-04 -7.93
585305 rs11021233 11_54 0.957 31.60 9.6e-05 -5.56
741309 rs78547885 16_42 0.957 25.50 7.8e-05 -4.73
135009 rs11563208 2_137 0.954 27.33 8.3e-05 -4.37
793255 rs351988 19_2 0.954 27.01 8.2e-05 -4.92
350130 rs182595 6_58 0.953 27.51 8.3e-05 -5.13
599112 rs11063183 12_5 0.952 713.42 2.2e-03 31.75
431423 rs11373012 8_11 0.941 40.55 1.2e-04 7.34
599097 rs11063147 12_5 0.941 250.36 7.5e-04 0.10
629689 rs34263583 12_65 0.939 25.54 7.6e-05 4.68
733044 rs17616063 16_27 0.936 25.93 7.7e-05 4.84
36019 rs2990245 1_76 0.935 60.98 1.8e-04 8.48
375749 rs6910424 6_110 0.934 24.78 7.4e-05 -4.39
625959 rs10777805 12_57 0.933 24.83 7.4e-05 -4.65
823750 rs117184383 20_25 0.932 25.82 7.6e-05 -5.09
827699 rs2585441 20_32 0.932 25.16 7.5e-05 5.04
409471 rs7809629 7_63 0.929 38.32 1.1e-04 6.08
501766 rs4743022 9_54 0.926 24.81 7.3e-05 4.42
364988 rs4504488 6_90 0.920 34.71 1.0e-04 -6.87
739160 rs12931591 16_39 0.914 24.61 7.1e-05 3.34
536809 rs768395 10_50 0.910 25.40 7.3e-05 4.70
692956 rs143811999 14_45 0.906 24.22 7.0e-05 -4.52
803992 rs886136 19_24 0.903 33.97 9.8e-05 -6.15
8133 rs6694059 1_18 0.902 27.16 7.8e-05 -5.27
453200 rs2943540 8_56 0.901 53.86 1.5e-04 7.40
229880 rs2869683 4_59 0.899 70.96 2.0e-04 -7.90
614215 rs150158762 12_33 0.894 24.30 6.9e-05 4.31
794132 rs375325 19_3 0.890 24.01 6.8e-05 -4.47
608685 rs7302975 12_21 0.886 25.16 7.1e-05 4.74
670721 rs7143633 14_1 0.885 23.99 6.7e-05 4.40
515034 rs17486892 10_9 0.878 29.30 8.2e-05 4.88
658271 rs1330025 13_40 0.877 32.72 9.1e-05 -5.62
881975 rs55919316 6_27 0.871 29.86 8.3e-05 -4.58
363757 rs3843996 6_87 0.866 251.70 6.9e-04 -6.79
391893 rs62449682 7_30 0.862 25.16 6.9e-05 -4.59
761279 rs2097730 17_33 0.856 24.28 6.6e-05 -4.40
241098 rs1993191 4_80 0.855 27.21 7.4e-05 -4.40
37279 rs12731846 1_78 0.854 95.82 2.6e-04 11.21
482284 rs776756 9_14 0.851 28.13 7.6e-05 -5.09
147654 rs13325064 3_18 0.846 47.70 1.3e-04 -7.06
383549 rs542176135 7_17 0.846 30.00 8.1e-05 5.28
68704 rs896788 2_5 0.845 47.72 1.3e-04 7.90
312932 rs11135047 5_93 0.841 33.73 9.0e-05 5.34
423729 rs10282631 7_95 0.832 25.81 6.8e-05 4.82
515048 rs72786681 10_9 0.832 81.65 2.2e-04 -9.64
732628 rs11076504 16_27 0.824 36.88 9.7e-05 5.94
820982 rs6083843 20_19 0.811 26.44 6.8e-05 -4.54
93929 rs6745529 2_54 0.809 36.48 9.4e-05 6.05
365987 rs5880424 6_92 0.804 25.62 6.5e-05 4.76
548015 rs758211 10_71 0.801 44.25 1.1e-04 -6.61
555171 rs72636980 11_1 0.801 24.31 6.2e-05 -4.07SNPs 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
364935 rs199804242 6_89 1.000 126187.85 4.0e-01 11.21
364951 rs6923513 6_89 0.987 125824.41 3.9e-01 11.61
364934 rs2327654 6_89 0.970 125820.92 3.9e-01 11.60
364938 rs113527452 6_89 0.000 125174.57 0.0e+00 11.55
364943 rs200796875 6_89 0.000 124450.01 0.0e+00 11.55
364956 rs7756915 6_89 0.000 123651.62 0.0e+00 11.59
364949 rs6570040 6_89 0.000 118660.21 0.0e+00 11.37
364936 rs6570031 6_89 0.000 118358.15 0.0e+00 11.16
364937 rs9389323 6_89 0.000 118322.80 0.0e+00 11.21
364953 rs9321531 6_89 0.000 103847.28 0.0e+00 10.63
364926 rs9321528 6_89 0.000 102477.29 0.0e+00 10.40
364954 rs9494389 6_89 0.000 97542.91 0.0e+00 10.53
364958 rs5880262 6_89 0.000 97295.68 0.0e+00 9.90
364931 rs1033755 6_89 0.000 94099.59 0.0e+00 10.14
364932 rs2208574 6_89 0.000 94088.62 0.0e+00 9.84
364940 rs9494377 6_89 0.000 92369.93 0.0e+00 9.71
364929 rs2038551 6_89 0.000 92266.01 0.0e+00 9.48
364927 rs2038550 6_89 0.000 92026.32 0.0e+00 9.50
364916 rs6570026 6_89 0.000 76433.03 0.0e+00 9.69
364912 rs6926161 6_89 0.000 75437.19 0.0e+00 9.39
364921 rs6930773 6_89 0.000 74000.10 0.0e+00 8.67
364908 rs6925959 6_89 0.000 63812.73 0.0e+00 8.69
364913 rs6917005 6_89 0.000 61658.92 0.0e+00 8.59
977067 rs67479069 17_27 1.000 55245.58 1.8e-01 -9.47
976374 rs79406732 17_27 0.999 54818.62 1.7e-01 -9.37
977907 rs140412994 17_27 0.999 54794.05 1.7e-01 -9.32
977049 rs62063305 17_27 0.000 54745.28 1.2e-10 -9.38
977062 rs17651213 17_27 0.000 54745.19 1.1e-10 -9.38
977063 rs17572361 17_27 0.000 54745.19 1.1e-10 -9.38
977047 rs17572248 17_27 0.000 54745.18 1.1e-10 -9.38
977089 rs17572627 17_27 0.000 54745.18 2.5e-10 -9.39
977050 rs17651134 17_27 0.000 54745.16 1.0e-10 -9.38
977051 rs62063306 17_27 0.000 54745.16 1.0e-10 -9.38
977052 rs62063774 17_27 0.000 54745.16 1.0e-10 -9.38
977053 rs62063775 17_27 0.000 54745.16 1.0e-10 -9.38
977055 rs62063776 17_27 0.000 54745.16 1.0e-10 -9.38
977066 rs62063777 17_27 0.000 54745.16 1.0e-10 -9.38
977068 rs2217394 17_27 0.000 54745.16 1.0e-10 -9.38
977064 rs17651243 17_27 0.000 54745.15 1.0e-10 -9.38
977070 rs17651285 17_27 0.000 54745.13 1.0e-10 -9.38
977069 rs62063778 17_27 0.000 54745.11 9.5e-11 -9.38
977079 rs62063780 17_27 0.000 54745.11 2.1e-10 -9.38
977071 rs17572467 17_27 0.000 54745.10 1.0e-10 -9.38
977073 rs17572495 17_27 0.000 54744.93 1.0e-10 -9.38
977080 rs2163129 17_27 0.000 54744.93 1.6e-10 -9.38
977086 rs17572613 17_27 0.000 54744.85 2.0e-10 -9.38
977087 rs77527347 17_27 0.000 54744.85 1.7e-10 -9.38
977074 rs62063779 17_27 0.000 54744.79 9.8e-11 -9.38
977141 rs372810927 17_27 0.000 54744.71 1.3e-10 -9.38
977093 rs17651483 17_27 0.000 54744.69 1.2e-10 -9.38SNPs 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
364935 rs199804242 6_89 1.000 126187.85 0.40000 11.21
364934 rs2327654 6_89 0.970 125820.92 0.39000 11.60
364951 rs6923513 6_89 0.987 125824.41 0.39000 11.61
977067 rs67479069 17_27 1.000 55245.58 0.18000 -9.47
976374 rs79406732 17_27 0.999 54818.62 0.17000 -9.37
977907 rs140412994 17_27 0.999 54794.05 0.17000 -9.32
977405 rs62062325 17_27 0.459 54476.54 0.08000 -9.47
977404 rs62062324 17_27 0.451 54476.37 0.07800 -9.47
977406 rs62062133 17_27 0.089 54575.84 0.01600 -9.53
6832 rs1780323 1_15 1.000 1267.71 0.00400 34.33
599112 rs11063183 12_5 0.952 713.42 0.00220 31.75
599118 rs2970810 12_5 0.961 720.82 0.00220 32.03
333790 rs16869466 6_28 0.347 1274.66 0.00140 -39.18
333791 rs73743336 6_28 0.336 1274.58 0.00140 -39.17
599078 rs61909253 12_5 1.000 430.40 0.00140 17.10
6881 rs4654973 1_15 1.000 417.00 0.00130 -15.10
333788 rs73743328 6_28 0.317 1274.52 0.00130 -39.17
599097 rs11063147 12_5 0.941 250.36 0.00075 0.10
363757 rs3843996 6_87 0.866 251.70 0.00069 -6.79
6834 rs192212396 1_15 1.000 206.66 0.00066 -11.43
6876 rs34745509 1_15 1.000 195.56 0.00062 6.73
6851 rs4654961 1_15 1.000 171.12 0.00054 -15.27
6758 rs72654902 1_14 1.000 155.97 0.00050 11.58
878355 rs10051765 5_106 0.775 202.12 0.00050 -19.91
6765 rs77025042 1_14 0.979 132.88 0.00041 8.34
363769 rs111394979 6_87 0.342 342.10 0.00037 17.75
748219 rs11078597 17_2 1.000 117.46 0.00037 10.97
112390 rs4664862 2_94 1.000 104.31 0.00033 8.98
630511 rs653178 12_67 0.975 107.90 0.00033 -10.17
6787 rs72657133 1_14 0.561 181.57 0.00032 12.62
6797 rs2873296 1_14 1.000 97.04 0.00031 -2.62
37286 rs28383573 1_78 1.000 96.74 0.00031 -9.56
134829 rs7584554 2_137 0.992 98.47 0.00031 11.53
363768 rs34976633 6_87 0.540 182.06 0.00031 14.48
175829 rs76947531 3_76 1.000 90.11 0.00029 9.65
334013 rs1187117 6_28 0.995 90.71 0.00029 8.77
357236 rs9285397 6_73 0.967 91.69 0.00028 -10.11
965197 rs34933034 15_35 0.997 86.96 0.00028 -10.60
6724 rs3026894 1_14 0.986 85.98 0.00027 -1.50
37279 rs12731846 1_78 0.854 95.82 0.00026 11.21
363770 rs3911914 6_87 0.450 181.71 0.00026 14.44
601797 rs10743976 12_11 0.617 134.87 0.00026 12.20
196410 rs56328339 3_115 1.000 79.61 0.00025 9.80
6780 rs531429510 1_14 0.415 180.59 0.00024 12.56
369648 rs3020333 6_99 0.978 76.75 0.00024 -10.38
175845 rs13085498 3_76 0.514 143.68 0.00023 -13.95
977929 rs113667149 17_27 0.001 54411.82 0.00023 -9.56
175847 rs10934583 3_76 0.486 143.86 0.00022 -13.96
515048 rs72786681 10_9 0.832 81.65 0.00022 -9.64
883936 rs114024083 6_27 1.000 68.21 0.00022 -8.23SNPs 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
333790 rs16869466 6_28 0.347 1274.66 1.4e-03 -39.18
333788 rs73743328 6_28 0.317 1274.52 1.3e-03 -39.17
333791 rs73743336 6_28 0.336 1274.58 1.4e-03 -39.17
6832 rs1780323 1_15 1.000 1267.71 4.0e-03 34.33
333772 rs73743305 6_28 0.000 944.98 6.1e-08 -33.79
599118 rs2970810 12_5 0.961 720.82 2.2e-03 32.03
599122 rs10849086 12_5 0.042 716.18 9.5e-05 32.01
599120 rs2907498 12_5 0.029 714.90 6.5e-05 31.99
599119 rs11503123 12_5 0.006 709.32 1.4e-05 31.90
599112 rs11063183 12_5 0.952 713.42 2.2e-03 31.75
599115 rs11063194 12_5 0.005 704.09 1.0e-05 31.53
6812 rs12047493 1_15 0.000 1273.11 0.0e+00 31.39
6825 rs3820292 1_15 0.000 1238.46 0.0e+00 31.24
599109 rs10849077 12_5 0.000 443.64 7.7e-13 25.58
333823 rs767893 6_28 0.000 420.69 5.4e-09 -22.48
599116 rs11063200 12_5 0.000 323.40 1.7e-10 22.46
599110 rs7965800 12_5 0.000 317.39 1.5e-10 22.33
333786 rs568901 6_28 0.006 400.30 8.1e-06 -22.14
599125 rs7295624 12_5 0.000 276.41 1.2e-11 21.86
6822 rs7554122 1_15 0.000 479.30 0.0e+00 21.70
6817 rs10632199 1_15 0.000 478.35 0.0e+00 21.69
6828 rs2047653 1_15 0.000 446.16 0.0e+00 21.31
6826 rs12083322 1_15 0.000 709.64 5.0e-19 -21.09
878355 rs10051765 5_106 0.775 202.12 5.0e-04 -19.91
878351 rs56235845 5_106 0.225 199.63 1.4e-04 -19.87
6835 rs1780328 1_15 0.000 544.68 0.0e+00 19.76
6843 rs141227415 1_15 0.000 523.56 0.0e+00 19.38
333826 rs73747331 6_28 0.000 300.45 2.3e-09 -19.15
878357 rs11748297 5_106 0.000 158.26 8.9e-09 -18.80
6880 rs148785605 1_15 0.000 686.30 4.3e-11 18.77
878350 rs4074995 5_106 0.000 155.53 8.6e-09 -18.74
878356 rs11748165 5_106 0.000 151.80 8.3e-09 -18.55
6830 rs75824579 1_15 0.000 459.41 0.0e+00 -18.25
878352 rs11746443 5_106 0.000 146.06 7.9e-09 -18.21
878346 rs11741640 5_106 0.000 141.75 7.6e-09 -18.15
878347 rs12654812 5_106 0.000 157.10 8.6e-09 -18.12
878332 rs4075958 5_106 0.000 140.09 7.5e-09 -18.08
599084 rs16931344 12_5 0.000 301.22 6.2e-13 -18.04
599081 rs3812822 12_5 0.000 300.92 6.1e-13 -18.01
599088 rs11063130 12_5 0.000 298.73 5.9e-13 -17.97
363769 rs111394979 6_87 0.342 342.10 3.7e-04 17.75
363772 rs79557158 6_87 0.134 339.75 1.4e-04 17.72
363776 rs113653360 6_87 0.128 339.64 1.4e-04 17.72
363771 rs117404455 6_87 0.135 339.77 1.5e-04 17.71
363773 rs4897538 6_87 0.131 339.69 1.4e-04 17.71
363774 rs4897539 6_87 0.131 339.69 1.4e-04 17.71
175863 rs111928357 3_76 0.277 142.49 1.3e-04 -17.70
175869 rs112779893 3_76 0.257 142.20 1.2e-04 -17.69
175867 rs73186059 3_76 0.189 141.11 8.5e-05 -17.66
878345 rs11135015 5_106 0.000 137.18 7.2e-09 -17.50Gene 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] 23if (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 insulin-like growth factor receptor signaling pathway (GO:0048009)
2 positive regulation of cyclin-dependent protein serine/threonine kinase activity (GO:0045737)
3 positive regulation of cyclin-dependent protein kinase activity (GO:1904031)
4 positive regulation of ERBB signaling pathway (GO:1901186)
5 positive regulation of G1/S transition of mitotic cell cycle (GO:1900087)
6 positive regulation of epidermal growth factor receptor signaling pathway (GO:0045742)
7 positive regulation of cell cycle G1/S phase transition (GO:1902808)
Overlap Adjusted.P.value Genes
1 2/11 0.02482067 GIGYF1;AKT1
2 2/17 0.02840221 CCND2;AKT1
3 2/20 0.02840221 CCND2;AKT1
4 2/25 0.02902766 AKT1;CBL
5 2/26 0.02902766 CCND2;AKT1
6 2/34 0.03772124 AKT1;CBL
7 2/35 0.03772124 CCND2;AKT1
[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-219B17.3 gene(s) from the input list not found in DisGeNET CURATEDPALM3 gene(s) from the input list not found in DisGeNET CURATEDSUPT3H gene(s) from the input list not found in DisGeNET CURATEDGRHL1 gene(s) from the input list not found in DisGeNET CURATEDPRIM1 gene(s) from the input list not found in DisGeNET CURATEDC11orf63 gene(s) from the input list not found in DisGeNET CURATEDSYTL1 gene(s) from the input list not found in DisGeNET CURATEDMIATNB gene(s) from the input list not found in DisGeNET CURATEDTMEM176A gene(s) from the input list not found in DisGeNET CURATEDRP11-1109F11.3 gene(s) from the input list not found in DisGeNET CURATEDELP5 gene(s) from the input list not found in DisGeNET CURATEDGIGYF1 gene(s) from the input list not found in DisGeNET CURATED Description
25 Fibrosis
90 Microcornea
134 ovarian neoplasm
135 Malignant neoplasm of ovary
155 Cirrhosis
168 NOONAN SYNDROME-LIKE DISORDER WITH OR WITHOUT JUVENILE MYELOMONOCYTIC LEUKEMIA
171 COWDEN SYNDROME 6
175 MEGALENCEPHALY-POLYMICROGYRIA-POLYDACTYLY-HYDROCEPHALUS SYNDROME 3
179 Fetal hydrops (in some patients)
183 INTELLECTUAL DEVELOPMENTAL DISORDER, AUTOSOMAL RECESSIVE 68
FDR Ratio BgRatio
25 0.02543853 2/11 50/9703
90 0.02543853 1/11 1/9703
134 0.02543853 3/11 134/9703
135 0.02543853 3/11 137/9703
155 0.02543853 2/11 50/9703
168 0.02543853 1/11 1/9703
171 0.02543853 1/11 1/9703
175 0.02543853 1/11 1/9703
179 0.02543853 1/11 1/9703
183 0.02543853 1/11 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