Glucose (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
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_Liver.Rmd) and HTML (docs/ukb-d-30740_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 Glucose (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-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 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.0091025951 0.0001641765 #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
12.62894 10.12580 #report sample size
print(sample_size)[1] 314916#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.003979272 0.045912578 #compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))[1] 0.01661887 0.30199457Genes 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
3212 CCND2 12_4 0.993 120.36 3.8e-04 -11.04
5991 FADS1 11_34 0.989 62.63 2.0e-04 8.16
8531 TNKS 8_12 0.982 55.36 1.7e-04 -9.51
77 ABCC8 11_12 0.952 22.91 6.9e-05 4.16
5632 CAND2 3_9 0.938 33.28 9.9e-05 -5.66
189 QPCTL 19_32 0.923 37.75 1.1e-04 -6.51
3716 PPDPF 20_37 0.798 23.11 5.9e-05 -4.60
3175 SGIP1 1_42 0.775 26.28 6.5e-05 -4.94
5822 GIGYF1 7_62 0.768 42.08 1.0e-04 -6.95
7394 TP53INP1 8_66 0.764 35.82 8.7e-05 -6.53
9140 TH 11_2 0.757 36.51 8.8e-05 1.67
8716 ARHGAP1 11_28 0.753 40.28 9.6e-05 -6.27
5465 PRCC 1_77 0.739 25.99 6.1e-05 4.88
6702 ACE 17_37 0.720 22.66 5.2e-05 -4.54
3666 KCNK17 6_30 0.695 23.62 5.2e-05 3.94
2288 DNAJC12 10_44 0.679 27.55 5.9e-05 5.15
8345 RGS19 20_38 0.633 24.92 5.0e-05 -3.86
5174 WARS 14_52 0.628 25.43 5.1e-05 -4.88
11727 LINC01863 5_104 0.609 22.62 4.4e-05 3.83
8434 DNAJB7 22_17 0.597 25.49 4.8e-05 -4.33Genes 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
7118 SLC2A2 3_104 0.029 140.93 1.3e-05 16.86
3212 CCND2 12_4 0.993 120.36 3.8e-04 -11.04
500 CAMK2B 7_32 0.000 111.74 7.0e-08 7.60
2887 NRBP1 2_16 0.041 110.50 1.4e-05 -10.47
7119 RPL22L1 3_104 0.036 107.45 1.2e-05 10.33
2183 AEBP1 7_32 0.000 105.98 5.4e-10 -4.21
6177 SPC25 2_102 0.000 93.08 3.9e-10 -16.98
12661 LINC01126 2_27 0.074 85.95 2.0e-05 -9.57
2185 MYL7 7_32 0.000 84.18 4.2e-11 12.48
2891 SNX17 2_16 0.094 77.31 2.3e-05 8.04
11916 LINC00261 20_16 0.141 70.80 3.2e-05 8.50
5991 FADS1 11_34 0.989 62.63 2.0e-04 8.16
3747 FOXA2 20_16 0.016 59.56 3.1e-06 -7.98
2186 GCK 7_32 0.000 58.74 1.4e-08 -2.52
8531 TNKS 8_12 0.982 55.36 1.7e-04 -9.51
6046 ADRA2A 10_70 0.000 53.96 5.5e-08 -3.47
11738 RP11-115J16.2 8_12 0.019 53.82 3.2e-06 -8.49
9887 NCR3LG1 11_12 0.592 53.29 1.0e-04 -7.16
721 WIPI1 17_39 0.424 48.77 6.6e-05 -6.92
4507 FADS2 11_34 0.010 45.88 1.4e-06 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
3212 CCND2 12_4 0.993 120.36 3.8e-04 -11.04
5991 FADS1 11_34 0.989 62.63 2.0e-04 8.16
8531 TNKS 8_12 0.982 55.36 1.7e-04 -9.51
189 QPCTL 19_32 0.923 37.75 1.1e-04 -6.51
5822 GIGYF1 7_62 0.768 42.08 1.0e-04 -6.95
9887 NCR3LG1 11_12 0.592 53.29 1.0e-04 -7.16
5632 CAND2 3_9 0.938 33.28 9.9e-05 -5.66
8716 ARHGAP1 11_28 0.753 40.28 9.6e-05 -6.27
9140 TH 11_2 0.757 36.51 8.8e-05 1.67
7394 TP53INP1 8_66 0.764 35.82 8.7e-05 -6.53
77 ABCC8 11_12 0.952 22.91 6.9e-05 4.16
721 WIPI1 17_39 0.424 48.77 6.6e-05 -6.92
3175 SGIP1 1_42 0.775 26.28 6.5e-05 -4.94
5465 PRCC 1_77 0.739 25.99 6.1e-05 4.88
2288 DNAJC12 10_44 0.679 27.55 5.9e-05 5.15
7862 FAM234A 16_1 0.557 33.33 5.9e-05 7.59
3716 PPDPF 20_37 0.798 23.11 5.9e-05 -4.60
3684 GOT2 16_31 0.588 28.01 5.2e-05 3.93
6702 ACE 17_37 0.720 22.66 5.2e-05 -4.54
3666 KCNK17 6_30 0.695 23.62 5.2e-05 3.94Genes 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
6177 SPC25 2_102 0.000 93.08 3.9e-10 -16.98
7118 SLC2A2 3_104 0.029 140.93 1.3e-05 16.86
2185 MYL7 7_32 0.000 84.18 4.2e-11 12.48
3212 CCND2 12_4 0.993 120.36 3.8e-04 -11.04
2887 NRBP1 2_16 0.041 110.50 1.4e-05 -10.47
7119 RPL22L1 3_104 0.036 107.45 1.2e-05 10.33
12661 LINC01126 2_27 0.074 85.95 2.0e-05 -9.57
8531 TNKS 8_12 0.982 55.36 1.7e-04 -9.51
11916 LINC00261 20_16 0.141 70.80 3.2e-05 8.50
11738 RP11-115J16.2 8_12 0.019 53.82 3.2e-06 -8.49
5991 FADS1 11_34 0.989 62.63 2.0e-04 8.16
2891 SNX17 2_16 0.094 77.31 2.3e-05 8.04
3747 FOXA2 20_16 0.016 59.56 3.1e-06 -7.98
500 CAMK2B 7_32 0.000 111.74 7.0e-08 7.60
7862 FAM234A 16_1 0.557 33.33 5.9e-05 7.59
116 LUC7L 16_1 0.326 31.02 3.2e-05 7.44
4486 ACP2 11_29 0.023 39.57 2.9e-06 -7.41
623 SPI1 11_29 0.011 41.61 1.5e-06 -7.31
9887 NCR3LG1 11_12 0.592 53.29 1.0e-04 -7.16
3551 KBTBD4 11_29 0.022 40.33 2.8e-06 -7.05Comparing 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.006513164#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
6177 SPC25 2_102 0.000 93.08 3.9e-10 -16.98
7118 SLC2A2 3_104 0.029 140.93 1.3e-05 16.86
2185 MYL7 7_32 0.000 84.18 4.2e-11 12.48
3212 CCND2 12_4 0.993 120.36 3.8e-04 -11.04
2887 NRBP1 2_16 0.041 110.50 1.4e-05 -10.47
7119 RPL22L1 3_104 0.036 107.45 1.2e-05 10.33
12661 LINC01126 2_27 0.074 85.95 2.0e-05 -9.57
8531 TNKS 8_12 0.982 55.36 1.7e-04 -9.51
11916 LINC00261 20_16 0.141 70.80 3.2e-05 8.50
11738 RP11-115J16.2 8_12 0.019 53.82 3.2e-06 -8.49
5991 FADS1 11_34 0.989 62.63 2.0e-04 8.16
2891 SNX17 2_16 0.094 77.31 2.3e-05 8.04
3747 FOXA2 20_16 0.016 59.56 3.1e-06 -7.98
500 CAMK2B 7_32 0.000 111.74 7.0e-08 7.60
7862 FAM234A 16_1 0.557 33.33 5.9e-05 7.59
116 LUC7L 16_1 0.326 31.02 3.2e-05 7.44
4486 ACP2 11_29 0.023 39.57 2.9e-06 -7.41
623 SPI1 11_29 0.011 41.61 1.5e-06 -7.31
9887 NCR3LG1 11_12 0.592 53.29 1.0e-04 -7.16
3551 KBTBD4 11_29 0.022 40.33 2.8e-06 -7.05Locus 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: 2_102"
genename region_tag susie_pip mu2 PVE z
6177 SPC25 2_102 0 93.08 3.9e-10 -16.98
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 3_104"
genename region_tag susie_pip mu2 PVE z
9514 ACTRT3 3_104 0.006 5.96 1.1e-07 -1.27
8352 LRRC34 3_104 0.009 8.78 2.4e-07 -0.85
2806 LRRC31 3_104 0.006 6.14 1.2e-07 0.98
143 SEC62 3_104 0.009 8.89 2.6e-07 0.70
8591 GPR160 3_104 0.005 4.98 8.6e-08 -0.47
8590 PHC3 3_104 0.006 6.05 1.2e-07 0.03
7112 PRKCI 3_104 0.006 7.47 1.3e-07 -1.95
4734 SKIL 3_104 0.025 18.61 1.5e-06 1.90
11479 RP11-469J4.3 3_104 0.006 5.64 1.0e-07 -0.80
7119 RPL22L1 3_104 0.036 107.45 1.2e-05 10.33
7117 EIF5A2 3_104 0.011 34.70 1.2e-06 6.78
7118 SLC2A2 3_104 0.029 140.93 1.3e-05 16.86
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 7_32"
genename region_tag susie_pip mu2 PVE z
7330 STK17A 7_32 0 10.85 7.1e-12 -2.06
2177 COA1 7_32 0 11.72 6.9e-12 -1.05
2178 BLVRA 7_32 0 5.67 1.3e-12 -1.26
541 MRPS24 7_32 0 12.45 7.5e-12 0.07
2179 URGCP 7_32 0 36.25 8.0e-12 6.47
927 UBE2D4 7_32 0 16.55 3.7e-11 -1.47
11147 AC004951.6 7_32 0 9.35 2.2e-12 -0.37
4706 DBNL 7_32 0 19.16 7.7e-12 2.51
3488 POLM 7_32 0 14.48 7.5e-12 1.68
2183 AEBP1 7_32 0 105.98 5.4e-10 -4.21
2184 POLD2 7_32 0 30.12 4.9e-10 -1.10
2185 MYL7 7_32 0 84.18 4.2e-11 12.48
2186 GCK 7_32 0 58.74 1.4e-08 -2.52
500 CAMK2B 7_32 0 111.74 7.0e-08 7.60
233 NPC1L1 7_32 0 14.94 1.2e-11 0.43
4704 DDX56 7_32 0 6.66 1.6e-12 -0.02
6619 TMED4 7_32 0 17.11 1.7e-11 1.45
2101 OGDH 7_32 0 26.37 9.2e-11 -2.16
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 12_4"
genename region_tag susie_pip mu2 PVE z
4041 CRACR2A 12_4 0.009 11.14 3.0e-07 1.25
2530 PARP11 12_4 0.004 4.96 6.7e-08 0.35
11823 RP11-320N7.2 12_4 0.005 5.83 8.7e-08 0.65
3212 CCND2 12_4 0.993 120.36 3.8e-04 -11.04
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 2_16"
genename region_tag susie_pip mu2 PVE z
2881 CENPA 2_16 0.024 19.44 1.5e-06 -4.54
11149 OST4 2_16 0.015 12.93 6.2e-07 -2.90
4939 EMILIN1 2_16 0.016 8.29 4.3e-07 2.83
4927 KHK 2_16 0.025 15.33 1.2e-06 -3.75
4935 PREB 2_16 0.012 14.11 5.4e-07 -4.44
4941 ATRAID 2_16 0.012 24.97 9.4e-07 5.05
4936 SLC5A6 2_16 0.012 25.59 9.6e-07 -5.16
1060 CAD 2_16 0.014 16.66 7.6e-07 -3.47
2885 SLC30A3 2_16 0.012 22.90 8.8e-07 -4.27
7169 UCN 2_16 0.012 9.74 3.8e-07 -2.74
2891 SNX17 2_16 0.094 77.31 2.3e-05 8.04
7170 ZNF513 2_16 0.012 14.93 5.6e-07 -3.39
2887 NRBP1 2_16 0.041 110.50 1.4e-05 -10.47
4925 IFT172 2_16 0.016 15.98 8.2e-07 -3.49
1058 GCKR 2_16 0.022 43.99 3.1e-06 6.95
10987 C2orf16 2_16 0.022 43.99 3.1e-06 6.95
10407 GPN1 2_16 0.012 22.20 8.3e-07 -4.03
8847 CCDC121 2_16 0.022 9.31 6.6e-07 1.10
6575 BRE 2_16 0.023 18.57 1.3e-06 3.46
8284 RBKS 2_16 0.011 39.56 1.4e-06 -6.85
| 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
54259 rs79687284 1_108 1.000 109.90 3.5e-04 12.08
75333 rs780093 2_16 1.000 185.69 5.9e-04 14.95
81659 rs2121564 2_28 1.000 59.62 1.9e-04 8.01
114461 rs12692596 2_96 1.000 45.75 1.5e-04 7.24
116753 rs11396827 2_102 1.000 207.52 6.6e-04 20.35
116755 rs71397673 2_102 1.000 475.72 1.5e-03 34.50
116756 rs537183 2_102 1.000 874.00 2.8e-03 52.66
116763 rs853789 2_102 1.000 902.00 2.9e-03 53.04
166475 rs148685409 3_57 1.000 822.42 2.6e-03 2.99
176589 rs72964564 3_76 1.000 233.33 7.4e-04 -16.88
196639 rs4234603 3_115 1.000 38.93 1.2e-04 5.05
323494 rs76623841 6_7 1.000 56.94 1.8e-04 -6.80
385865 rs10225316 7_15 1.000 44.78 1.4e-04 7.51
396102 rs138917529 7_32 1.000 77.85 2.5e-04 -10.45
464690 rs146191002 8_70 1.000 559.76 1.8e-03 -0.15
531528 rs61856594 10_33 1.000 38.02 1.2e-04 -6.23
551341 rs12244851 10_70 1.000 290.24 9.2e-04 16.46
561396 rs3750952 11_6 1.000 52.32 1.7e-04 -7.40
573597 rs2596407 11_29 1.000 74.93 2.4e-04 9.98
647341 rs576674 13_10 1.000 103.82 3.3e-04 -10.82
696339 rs35889227 14_45 1.000 114.94 3.6e-04 -11.34
878465 rs10305492 6_30 1.000 94.37 3.0e-04 -10.36
479647 rs10758593 9_4 0.999 70.84 2.2e-04 8.64
54268 rs3754140 1_108 0.998 58.39 1.9e-04 -10.01
116735 rs11676084 2_102 0.998 127.22 4.0e-04 -23.20
288029 rs12189028 5_45 0.997 32.10 1.0e-04 -5.08
562818 rs34718245 11_9 0.997 31.60 1.0e-04 -5.37
759053 rs28489441 17_15 0.997 32.11 1.0e-04 -5.83
513273 rs115478735 9_70 0.996 81.46 2.6e-04 9.52
323513 rs55792466 6_7 0.994 96.62 3.0e-04 -9.67
486889 rs572168822 9_16 0.993 40.75 1.3e-04 -6.61
661989 rs1327315 13_40 0.991 34.59 1.1e-04 -7.02
638459 rs4765221 12_76 0.989 33.35 1.0e-04 5.80
396082 rs17769733 7_32 0.986 131.06 4.1e-04 -7.76
550793 rs11195508 10_70 0.986 60.57 1.9e-04 -10.76
413082 rs4729755 7_63 0.982 26.17 8.2e-05 4.96
707339 rs12912777 15_13 0.976 24.96 7.7e-05 3.77
486884 rs1333045 9_16 0.973 40.43 1.2e-04 6.62
629957 rs6538805 12_58 0.973 31.83 9.8e-05 -6.74
755 rs60330317 1_2 0.969 36.46 1.1e-04 -6.18
893567 rs231362 11_2 0.969 47.55 1.5e-04 6.83
191480 rs10653660 3_104 0.967 346.22 1.1e-03 -23.54
645651 rs60353775 13_7 0.967 81.75 2.5e-04 9.78
174622 rs9875598 3_73 0.963 27.01 8.3e-05 -5.10
673919 rs80081165 13_62 0.961 25.07 7.7e-05 4.82
468351 rs4433184 8_78 0.959 51.97 1.6e-04 4.78
396094 rs10259649 7_32 0.949 395.43 1.2e-03 28.65
235698 rs11728350 4_69 0.945 42.30 1.3e-04 6.68
169393 rs62276527 3_63 0.940 32.85 9.8e-05 5.85
514563 rs28624681 9_73 0.940 51.48 1.5e-04 7.55
801045 rs10410896 19_4 0.940 26.39 7.9e-05 5.04
703816 rs35767992 15_4 0.937 24.48 7.3e-05 4.72
573300 rs117396352 11_28 0.930 27.20 8.0e-05 4.93
760562 rs543720569 17_18 0.922 43.58 1.3e-04 -7.03
154904 rs3172494 3_35 0.919 26.75 7.8e-05 -5.09
577385 rs7941126 11_36 0.913 30.71 8.9e-05 -5.63
34661 rs893230 1_72 0.910 39.85 1.2e-04 -7.56
359143 rs118126621 6_73 0.905 24.74 7.1e-05 4.67
456501 rs10957704 8_54 0.904 24.91 7.1e-05 4.68
365623 rs112388031 6_87 0.902 24.11 6.9e-05 -4.65
575201 rs182512331 11_31 0.900 27.43 7.8e-05 -5.05
563750 rs117720468 11_11 0.890 44.19 1.2e-04 6.82
797828 rs41404946 18_44 0.880 23.82 6.7e-05 4.56
543881 rs11202472 10_56 0.856 30.78 8.4e-05 -4.88
574783 rs139913257 11_30 0.856 31.24 8.5e-05 -5.63
116811 rs112308555 2_103 0.846 24.06 6.5e-05 4.32
662340 rs79317015 13_40 0.840 24.23 6.5e-05 4.40
257834 rs62332172 4_113 0.833 24.60 6.5e-05 4.54
288090 rs6887019 5_45 0.830 27.96 7.4e-05 5.23
398438 rs11763778 7_36 0.817 44.49 1.2e-04 -7.61
641491 rs1882297 12_82 0.816 37.94 9.8e-05 6.47
196179 rs73185688 3_114 0.814 24.89 6.4e-05 4.69
607769 rs12582270 12_17 0.812 25.41 6.6e-05 -4.72
819032 rs6114190 20_1 0.806 25.10 6.4e-05 4.36
579415 rs11603349 11_41 0.804 44.44 1.1e-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
116763 rs853789 2_102 1.000 902.00 2.9e-03 53.04
116756 rs537183 2_102 1.000 874.00 2.8e-03 52.66
166478 rs7619398 3_57 0.515 844.12 1.4e-03 -3.08
166477 rs1436648 3_57 0.503 843.92 1.3e-03 -3.09
166479 rs2256473 3_57 0.367 843.48 9.8e-04 -3.06
166476 rs2575789 3_57 0.316 843.35 8.5e-04 -3.06
116757 rs518598 2_102 0.000 827.73 3.9e-08 52.13
166485 rs7620974 3_57 0.032 822.66 8.5e-05 2.95
166475 rs148685409 3_57 1.000 822.42 2.6e-03 2.99
166486 rs2575752 3_57 0.014 813.47 3.7e-05 -3.11
116759 rs485094 2_102 0.000 724.00 1.3e-10 51.04
396107 rs917793 7_32 0.581 701.86 1.3e-03 33.28
396111 rs2908282 7_32 0.258 699.79 5.7e-04 33.26
396101 rs4607517 7_32 0.160 699.76 3.6e-04 33.23
396113 rs732360 7_32 0.000 688.87 9.2e-07 32.89
166488 rs2575762 3_57 0.000 572.28 0.0e+00 -3.39
116761 rs2544360 2_102 0.000 569.09 1.6e-10 39.60
166500 rs9284813 3_57 0.000 560.31 0.0e+00 3.52
464690 rs146191002 8_70 1.000 559.76 1.8e-03 -0.15
464681 rs72681356 8_70 0.777 555.48 1.4e-03 4.78
464683 rs72681364 8_70 0.486 554.02 8.5e-04 4.74
116762 rs853791 2_102 0.000 539.25 7.4e-11 39.21
464691 rs17238125 8_70 0.036 535.10 6.1e-05 4.56
464707 rs72681393 8_70 0.036 534.99 6.1e-05 4.56
464693 rs17817023 8_70 0.032 534.88 5.5e-05 4.55
464698 rs72681385 8_70 0.033 534.88 5.6e-05 4.55
464700 rs11783512 8_70 0.033 534.88 5.5e-05 4.55
464702 rs72681389 8_70 0.033 534.88 5.5e-05 4.55
464709 rs55933208 8_70 0.033 534.87 5.6e-05 4.55
464694 rs17817071 8_70 0.034 534.86 5.8e-05 4.56
464701 rs11784228 8_70 0.031 534.76 5.3e-05 4.55
464713 rs56317008 8_70 0.042 534.67 7.2e-05 4.57
464715 rs72683107 8_70 0.036 534.34 6.1e-05 4.55
464704 rs56394923 8_70 0.032 534.18 5.4e-05 4.55
464714 rs55927269 8_70 0.027 533.41 4.7e-05 4.53
464717 rs2003 8_70 0.011 530.33 1.9e-05 4.44
464677 rs13275964 8_70 0.001 490.12 2.2e-06 4.06
464664 rs13252571 8_70 0.001 490.08 2.2e-06 4.05
166505 rs6775215 3_57 0.000 482.05 0.0e+00 3.37
116755 rs71397673 2_102 1.000 475.72 1.5e-03 34.50
166503 rs12631713 3_57 0.000 452.44 0.0e+00 3.29
166497 rs6791960 3_57 0.000 438.30 0.0e+00 3.17
166490 rs1036446 3_57 0.000 437.77 0.0e+00 -3.12
116758 rs579275 2_102 0.000 412.80 1.3e-11 36.60
116765 rs853785 2_102 0.000 409.63 1.3e-11 37.45
464651 rs142932284 8_70 0.001 407.81 1.7e-06 3.62
464646 rs28794838 8_70 0.001 407.34 1.7e-06 3.62
464628 rs34322085 8_70 0.001 406.63 1.7e-06 3.63
464627 rs13255888 8_70 0.001 406.39 1.7e-06 3.63
464620 rs13254156 8_70 0.001 405.71 1.7e-06 3.66SNPs 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
116763 rs853789 2_102 1.000 902.00 0.00290 53.04
116756 rs537183 2_102 1.000 874.00 0.00280 52.66
166475 rs148685409 3_57 1.000 822.42 0.00260 2.99
464690 rs146191002 8_70 1.000 559.76 0.00180 -0.15
116755 rs71397673 2_102 1.000 475.72 0.00150 34.50
166478 rs7619398 3_57 0.515 844.12 0.00140 -3.08
464681 rs72681356 8_70 0.777 555.48 0.00140 4.78
166477 rs1436648 3_57 0.503 843.92 0.00130 -3.09
396107 rs917793 7_32 0.581 701.86 0.00130 33.28
396094 rs10259649 7_32 0.949 395.43 0.00120 28.65
191480 rs10653660 3_104 0.967 346.22 0.00110 -23.54
166479 rs2256473 3_57 0.367 843.48 0.00098 -3.06
551341 rs12244851 10_70 1.000 290.24 0.00092 16.46
166476 rs2575789 3_57 0.316 843.35 0.00085 -3.06
464683 rs72681364 8_70 0.486 554.02 0.00085 4.74
176589 rs72964564 3_76 1.000 233.33 0.00074 -16.88
116753 rs11396827 2_102 1.000 207.52 0.00066 20.35
75333 rs780093 2_16 1.000 185.69 0.00059 14.95
396111 rs2908282 7_32 0.258 699.79 0.00057 33.26
385925 rs1974619 7_15 0.605 239.45 0.00046 16.89
396082 rs17769733 7_32 0.986 131.06 0.00041 -7.76
116735 rs11676084 2_102 0.998 127.22 0.00040 -23.20
396101 rs4607517 7_32 0.160 699.76 0.00036 33.23
696339 rs35889227 14_45 1.000 114.94 0.00036 -11.34
54259 rs79687284 1_108 1.000 109.90 0.00035 12.08
647341 rs576674 13_10 1.000 103.82 0.00033 -10.82
323513 rs55792466 6_7 0.994 96.62 0.00030 -9.67
878465 rs10305492 6_30 1.000 94.37 0.00030 -10.36
513273 rs115478735 9_70 0.996 81.46 0.00026 9.52
396102 rs138917529 7_32 1.000 77.85 0.00025 -10.45
645651 rs60353775 13_7 0.967 81.75 0.00025 9.78
573597 rs2596407 11_29 1.000 74.93 0.00024 9.98
479647 rs10758593 9_4 0.999 70.84 0.00022 8.64
293293 rs17085655 5_57 0.546 118.06 0.00020 -11.48
54268 rs3754140 1_108 0.998 58.39 0.00019 -10.01
81659 rs2121564 2_28 1.000 59.62 0.00019 8.01
550793 rs11195508 10_70 0.986 60.57 0.00019 -10.76
323494 rs76623841 6_7 1.000 56.94 0.00018 -6.80
385923 rs10228796 7_15 0.222 236.94 0.00017 16.82
551335 rs117764423 10_70 0.725 75.34 0.00017 -5.24
561396 rs3750952 11_6 1.000 52.32 0.00017 -7.40
891663 rs3842753 11_2 0.487 106.86 0.00017 -9.28
891668 rs689 11_2 0.508 106.96 0.00017 -9.28
468351 rs4433184 8_78 0.959 51.97 0.00016 4.78
114461 rs12692596 2_96 1.000 45.75 0.00015 7.24
261918 rs12499544 4_119 0.720 67.46 0.00015 -8.39
328105 rs2206734 6_15 0.744 64.61 0.00015 8.26
468375 rs28529793 8_78 0.687 69.35 0.00015 -6.80
514563 rs28624681 9_73 0.940 51.48 0.00015 7.55
893567 rs231362 11_2 0.969 47.55 0.00015 6.83SNPs 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
116763 rs853789 2_102 1.000 902.00 2.9e-03 53.04
116756 rs537183 2_102 1.000 874.00 2.8e-03 52.66
116757 rs518598 2_102 0.000 827.73 3.9e-08 52.13
116759 rs485094 2_102 0.000 724.00 1.3e-10 51.04
116761 rs2544360 2_102 0.000 569.09 1.6e-10 39.60
116762 rs853791 2_102 0.000 539.25 7.4e-11 39.21
116765 rs853785 2_102 0.000 409.63 1.3e-11 37.45
116764 rs860510 2_102 0.000 397.15 1.3e-11 37.03
116758 rs579275 2_102 0.000 412.80 1.3e-11 36.60
116755 rs71397673 2_102 1.000 475.72 1.5e-03 34.50
396107 rs917793 7_32 0.581 701.86 1.3e-03 33.28
396111 rs2908282 7_32 0.258 699.79 5.7e-04 33.26
396101 rs4607517 7_32 0.160 699.76 3.6e-04 33.23
396113 rs732360 7_32 0.000 688.87 9.2e-07 32.89
396094 rs10259649 7_32 0.949 395.43 1.2e-03 28.65
396092 rs2908294 7_32 0.051 385.68 6.2e-05 28.27
116749 rs62176784 2_102 0.001 136.26 3.8e-07 -25.09
116751 rs549410 2_102 0.000 189.35 3.4e-08 -23.64
191480 rs10653660 3_104 0.967 346.22 1.1e-03 -23.54
191490 rs8192675 3_104 0.037 339.11 4.0e-05 -23.31
116735 rs11676084 2_102 0.998 127.22 4.0e-04 -23.20
116753 rs11396827 2_102 1.000 207.52 6.6e-04 20.35
116736 rs2140046 2_102 0.000 66.29 3.2e-09 -18.78
191488 rs11920090 3_104 0.229 160.37 1.2e-04 -18.33
191474 rs12492910 3_104 0.179 159.41 9.1e-05 -18.32
191487 rs11923694 3_104 0.158 158.92 8.0e-05 -18.31
191477 rs12496506 3_104 0.168 159.03 8.5e-05 -18.30
191491 rs11928798 3_104 0.118 157.40 5.9e-05 -18.26
191492 rs6785803 3_104 0.121 157.40 6.1e-05 -18.25
191467 rs56351320 3_104 0.002 210.96 1.3e-06 -17.54
191452 rs6792607 3_104 0.006 198.65 4.0e-06 -17.19
551349 rs12260037 10_70 0.000 249.92 1.0e-09 17.19
116754 rs13430620 2_102 0.000 89.87 7.2e-11 -16.98
385925 rs1974619 7_15 0.605 239.45 4.6e-04 16.89
176589 rs72964564 3_76 1.000 233.33 7.4e-04 -16.88
385923 rs10228796 7_15 0.222 236.94 1.7e-04 16.82
385924 rs2191349 7_15 0.174 236.46 1.3e-04 16.80
116760 rs114932341 2_102 0.000 184.89 3.7e-10 16.61
385926 rs188745922 7_15 0.011 230.65 7.8e-06 16.59
551341 rs12244851 10_70 1.000 290.24 9.2e-04 16.46
191444 rs11919048 3_104 0.013 135.07 5.6e-06 -16.13
176587 rs34642857 3_76 0.007 212.52 5.0e-06 -16.11
191459 rs79560566 3_104 0.001 126.20 4.4e-07 -15.48
385921 rs6461153 7_15 0.003 203.80 1.9e-06 15.42
385920 rs10266209 7_15 0.003 203.43 1.8e-06 15.41
385919 rs10249299 7_15 0.005 204.26 3.2e-06 -15.31
75333 rs780093 2_16 1.000 185.69 5.9e-04 14.95
191482 rs143791579 3_104 0.001 96.41 3.4e-07 -14.91
385914 rs4721398 7_15 0.003 187.44 1.5e-06 -14.88
191475 rs73167792 3_104 0.001 93.46 3.5e-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] 6if (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 SCF complex assembly (GO:0010265)
2 unsaturated fatty acid biosynthetic process (GO:0006636)
3 protein poly-ADP-ribosylation (GO:0070212)
4 protein auto-ADP-ribosylation (GO:0070213)
5 inorganic ion transmembrane transport (GO:0098660)
6 carboxylic acid biosynthetic process (GO:0046394)
7 alpha-linolenic acid metabolic process (GO:0036109)
8 icosanoid biosynthetic process (GO:0046456)
9 protein localization to chromosome, telomeric region (GO:0070198)
10 negative regulation of telomere maintenance (GO:0032205)
11 positive regulation of telomere capping (GO:1904355)
12 positive regulation of cyclin-dependent protein serine/threonine kinase activity (GO:0045737)
13 regulation of telomere maintenance via telomere lengthening (GO:1904356)
14 positive regulation of cyclin-dependent protein kinase activity (GO:1904031)
15 linoleic acid metabolic process (GO:0043651)
16 regulation of telomere capping (GO:1904353)
17 positive regulation of G1/S transition of mitotic cell cycle (GO:1900087)
18 negative regulation of telomere maintenance via telomere lengthening (GO:1904357)
19 protein ADP-ribosylation (GO:0006471)
20 regulation of DNA biosynthetic process (GO:2000278)
21 long-chain fatty acid biosynthetic process (GO:0042759)
22 negative regulation of DNA binding (GO:0043392)
23 positive regulation of telomere maintenance (GO:0032206)
24 protein ubiquitination (GO:0016567)
25 positive regulation of telomere maintenance via telomerase (GO:0032212)
26 positive regulation of telomerase activity (GO:0051973)
27 protein localization to chromosome (GO:0034502)
28 positive regulation of cell cycle G1/S phase transition (GO:1902808)
29 positive regulation of telomere maintenance via telomere lengthening (GO:1904358)
30 phospholipid biosynthetic process (GO:0008654)
31 organophosphate biosynthetic process (GO:0090407)
32 cation transmembrane transport (GO:0098655)
33 regulation of telomerase activity (GO:0051972)
34 regulation of telomere maintenance via telomerase (GO:0032210)
35 regulation of cyclin-dependent protein kinase activity (GO:1904029)
36 unsaturated fatty acid metabolic process (GO:0033559)
37 icosanoid metabolic process (GO:0006690)
38 positive regulation of mitotic cell cycle phase transition (GO:1901992)
39 peptidyl-threonine phosphorylation (GO:0018107)
40 positive regulation of DNA biosynthetic process (GO:2000573)
41 cellular response to nutrient levels (GO:0031669)
42 positive regulation of cell cycle (GO:0045787)
43 peptidyl-threonine modification (GO:0018210)
44 regulation of G1/S transition of mitotic cell cycle (GO:2000045)
45 fatty acid biosynthetic process (GO:0006633)
46 regulation of peptide hormone secretion (GO:0090276)
47 phospholipid metabolic process (GO:0006644)
48 lipid biosynthetic process (GO:0008610)
49 regulation of cyclin-dependent protein serine/threonine kinase activity (GO:0000079)
50 long-chain fatty acid metabolic process (GO:0001676)
51 spindle assembly (GO:0051225)
Overlap Adjusted.P.value Genes
1 1/6 0.03462351 CAND2
2 1/9 0.03462351 FADS1
3 1/10 0.03462351 TNKS
4 1/12 0.03462351 TNKS
5 1/12 0.03462351 ABCC8
6 1/13 0.03462351 FADS1
7 1/13 0.03462351 FADS1
8 1/15 0.03462351 FADS1
9 1/15 0.03462351 TNKS
10 1/16 0.03462351 TNKS
11 1/17 0.03462351 TNKS
12 1/17 0.03462351 CCND2
13 1/18 0.03462351 TNKS
14 1/20 0.03462351 CCND2
15 1/21 0.03462351 FADS1
16 1/23 0.03462351 TNKS
17 1/26 0.03462351 CCND2
18 1/26 0.03462351 TNKS
19 1/29 0.03462351 TNKS
20 1/29 0.03462351 TNKS
21 1/30 0.03462351 FADS1
22 1/31 0.03462351 TNKS
23 1/32 0.03462351 TNKS
24 2/525 0.03462351 CAND2;TNKS
25 1/33 0.03462351 TNKS
26 1/34 0.03462351 TNKS
27 1/34 0.03462351 TNKS
28 1/35 0.03462351 CCND2
29 1/36 0.03462351 TNKS
30 1/37 0.03462351 FADS1
31 1/39 0.03530897 FADS1
32 1/44 0.03856680 ABCC8
33 1/46 0.03908826 TNKS
34 1/52 0.04202037 TNKS
35 1/54 0.04202037 CCND2
36 1/54 0.04202037 FADS1
37 1/57 0.04268337 FADS1
38 1/58 0.04268337 CCND2
39 1/60 0.04268337 TNKS
40 1/61 0.04268337 TNKS
41 1/66 0.04357825 FADS1
42 1/66 0.04357825 CCND2
43 1/67 0.04357825 TNKS
44 1/71 0.04410546 CCND2
45 1/71 0.04410546 FADS1
46 1/74 0.04495290 ABCC8
47 1/76 0.04517426 FADS1
48 1/80 0.04596746 FADS1
49 1/82 0.04596746 CCND2
50 1/83 0.04596746 FADS1
51 1/84 0.04596746 TNKS
[1] "GO_Cellular_Component_2021"
Term Overlap
1 nuclear membrane (GO:0031965) 2/204
2 cation-transporting ATPase complex (GO:0090533) 1/16
3 cyclin-dependent protein kinase holoenzyme complex (GO:0000307) 1/30
4 serine/threonine protein kinase complex (GO:1902554) 1/37
Adjusted.P.value Genes
1 0.01813961 CCND2;TNKS
2 0.02874565 ABCC8
3 0.03315016 CCND2
4 0.03315016 CCND2
[1] "GO_Molecular_Function_2021"
Term
1 ATP-activated inward rectifier potassium channel activity (GO:0015272)
2 zinc ion binding (GO:0008270)
3 transition metal ion binding (GO:0046914)
4 inward rectifier potassium channel activity (GO:0005242)
5 NAD+ ADP-ribosyltransferase activity (GO:0003950)
6 potassium ion transmembrane transporter activity (GO:0015079)
7 pentosyltransferase activity (GO:0016763)
8 cyclin-dependent protein serine/threonine kinase regulator activity (GO:0016538)
9 potassium channel activity (GO:0005267)
10 protein kinase regulator activity (GO:0019887)
11 cation channel activity (GO:0005261)
Overlap Adjusted.P.value Genes
1 1/6 0.02021649 ABCC8
2 2/336 0.02021649 TNKS;QPCTL
3 2/445 0.02021649 TNKS;QPCTL
4 1/25 0.02021649 ABCC8
5 1/26 0.02021649 TNKS
6 1/32 0.02071934 ABCC8
7 1/38 0.02107353 TNKS
8 1/44 0.02133483 CCND2
9 1/91 0.03432662 ABCC8
10 1/98 0.03432662 CCND2
11 1/98 0.03432662 ABCC8QPCTL gene(s) from the input list not found in DisGeNET CURATEDTNKS gene(s) from the input list not found in DisGeNET CURATED Description
6 Colorectal Neoplasms
31 Hypoglycemia, leucine-induced
48 DIABETES MELLITUS, TRANSIENT NEONATAL, 2 (disorder)
67 MEGALENCEPHALY-POLYMICROGYRIA-POLYDACTYLY-HYDROCEPHALUS SYNDROME 3
68 Autosomal dominant hyperinsulinism due to SUR1 deficiency
5 Colorectal Carcinoma
26 POLYDACTYLY, POSTAXIAL
49 Perisylvian syndrome
50 Developmental Delay, Epilepsy, and Neonatal Diabetes
51 Megalanecephaly Polymicrogyria-Polydactyly Hydrocephalus Syndrome
FDR Ratio BgRatio
6 0.005936920 3/4 277/9703
31 0.005936920 1/4 1/9703
48 0.005936920 1/4 1/9703
67 0.005936920 1/4 1/9703
68 0.005936920 1/4 1/9703
5 0.006981373 3/4 702/9703
26 0.006981373 1/4 4/9703
49 0.006981373 1/4 4/9703
50 0.006981373 1/4 2/9703
51 0.006981373 1/4 4/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