Glucose (quantile) - Whole_Blood
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-30500_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30600_irnt_Liver.Rmd
Modified: analysis/ukb-d-30600_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30610_irnt_Liver.Rmd
Modified: analysis/ukb-d-30610_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30620_irnt_Liver.Rmd
Modified: analysis/ukb-d-30620_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30630_irnt_Liver.Rmd
Modified: analysis/ukb-d-30630_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30640_irnt_Liver.Rmd
Modified: analysis/ukb-d-30640_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30650_irnt_Liver.Rmd
Modified: analysis/ukb-d-30650_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30660_irnt_Liver.Rmd
Modified: analysis/ukb-d-30660_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30670_irnt_Liver.Rmd
Modified: analysis/ukb-d-30670_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30680_irnt_Liver.Rmd
Modified: analysis/ukb-d-30690_irnt_Liver.Rmd
Modified: analysis/ukb-d-30690_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30700_irnt_Liver.Rmd
Modified: analysis/ukb-d-30700_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30710_irnt_Liver.Rmd
Modified: analysis/ukb-d-30710_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30720_irnt_Liver.Rmd
Modified: analysis/ukb-d-30720_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30730_irnt_Liver.Rmd
Modified: analysis/ukb-d-30740_irnt_Liver.Rmd
Modified: analysis/ukb-d-30740_irnt_Whole_Blood.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
Modified: analysis/ukb-d-30820_irnt_Liver.Rmd
Modified: analysis/ukb-d-30830_irnt_Liver.Rmd
Modified: analysis/ukb-d-30840_irnt_Liver.Rmd
Modified: analysis/ukb-d-30850_irnt_Liver.Rmd
Modified: analysis/ukb-d-30860_irnt_Liver.Rmd
Modified: analysis/ukb-d-30870_irnt_Liver.Rmd
Modified: analysis/ukb-d-30880_irnt_Liver.Rmd
Modified: analysis/ukb-d-30890_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-30740_irnt_Whole_Blood.Rmd) and HTML (docs/ukb-d-30740_irnt_Whole_Blood.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 | 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 Glucose (quantile) using Whole_Blood 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-30740_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 Whole_Blood 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] 11095#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
1129 747 624 400 479 621 560 383 404 430 682 652 192 362 331
16 17 18 19 20 21 22
551 725 159 911 313 130 310 #proportion of imputed weights without missing variants
mean(qclist_all$nmiss==0)[1] 0.762776Load 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.010986634 0.000147187 #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
8.491741 11.082572 #report sample size
print(sample_size)[1] 314916#report group size
group_size <- c(nrow(ctwas_gene_res), n_snps)
print(group_size)[1] 11095 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.003286957 0.045050671 #compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))[1] 0.01977635 0.34170116Genes 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
2550 MADD 11_29 0.991 42.25 1.3e-04 9.18
6290 ZFP36L2 2_27 0.921 69.43 2.0e-04 10.18
7872 IGF2 11_2 0.919 101.36 3.0e-04 -9.28
10110 C15orf52 15_14 0.872 25.19 7.0e-05 -4.80
1693 PTK6 20_37 0.798 21.59 5.5e-05 -4.58
3262 SGIP1 1_42 0.790 25.11 6.3e-05 -4.95
11139 NPEPL1 20_34 0.788 20.12 5.0e-05 3.38
9643 AIFM3 22_4 0.788 22.53 5.6e-05 -4.70
1667 PABPC1L 20_28 0.785 43.46 1.1e-04 -6.88
2787 NNT 5_28 0.775 20.93 5.2e-05 -4.13
8719 CHD2 15_43 0.775 21.68 5.3e-05 -4.23
5922 GIGYF1 7_62 0.711 43.06 9.7e-05 -6.95
4397 H3F3B 17_42 0.628 21.96 4.4e-05 4.20
6089 FADS1 11_34 0.576 50.73 9.3e-05 -7.62
7996 FAM234A 16_1 0.546 49.05 8.5e-05 8.01
2989 USP34 2_40 0.526 23.05 3.9e-05 3.75
1366 CWF19L1 10_64 0.514 19.13 3.1e-05 -3.90
4963 LRRC1 6_40 0.495 23.63 3.7e-05 4.61
10164 RNFT1 17_35 0.489 20.20 3.1e-05 -4.24
6952 LRWD1 7_63 0.479 24.60 3.7e-05 3.83Genes 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
12599 HCP5B 6_26 0 2027.17 0 4.41
10848 TRIM10 6_26 0 1340.43 0 -3.71
10855 HLA-G 6_26 0 1290.24 0 -5.57
10853 HCG9 6_26 0 820.97 0 0.99
10968 HLA-A 6_26 0 705.76 0 1.28
10844 HLA-E 6_26 0 570.91 0 1.29
11418 TRIM26 6_26 0 563.57 0 -0.45
11120 LINC00243 6_26 0 552.63 0 4.13
5868 PPP1R18 6_26 0 393.37 0 -0.54
10841 MRPS18B 6_26 0 301.95 0 -0.03
11652 C4A 6_26 0 150.35 0 -6.04
11218 C4B 6_26 0 147.72 0 5.98
7712 C2 6_26 0 146.07 0 5.95
10847 TRIM15 6_26 0 143.95 0 -1.81
11047 CLIC1 6_26 0 142.44 0 -5.87
10808 NEU1 6_26 0 139.47 0 -5.81
10825 APOM 6_26 0 133.32 0 -5.67
4971 IER3 6_26 0 125.90 0 0.96
10789 PBX2 6_26 0 116.23 0 -5.28
11374 CYP21A2 6_26 0 103.60 0 4.98Genes with highest PVE
#genes with 20 highest pve
head(ctwas_gene_res[order(-ctwas_gene_res$PVE),report_cols],20) genename region_tag susie_pip mu2 PVE z
7872 IGF2 11_2 0.919 101.36 3.0e-04 -9.28
6290 ZFP36L2 2_27 0.921 69.43 2.0e-04 10.18
2550 MADD 11_29 0.991 42.25 1.3e-04 9.18
1667 PABPC1L 20_28 0.785 43.46 1.1e-04 -6.88
5922 GIGYF1 7_62 0.711 43.06 9.7e-05 -6.95
6089 FADS1 11_34 0.576 50.73 9.3e-05 -7.62
7996 FAM234A 16_1 0.546 49.05 8.5e-05 8.01
10110 C15orf52 15_14 0.872 25.19 7.0e-05 -4.80
3262 SGIP1 1_42 0.790 25.11 6.3e-05 -4.95
9643 AIFM3 22_4 0.788 22.53 5.6e-05 -4.70
1693 PTK6 20_37 0.798 21.59 5.5e-05 -4.58
8719 CHD2 15_43 0.775 21.68 5.3e-05 -4.23
2787 NNT 5_28 0.775 20.93 5.2e-05 -4.13
11139 NPEPL1 20_34 0.788 20.12 5.0e-05 3.38
4397 H3F3B 17_42 0.628 21.96 4.4e-05 4.20
6761 UBE2Z 17_28 0.424 32.25 4.3e-05 -5.67
9050 FBXO46 19_32 0.460 28.27 4.1e-05 -5.98
2861 OGG1 3_8 0.421 29.90 4.0e-05 3.43
2989 USP34 2_40 0.526 23.05 3.9e-05 3.75
3804 OPRL1 20_38 0.372 32.31 3.8e-05 3.65Genes 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
2236 YKT6 7_32 0.004 41.56 4.8e-07 -11.82
6290 ZFP36L2 2_27 0.921 69.43 2.0e-04 10.18
2953 NRBP1 2_16 0.043 97.49 1.3e-05 -10.09
7872 IGF2 11_2 0.919 101.36 3.0e-04 -9.28
2550 MADD 11_29 0.991 42.25 1.3e-04 9.18
2956 SNX17 2_16 0.039 83.76 1.0e-05 -9.16
3490 SEC22A 3_76 0.038 60.53 7.3e-06 8.18
7996 FAM234A 16_1 0.546 49.05 8.5e-05 8.01
7255 EIF5A2 3_104 0.016 39.41 2.0e-06 7.88
9564 CAMK1D 10_10 0.018 50.84 2.9e-06 7.87
6089 FADS1 11_34 0.576 50.73 9.3e-05 -7.62
3631 KBTBD4 11_29 0.024 28.80 2.2e-06 -7.05
5922 GIGYF1 7_62 0.711 43.06 9.7e-05 -6.95
1667 PABPC1L 20_28 0.785 43.46 1.1e-04 -6.88
970 UBE2D4 7_32 0.000 37.96 1.6e-11 -6.85
4610 ACP2 11_29 0.055 34.33 6.0e-06 -6.79
2497 FNBP4 11_29 0.019 24.80 1.5e-06 6.50
7304 ZNF513 2_16 0.017 37.72 2.0e-06 -6.27
7821 YWHAB 20_28 0.037 36.75 4.3e-06 6.26
8876 ARHGAP1 11_28 0.166 43.25 2.3e-05 -6.22Comparing 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.006219018#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
2236 YKT6 7_32 0.004 41.56 4.8e-07 -11.82
6290 ZFP36L2 2_27 0.921 69.43 2.0e-04 10.18
2953 NRBP1 2_16 0.043 97.49 1.3e-05 -10.09
7872 IGF2 11_2 0.919 101.36 3.0e-04 -9.28
2550 MADD 11_29 0.991 42.25 1.3e-04 9.18
2956 SNX17 2_16 0.039 83.76 1.0e-05 -9.16
3490 SEC22A 3_76 0.038 60.53 7.3e-06 8.18
7996 FAM234A 16_1 0.546 49.05 8.5e-05 8.01
7255 EIF5A2 3_104 0.016 39.41 2.0e-06 7.88
9564 CAMK1D 10_10 0.018 50.84 2.9e-06 7.87
6089 FADS1 11_34 0.576 50.73 9.3e-05 -7.62
3631 KBTBD4 11_29 0.024 28.80 2.2e-06 -7.05
5922 GIGYF1 7_62 0.711 43.06 9.7e-05 -6.95
1667 PABPC1L 20_28 0.785 43.46 1.1e-04 -6.88
970 UBE2D4 7_32 0.000 37.96 1.6e-11 -6.85
4610 ACP2 11_29 0.055 34.33 6.0e-06 -6.79
2497 FNBP4 11_29 0.019 24.80 1.5e-06 6.50
7304 ZNF513 2_16 0.017 37.72 2.0e-06 -6.27
7821 YWHAB 20_28 0.037 36.75 4.3e-06 6.26
8876 ARHGAP1 11_28 0.166 43.25 2.3e-05 -6.22Locus 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: 7_32"
genename region_tag susie_pip mu2 PVE z
2226 COA1 7_32 0.000 9.20 4.3e-12 0.06
2227 BLVRA 7_32 0.000 13.04 1.2e-11 1.08
565 MRPS24 7_32 0.000 11.89 1.4e-11 0.07
970 UBE2D4 7_32 0.000 37.96 1.6e-11 -6.85
2228 URGCP 7_32 0.000 29.68 1.3e-11 5.98
11314 AC004951.6 7_32 0.000 8.88 4.2e-12 -0.37
4841 DBNL 7_32 0.000 24.69 2.4e-11 -1.07
3572 POLM 7_32 0.000 10.02 8.4e-12 -0.57
2232 AEBP1 7_32 0.000 75.46 1.9e-10 -1.82
2233 POLD2 7_32 0.000 27.84 5.8e-10 0.92
2236 YKT6 7_32 0.004 41.56 4.8e-07 -11.82
523 CAMK2B 7_32 0.000 69.97 3.2e-08 5.10
4838 DDX56 7_32 0.000 28.29 1.9e-10 2.87
6708 TMED4 7_32 0.000 28.29 1.9e-10 2.87
2151 OGDH 7_32 0.000 20.23 6.5e-11 -1.77
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 2_27"
genename region_tag susie_pip mu2 PVE z
11353 AC093609.1 2_27 0.011 4.69 1.6e-07 0.51
12581 LINC01126 2_27 0.084 30.16 8.0e-06 -4.46
6290 ZFP36L2 2_27 0.921 69.43 2.0e-04 10.18
3047 THADA 2_27 0.046 30.66 4.5e-06 5.53
6292 PLEKHH2 2_27 0.120 30.42 1.2e-05 4.07
5065 DYNC2LI1 2_27 0.022 12.27 8.5e-07 -1.86
5077 LRPPRC 2_27 0.012 6.50 2.5e-07 -1.13
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 2_16"
genename region_tag susie_pip mu2 PVE z
11045 SLC35F6 2_16 0.022 14.70 1.0e-06 4.06
3366 TMEM214 2_16 0.037 14.82 1.8e-06 -3.75
5074 EMILIN1 2_16 0.020 10.91 7.0e-07 -4.40
5061 KHK 2_16 0.033 14.43 1.5e-06 3.75
5059 CGREF1 2_16 0.020 11.55 7.4e-07 2.81
5070 PREB 2_16 0.023 9.40 6.9e-07 0.81
5076 ATRAID 2_16 0.018 22.16 1.3e-06 4.81
1090 CAD 2_16 0.022 8.16 5.8e-07 1.60
5071 SLC5A6 2_16 0.017 10.47 5.6e-07 -2.97
7303 UCN 2_16 0.017 9.04 5.0e-07 2.74
2952 GTF3C2 2_16 0.017 9.09 5.0e-07 -2.75
2956 SNX17 2_16 0.039 83.76 1.0e-05 -9.16
7304 ZNF513 2_16 0.017 37.72 2.0e-06 -6.27
2953 NRBP1 2_16 0.043 97.49 1.3e-05 -10.09
5057 IFT172 2_16 0.022 14.62 1.0e-06 3.49
1087 GCKR 2_16 0.023 15.17 1.1e-06 -3.55
10613 GPN1 2_16 0.019 14.69 8.7e-07 -2.91
9018 CCDC121 2_16 0.036 9.35 1.1e-06 0.10
6660 BRE 2_16 0.022 11.99 8.6e-07 -2.77
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 11_2"
genename region_tag susie_pip mu2 PVE z
969 TOLLIP 11_2 0.002 6.93 3.5e-08 0.78
9483 MOB2 11_2 0.003 10.51 8.6e-08 -1.26
9709 DUSP8 11_2 0.002 9.91 7.3e-08 -1.15
11638 IFITM10 11_2 0.001 4.62 1.8e-08 -0.32
3233 CTSD 11_2 0.001 5.91 2.7e-08 -0.71
4206 TNNI2 11_2 0.002 7.38 3.8e-08 0.99
4204 LSP1 11_2 0.001 5.14 2.2e-08 0.14
4205 TNNT3 11_2 0.002 8.50 5.5e-08 -0.79
12620 PRR33 11_2 0.002 6.71 3.4e-08 0.63
11337 LINC01150 11_2 0.002 7.83 4.7e-08 0.56
11078 MRPL23 11_2 0.005 16.89 2.5e-07 -1.96
7872 IGF2 11_2 0.919 101.36 3.0e-04 -9.28
9638 ASCL2 11_2 0.001 4.63 1.8e-08 -0.64
2557 C11orf21 11_2 0.001 5.18 2.2e-08 -0.16
588 TSPAN32 11_2 0.001 5.86 2.6e-08 -0.84
2555 CD81 11_2 0.002 7.84 4.4e-08 -0.94
9683 TSSC4 11_2 0.002 9.15 5.6e-08 1.26
4133 CDKN1C 11_2 0.001 5.36 2.1e-08 -0.87
2554 SLC22A18 11_2 0.078 32.29 8.0e-06 2.28
11813 SLC22A18AS 11_2 0.041 27.65 3.6e-06 2.14
9430 PHLDA2 11_2 0.001 4.77 1.9e-08 -0.11
10926 NAP1L4 11_2 0.001 4.79 1.9e-08 0.43
2553 CARS 11_2 0.004 14.78 1.8e-07 1.48
272 OSBPL5 11_2 0.001 5.63 2.4e-08 0.87
9688 MRGPRE 11_2 0.001 5.11 2.1e-08 -0.45
76 ZNF195 11_2 0.008 21.18 5.6e-07 -2.12
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 11_29"
genename region_tag susie_pip mu2 PVE z
6066 ARFGAP2 11_29 0.024 12.36 9.5e-07 3.25
300 NR1H3 11_29 0.051 18.49 3.0e-06 -3.60
4610 ACP2 11_29 0.055 34.33 6.0e-06 -6.79
2550 MADD 11_29 0.991 42.25 1.3e-04 9.18
4609 MYBPC3 11_29 0.022 5.90 4.1e-07 0.19
7654 PSMC3 11_29 0.039 14.63 1.8e-06 -3.30
7653 SLC39A13 11_29 0.028 24.42 2.2e-06 -6.20
7655 RAPSN 11_29 0.041 15.78 2.0e-06 -3.30
2551 PTPMT1 11_29 0.018 6.68 3.8e-07 -1.73
3631 KBTBD4 11_29 0.024 28.80 2.2e-06 -7.05
8552 C1QTNF4 11_29 0.040 13.68 1.7e-06 -1.71
7656 AGBL2 11_29 0.021 5.78 3.9e-07 0.26
2497 FNBP4 11_29 0.019 24.80 1.5e-06 6.50
324 NUP160 11_29 0.027 9.29 7.9e-07 -2.16
6064 PTPRJ 11_29 0.021 7.24 4.9e-07 2.04
| 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
54771 rs79687284 1_108 1.000 111.12 3.5e-04 12.08
75845 rs780093 2_16 1.000 190.20 6.0e-04 14.95
81784 rs2121564 2_28 1.000 60.29 1.9e-04 8.01
114586 rs12692596 2_96 1.000 46.08 1.5e-04 7.24
116878 rs11396827 2_102 1.000 209.10 6.6e-04 20.35
116880 rs71397673 2_102 1.000 477.91 1.5e-03 34.50
116881 rs537183 2_102 1.000 875.78 2.8e-03 52.66
116888 rs853789 2_102 1.000 903.71 2.9e-03 53.04
167265 rs148685409 3_57 1.000 907.17 2.9e-03 2.99
177379 rs72964564 3_76 1.000 237.79 7.6e-04 -16.88
197429 rs4234603 3_115 1.000 39.59 1.3e-04 5.05
324284 rs76623841 6_7 1.000 57.90 1.8e-04 -6.80
383850 rs10225316 7_15 1.000 45.31 1.4e-04 7.51
394087 rs138917529 7_32 1.000 78.46 2.5e-04 -10.45
433869 rs7012814 8_12 1.000 62.91 2.0e-04 -9.18
463448 rs146191002 8_70 1.000 641.34 2.0e-03 -0.15
529445 rs61856594 10_33 1.000 38.41 1.2e-04 -6.23
549258 rs12244851 10_70 1.000 294.45 9.4e-04 16.46
559313 rs3750952 11_7 1.000 52.32 1.7e-04 -7.40
645076 rs576674 13_10 1.000 105.10 3.3e-04 -10.82
694074 rs35889227 14_45 1.000 116.43 3.7e-04 -11.34
865929 rs1611236 6_26 1.000 4824.25 1.5e-02 3.25
288819 rs12189028 5_45 0.999 34.80 1.1e-04 -5.08
478405 rs10758593 9_4 0.999 71.27 2.3e-04 8.64
54780 rs3754140 1_108 0.998 58.70 1.9e-04 -10.01
116860 rs11676084 2_102 0.998 127.50 4.0e-04 -23.20
236488 rs11728350 4_69 0.997 42.57 1.3e-04 6.68
512031 rs115478735 9_70 0.997 82.62 2.6e-04 9.52
560735 rs34718245 11_9 0.997 32.29 1.0e-04 -5.37
757231 rs28489441 17_15 0.997 32.63 1.0e-04 -5.83
324303 rs55792466 6_7 0.994 98.16 3.1e-04 -9.67
485647 rs572168822 9_16 0.993 41.18 1.3e-04 -6.61
659724 rs1327315 13_40 0.992 34.84 1.1e-04 -7.02
636194 rs4765221 12_76 0.989 33.66 1.1e-04 5.80
394067 rs17769733 7_32 0.986 134.24 4.2e-04 -7.76
548710 rs11195508 10_70 0.986 60.91 1.9e-04 -10.76
836100 rs6026545 20_34 0.985 35.46 1.1e-04 5.72
411067 rs4729755 7_63 0.984 26.54 8.3e-05 4.96
476363 rs4977218 8_94 0.979 29.83 9.3e-05 -5.34
705074 rs12912777 15_13 0.977 25.47 7.9e-05 3.77
485642 rs1333045 9_16 0.974 40.85 1.3e-04 6.62
627692 rs6538805 12_58 0.973 32.06 9.9e-05 -6.74
900411 rs231362 11_2 0.971 48.07 1.5e-04 6.83
755 rs60330317 1_2 0.969 36.84 1.1e-04 -6.18
192270 rs10653660 3_104 0.969 351.02 1.1e-03 -23.54
643386 rs60353775 13_7 0.968 82.69 2.5e-04 9.78
175412 rs9875598 3_73 0.963 27.30 8.3e-05 -5.10
467109 rs4433184 8_78 0.962 53.14 1.6e-04 4.78
671654 rs80081165 13_62 0.956 25.21 7.7e-05 4.82
893019 rs11257655 10_10 0.950 72.24 2.2e-04 9.12
34661 rs893230 1_72 0.948 39.75 1.2e-04 -7.56
394079 rs10259649 7_32 0.948 397.09 1.2e-03 28.65
513321 rs28624681 9_73 0.942 52.20 1.6e-04 7.55
701551 rs35767992 15_4 0.940 24.63 7.4e-05 4.72
454712 rs10957704 8_54 0.931 24.45 7.2e-05 4.68
758740 rs543720569 17_18 0.931 45.59 1.3e-04 -7.03
334606 rs62396405 6_30 0.930 25.29 7.5e-05 -4.78
571670 rs117396352 11_28 0.930 27.25 8.1e-05 4.93
574874 rs7941126 11_36 0.912 31.50 9.1e-05 -5.63
357128 rs118126621 6_73 0.906 24.98 7.2e-05 4.67
572690 rs182512331 11_31 0.901 27.67 7.9e-05 -5.05
571867 rs7111517 11_28 0.899 39.18 1.1e-04 -6.64
170183 rs62276527 3_63 0.897 33.73 9.6e-05 5.85
363608 rs112388031 6_87 0.896 24.52 7.0e-05 -4.65
561667 rs117720468 11_11 0.892 44.58 1.3e-04 6.82
288880 rs6887019 5_45 0.884 26.65 7.5e-05 5.23
155694 rs3172494 3_34 0.881 26.79 7.5e-05 -5.09
796854 rs41404946 18_44 0.876 24.22 6.7e-05 4.56
639226 rs1882297 12_82 0.865 38.60 1.1e-04 6.47
572272 rs139913257 11_30 0.858 31.56 8.6e-05 -5.63
258624 rs62332172 4_113 0.852 24.65 6.7e-05 4.54
660075 rs79317015 13_40 0.849 24.17 6.5e-05 4.40
800071 rs10410896 19_4 0.825 28.22 7.4e-05 5.04
396423 rs11763778 7_36 0.823 44.87 1.2e-04 -7.61
116936 rs112308555 2_103 0.820 24.85 6.5e-05 4.32
196969 rs73185688 3_114 0.815 25.17 6.5e-05 4.69
328895 rs2206734 6_15 0.811 64.37 1.7e-04 8.26
712481 rs11637069 15_29 0.809 28.58 7.3e-05 5.00
576904 rs11603349 11_41 0.808 45.18 1.2e-04 -6.78SNPs 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
865929 rs1611236 6_26 1.000 4824.25 1.5e-02 3.25
865911 rs2508055 6_26 0.121 4795.98 1.8e-03 3.28
865914 rs111734624 6_26 0.120 4795.98 1.8e-03 3.28
865951 rs1611252 6_26 0.113 4795.93 1.7e-03 3.28
865965 rs1611260 6_26 0.111 4795.91 1.7e-03 3.28
865970 rs1611265 6_26 0.110 4795.90 1.7e-03 3.28
865947 rs1611248 6_26 0.106 4795.88 1.6e-03 3.28
865972 rs1611267 6_26 0.091 4795.69 1.4e-03 3.27
865908 rs1737020 6_26 0.094 4795.68 1.4e-03 3.27
865909 rs1737019 6_26 0.094 4795.68 1.4e-03 3.27
865870 rs2844838 6_26 0.094 4795.66 1.4e-03 3.27
865973 rs2394171 6_26 0.085 4795.64 1.3e-03 3.27
865916 rs1611228 6_26 0.086 4795.62 1.3e-03 3.27
865975 rs2893981 6_26 0.085 4795.62 1.3e-03 3.27
865874 rs1633032 6_26 0.101 4795.44 1.5e-03 3.28
865861 rs1633033 6_26 0.062 4795.36 9.5e-04 3.26
865906 rs1633018 6_26 0.094 4795.13 1.4e-03 3.27
865903 rs1633020 6_26 0.090 4795.12 1.4e-03 3.27
865815 rs1610726 6_26 0.128 4794.89 1.9e-03 3.29
865868 rs2844840 6_26 0.101 4794.24 1.5e-03 3.28
865927 rs1611234 6_26 0.040 4794.15 6.1e-04 3.24
866101 rs3129185 6_26 0.077 4793.96 1.2e-03 3.27
865945 rs1611246 6_26 0.154 4793.87 2.3e-03 3.30
866225 rs1632980 6_26 0.124 4793.55 1.9e-03 3.29
866108 rs1736999 6_26 0.055 4793.48 8.3e-04 3.26
866116 rs1633001 6_26 0.066 4793.27 1.0e-03 3.27
865893 rs1614309 6_26 0.061 4791.72 9.2e-04 3.27
865892 rs1633030 6_26 0.142 4789.05 2.2e-03 3.31
865983 rs9258382 6_26 0.026 4782.71 4.0e-04 3.26
865980 rs9258379 6_26 0.091 4776.89 1.4e-03 3.34
865939 rs1611241 6_26 0.095 4772.39 1.4e-03 3.37
865895 rs1633028 6_26 0.001 4762.19 1.8e-05 3.22
865941 rs1611244 6_26 0.000 4744.42 4.2e-06 3.23
865904 rs2735042 6_26 0.000 4736.29 2.6e-06 3.20
865971 rs1611266 6_26 0.052 4710.22 7.8e-04 3.56
865948 rs1611249 6_26 0.002 4688.22 2.4e-05 3.50
865920 rs1611230 6_26 0.002 4677.86 2.5e-05 3.53
865962 rs145043018 6_26 0.002 4676.89 2.3e-05 3.53
865969 rs147376303 6_26 0.002 4676.87 2.3e-05 3.53
865978 rs9258376 6_26 0.002 4676.81 2.3e-05 3.53
865984 rs1633016 6_26 0.002 4676.76 2.3e-05 3.53
865884 rs1618670 6_26 0.001 4676.25 2.1e-05 3.53
865856 rs1633035 6_26 0.001 4676.24 1.4e-05 3.52
866021 rs1633014 6_26 0.001 4676.10 1.6e-05 3.52
865905 rs1633019 6_26 0.001 4675.96 2.0e-05 3.53
866091 rs1633010 6_26 0.001 4674.93 1.9e-05 3.54
866134 rs1736991 6_26 0.002 4674.59 3.1e-05 3.56
866182 rs1610713 6_26 0.002 4674.48 2.6e-05 3.55
866157 rs909722 6_26 0.002 4674.43 2.3e-05 3.55
866175 rs1610653 6_26 0.002 4674.06 2.6e-05 3.55SNPs 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
865929 rs1611236 6_26 1.000 4824.25 0.01500 3.25
116888 rs853789 2_102 1.000 903.71 0.00290 53.04
167265 rs148685409 3_57 1.000 907.17 0.00290 2.99
116881 rs537183 2_102 1.000 875.78 0.00280 52.66
865945 rs1611246 6_26 0.154 4793.87 0.00230 3.30
865892 rs1633030 6_26 0.142 4789.05 0.00220 3.31
463448 rs146191002 8_70 1.000 641.34 0.00200 -0.15
865815 rs1610726 6_26 0.128 4794.89 0.00190 3.29
866225 rs1632980 6_26 0.124 4793.55 0.00190 3.29
865911 rs2508055 6_26 0.121 4795.98 0.00180 3.28
865914 rs111734624 6_26 0.120 4795.98 0.00180 3.28
865951 rs1611252 6_26 0.113 4795.93 0.00170 3.28
865965 rs1611260 6_26 0.111 4795.91 0.00170 3.28
865970 rs1611265 6_26 0.110 4795.90 0.00170 3.28
463439 rs72681356 8_70 0.796 635.51 0.00160 4.78
865947 rs1611248 6_26 0.106 4795.88 0.00160 3.28
116880 rs71397673 2_102 1.000 477.91 0.00150 34.50
167267 rs1436648 3_57 0.507 929.81 0.00150 -3.09
167268 rs7619398 3_57 0.522 930.03 0.00150 -3.08
865868 rs2844840 6_26 0.101 4794.24 0.00150 3.28
865874 rs1633032 6_26 0.101 4795.44 0.00150 3.28
865870 rs2844838 6_26 0.094 4795.66 0.00140 3.27
865903 rs1633020 6_26 0.090 4795.12 0.00140 3.27
865906 rs1633018 6_26 0.094 4795.13 0.00140 3.27
865908 rs1737020 6_26 0.094 4795.68 0.00140 3.27
865909 rs1737019 6_26 0.094 4795.68 0.00140 3.27
865939 rs1611241 6_26 0.095 4772.39 0.00140 3.37
865972 rs1611267 6_26 0.091 4795.69 0.00140 3.27
865980 rs9258379 6_26 0.091 4776.89 0.00140 3.34
394092 rs917793 7_32 0.583 708.36 0.00130 33.28
865916 rs1611228 6_26 0.086 4795.62 0.00130 3.27
865973 rs2394171 6_26 0.085 4795.64 0.00130 3.27
865975 rs2893981 6_26 0.085 4795.62 0.00130 3.27
394079 rs10259649 7_32 0.948 397.09 0.00120 28.65
866101 rs3129185 6_26 0.077 4793.96 0.00120 3.27
167269 rs2256473 3_57 0.366 929.35 0.00110 -3.06
192270 rs10653660 3_104 0.969 351.02 0.00110 -23.54
866116 rs1633001 6_26 0.066 4793.27 0.00100 3.27
463441 rs72681364 8_70 0.483 633.93 0.00097 4.74
865861 rs1633033 6_26 0.062 4795.36 0.00095 3.26
549258 rs12244851 10_70 1.000 294.45 0.00094 16.46
167266 rs2575789 3_57 0.311 929.24 0.00092 -3.06
865893 rs1614309 6_26 0.061 4791.72 0.00092 3.27
866108 rs1736999 6_26 0.055 4793.48 0.00083 3.26
865971 rs1611266 6_26 0.052 4710.22 0.00078 3.56
177379 rs72964564 3_76 1.000 237.79 0.00076 -16.88
116878 rs11396827 2_102 1.000 209.10 0.00066 20.35
865927 rs1611234 6_26 0.040 4794.15 0.00061 3.24
75845 rs780093 2_16 1.000 190.20 0.00060 14.95
394096 rs2908282 7_32 0.257 706.28 0.00058 33.26SNPs 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
116888 rs853789 2_102 1.000 903.71 2.9e-03 53.04
116881 rs537183 2_102 1.000 875.78 2.8e-03 52.66
116882 rs518598 2_102 0.000 829.26 3.8e-08 52.13
116884 rs485094 2_102 0.000 725.02 1.3e-10 51.04
116886 rs2544360 2_102 0.000 573.16 1.6e-10 39.60
116887 rs853791 2_102 0.000 543.15 7.9e-11 39.21
116890 rs853785 2_102 0.000 413.09 1.8e-11 37.45
116889 rs860510 2_102 0.000 400.58 1.9e-11 37.03
116883 rs579275 2_102 0.000 416.54 1.7e-11 36.60
116880 rs71397673 2_102 1.000 477.91 1.5e-03 34.50
394092 rs917793 7_32 0.583 708.36 1.3e-03 33.28
394096 rs2908282 7_32 0.257 706.28 5.8e-04 33.26
394086 rs4607517 7_32 0.159 706.24 3.6e-04 33.23
394098 rs732360 7_32 0.000 695.17 8.8e-07 32.89
394079 rs10259649 7_32 0.948 397.09 1.2e-03 28.65
394077 rs2908294 7_32 0.052 387.36 6.4e-05 28.27
116874 rs62176784 2_102 0.001 136.22 3.7e-07 -25.09
116876 rs549410 2_102 0.000 188.39 3.1e-08 -23.64
192270 rs10653660 3_104 0.969 351.02 1.1e-03 -23.54
192280 rs8192675 3_104 0.035 343.83 3.8e-05 -23.31
116860 rs11676084 2_102 0.998 127.50 4.0e-04 -23.20
116878 rs11396827 2_102 1.000 209.10 6.6e-04 20.35
116861 rs2140046 2_102 0.000 66.41 3.5e-09 -18.78
192278 rs11920090 3_104 0.233 161.93 1.2e-04 -18.33
192264 rs12492910 3_104 0.183 160.97 9.3e-05 -18.32
192277 rs11923694 3_104 0.161 160.50 8.2e-05 -18.31
192267 rs12496506 3_104 0.172 160.59 8.8e-05 -18.30
192281 rs11928798 3_104 0.121 158.98 6.1e-05 -18.26
192282 rs6785803 3_104 0.124 158.98 6.3e-05 -18.25
192257 rs56351320 3_104 0.002 214.32 1.2e-06 -17.54
192242 rs6792607 3_104 0.005 202.53 3.5e-06 -17.19
549266 rs12260037 10_70 0.000 251.89 7.5e-10 17.19
116879 rs13430620 2_102 0.000 90.17 8.3e-11 -16.98
383910 rs1974619 7_15 0.608 243.14 4.7e-04 16.89
177379 rs72964564 3_76 1.000 237.79 7.6e-04 -16.88
383908 rs10228796 7_15 0.220 240.65 1.7e-04 16.82
383909 rs2191349 7_15 0.172 240.16 1.3e-04 16.80
116885 rs114932341 2_102 0.000 185.76 4.1e-10 16.61
383911 rs188745922 7_15 0.009 234.23 6.8e-06 16.59
549258 rs12244851 10_70 1.000 294.45 9.4e-04 16.46
192234 rs11919048 3_104 0.015 136.31 6.4e-06 -16.13
177377 rs34642857 3_76 0.005 216.60 3.6e-06 -16.11
192249 rs79560566 3_104 0.001 127.99 4.3e-07 -15.48
383906 rs6461153 7_15 0.002 207.12 1.4e-06 15.42
383905 rs10266209 7_15 0.002 206.75 1.4e-06 15.41
383904 rs10249299 7_15 0.004 207.33 2.4e-06 -15.31
75845 rs780093 2_16 1.000 190.20 6.0e-04 14.95
192272 rs143791579 3_104 0.001 97.73 3.4e-07 -14.91
383899 rs4721398 7_15 0.002 190.52 1.2e-06 -14.88
192265 rs73167792 3_104 0.001 94.78 3.6e-07 -14.81Gene 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] 4if (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 cell differentiation (GO:0045596)
2 cellular response to granulocyte macrophage colony-stimulating factor stimulus (GO:0097011)
3 response to granulocyte macrophage colony-stimulating factor (GO:0097012)
4 regulation of Rab protein signal transduction (GO:0032483)
5 definitive hemopoiesis (GO:0060216)
6 regulation of glycogen (starch) synthase activity (GO:2000465)
7 negative regulation of cell cycle phase transition (GO:1901988)
8 negative regulation of muscle cell differentiation (GO:0051148)
9 genetic imprinting (GO:0071514)
10 regulation of gene expression by genetic imprinting (GO:0006349)
11 negative regulation of stem cell differentiation (GO:2000737)
12 positive regulation of vascular endothelial cell proliferation (GO:1905564)
13 positive regulation of insulin receptor signaling pathway (GO:0046628)
14 regulation of histone modification (GO:0031056)
15 positive regulation of glycogen biosynthetic process (GO:0045725)
16 regulation of chromatin organization (GO:1902275)
17 positive regulation of cellular response to insulin stimulus (GO:1900078)
18 T cell differentiation in thymus (GO:0033077)
19 positive regulation of glycogen metabolic process (GO:0070875)
20 positive regulation of nuclear-transcribed mRNA catabolic process, deadenylation-dependent decay (GO:1900153)
21 regulation of nuclear-transcribed mRNA catabolic process, deadenylation-dependent decay (GO:1900151)
22 3'-UTR-mediated mRNA destabilization (GO:0061158)
23 cellular response to corticosteroid stimulus (GO:0071384)
24 cellular response to glucocorticoid stimulus (GO:0071385)
25 regulation of B cell differentiation (GO:0045577)
26 regulation of vascular endothelial cell proliferation (GO:1905562)
27 positive regulation of activated T cell proliferation (GO:0042104)
28 regulation of lymphocyte differentiation (GO:0045619)
29 ERK1 and ERK2 cascade (GO:0070371)
30 in utero embryonic development (GO:0001701)
31 regulation of glycogen biosynthetic process (GO:0005979)
32 positive regulation of nuclear division (GO:0051785)
33 response to glucocorticoid (GO:0051384)
34 regulation of B cell activation (GO:0050864)
35 cellular response to epidermal growth factor stimulus (GO:0071364)
36 regulation of activated T cell proliferation (GO:0046006)
37 regulation of muscle cell differentiation (GO:0051147)
38 negative regulation of fat cell differentiation (GO:0045599)
39 positive regulation of mitotic nuclear division (GO:0045840)
40 mRNA destabilization (GO:0061157)
41 regulation of apoptotic signaling pathway (GO:2001233)
42 positive regulation of catalytic activity (GO:0043085)
43 T cell differentiation (GO:0030217)
44 regulation of extrinsic apoptotic signaling pathway via death domain receptors (GO:1902041)
45 mRNA catabolic process (GO:0006402)
46 positive regulation of mRNA catabolic process (GO:0061014)
47 regulation of insulin receptor signaling pathway (GO:0046626)
48 negative regulation of mitotic cell cycle (GO:0045930)
49 RNA catabolic process (GO:0006401)
50 regulation of protein modification process (GO:0031399)
51 regulation of mitotic nuclear division (GO:0007088)
52 chordate embryonic development (GO:0043009)
53 regulation of tumor necrosis factor-mediated signaling pathway (GO:0010803)
54 regulation of extrinsic apoptotic signaling pathway (GO:2001236)
55 positive regulation of T cell proliferation (GO:0042102)
56 insulin receptor signaling pathway (GO:0008286)
57 regulation of cytokine-mediated signaling pathway (GO:0001959)
58 positive regulation of endothelial cell proliferation (GO:0001938)
59 regulation of fat cell differentiation (GO:0045598)
60 regulation of gene expression, epigenetic (GO:0040029)
61 regulation of Ras protein signal transduction (GO:0046578)
62 regulation of stem cell differentiation (GO:2000736)
63 cellular response to fibroblast growth factor stimulus (GO:0044344)
64 negative regulation of mitotic cell cycle phase transition (GO:1901991)
65 regulation of peptidyl-tyrosine phosphorylation (GO:0050730)
66 hemopoiesis (GO:0030097)
67 mRNA metabolic process (GO:0016071)
68 positive regulation of cell cycle process (GO:0090068)
69 response to tumor necrosis factor (GO:0034612)
70 cellular response to transforming growth factor beta stimulus (GO:0071560)
71 regulation of mRNA catabolic process (GO:0061013)
72 platelet degranulation (GO:0002576)
73 cellular response to insulin stimulus (GO:0032869)
74 positive regulation of macromolecule biosynthetic process (GO:0010557)
75 positive regulation of peptidyl-tyrosine phosphorylation (GO:0050731)
76 regulation of mRNA stability (GO:0043488)
77 positive regulation of transferase activity (GO:0051347)
78 regulation of cell differentiation (GO:0045595)
79 cellular response to growth factor stimulus (GO:0071363)
80 skeletal system development (GO:0001501)
81 positive regulation of protein kinase B signaling (GO:0051897)
82 regulation of MAPK cascade (GO:0043408)
83 positive regulation of cellular biosynthetic process (GO:0031328)
84 regulated exocytosis (GO:0045055)
85 regulation of mitotic cell cycle phase transition (GO:1901990)
Overlap Adjusted.P.value Genes
1 2/191 0.01576426 IGF2;ZFP36L2
2 1/6 0.01576426 ZFP36L2
3 1/6 0.01576426 ZFP36L2
4 1/7 0.01576426 MADD
5 1/7 0.01576426 ZFP36L2
6 1/7 0.01576426 IGF2
7 1/7 0.01576426 ZFP36L2
8 1/8 0.01576426 IGF2
9 1/10 0.01576426 IGF2
10 1/11 0.01576426 IGF2
11 1/13 0.01576426 ZFP36L2
12 1/13 0.01576426 IGF2
13 1/13 0.01576426 IGF2
14 1/13 0.01576426 IGF2
15 1/14 0.01576426 IGF2
16 1/14 0.01576426 IGF2
17 1/14 0.01576426 IGF2
18 1/14 0.01576426 ZFP36L2
19 1/15 0.01576426 IGF2
20 1/15 0.01576426 ZFP36L2
21 1/15 0.01576426 ZFP36L2
22 1/16 0.01576426 ZFP36L2
23 1/16 0.01576426 ZFP36L2
24 1/18 0.01576426 ZFP36L2
25 1/18 0.01576426 ZFP36L2
26 1/18 0.01576426 IGF2
27 1/20 0.01626228 IGF2
28 1/20 0.01626228 ZFP36L2
29 1/23 0.01805269 ZFP36L2
30 1/25 0.01820286 IGF2
31 1/26 0.01820286 IGF2
32 1/27 0.01820286 IGF2
33 1/27 0.01820286 ZFP36L2
34 1/28 0.01820286 ZFP36L2
35 1/28 0.01820286 ZFP36L2
36 1/34 0.02099072 IGF2
37 1/35 0.02099072 IGF2
38 1/36 0.02099072 ZFP36L2
39 1/36 0.02099072 IGF2
40 1/38 0.02107289 ZFP36L2
41 1/38 0.02107289 MADD
42 1/40 0.02165060 IGF2
43 1/41 0.02167416 ZFP36L2
44 1/42 0.02169656 MADD
45 1/44 0.02173822 ZFP36L2
46 1/44 0.02173822 ZFP36L2
47 1/45 0.02175761 IGF2
48 1/48 0.02271781 ZFP36L2
49 1/49 0.02271781 ZFP36L2
50 1/51 0.02316869 IGF2
51 1/57 0.02484423 IGF2
52 1/58 0.02484423 IGF2
53 1/58 0.02484423 MADD
54 1/64 0.02689455 MADD
55 1/66 0.02722665 IGF2
56 1/73 0.02943806 IGF2
57 1/74 0.02943806 MADD
58 1/77 0.03009660 IGF2
59 1/80 0.03073229 ZFP36L2
60 1/82 0.03097094 IGF2
61 1/86 0.03193965 MADD
62 1/91 0.03205092 ZFP36L2
63 1/92 0.03205092 ZFP36L2
64 1/92 0.03205092 ZFP36L2
65 1/92 0.03205092 IGF2
66 1/94 0.03210089 ZFP36L2
67 1/95 0.03210089 ZFP36L2
68 1/101 0.03361129 IGF2
69 1/110 0.03605148 ZFP36L2
70 1/114 0.03681764 ZFP36L2
71 1/122 0.03882306 ZFP36L2
72 1/125 0.03921643 IGF2
73 1/129 0.03936571 IGF2
74 1/129 0.03936571 IGF2
75 1/134 0.04033116 IGF2
76 1/146 0.04332566 ZFP36L2
77 1/148 0.04334228 IGF2
78 1/156 0.04450223 IGF2
79 1/158 0.04450223 ZFP36L2
80 1/158 0.04450223 IGF2
81 1/161 0.04477728 IGF2
82 1/166 0.04558773 IGF2
83 1/180 0.04820481 IGF2
84 1/180 0.04820481 IGF2
85 1/188 0.04972504 ZFP36L2
[1] "GO_Cellular_Component_2021"
Term Overlap Adjusted.P.value Genes
1 platelet alpha granule lumen (GO:0031093) 1/67 0.04470013 IGF2
2 platelet alpha granule (GO:0031091) 1/90 0.04470013 IGF2
[1] "GO_Molecular_Function_2021"
Term Overlap
1 protein kinase activator activity (GO:0030295) 2/63
2 insulin-like growth factor receptor binding (GO:0005159) 1/14
3 death receptor binding (GO:0005123) 1/15
4 mRNA 3'-UTR AU-rich region binding (GO:0035925) 1/22
5 tumor necrosis factor receptor superfamily binding (GO:0032813) 1/28
6 kinase activator activity (GO:0019209) 1/31
7 protein serine/threonine kinase activator activity (GO:0043539) 1/37
8 mRNA 3'-UTR binding (GO:0003730) 1/85
9 growth factor activity (GO:0008083) 1/87
10 protein kinase regulator activity (GO:0019887) 1/98
11 guanyl-nucleotide exchange factor activity (GO:0005085) 1/149
Adjusted.P.value Genes
1 0.0008169314 MADD;IGF2
2 0.0139850795 IGF2
3 0.0139850795 MADD
4 0.0144339668 ZFP36L2
5 0.0144339668 MADD
6 0.0144339668 MADD
7 0.0147599130 IGF2
8 0.0268923548 ZFP36L2
9 0.0268923548 IGF2
10 0.0272407938 MADD
11 0.0375080843 MADDMADD gene(s) from the input list not found in DisGeNET CURATEDC15orf52 gene(s) from the input list not found in DisGeNET CURATED Description FDR Ratio
13 Colorectal Neoplasms 0.005679508 2/2
18 Polyhydramnios 0.005679508 1/2
26 Placenta Disorders 0.005679508 1/2
39 Congenital hemihypertrophy 0.005679508 1/2
42 Radiolabeled somatostatin analog study 0.005679508 1/2
53 MENTAL RETARDATION, X-LINKED, SNYDER-ROBINSON TYPE 0.005679508 1/2
55 Fetus Small for Gestational Age 0.005679508 1/2
56 HEMIHYPERPLASIA, ISOLATED 0.005679508 1/2
61 GROWTH RESTRICTION, SEVERE, WITH DISTINCTIVE FACIES 0.005679508 1/2
45 prenatal alcohol exposure 0.006389117 1/2
BgRatio
13 277/9703
18 1/9703
26 4/9703
39 4/9703
42 3/9703
53 4/9703
55 2/9703
56 4/9703
61 1/9703
45 5/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