Calcium (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
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-30680_irnt_Liver.Rmd) and HTML (docs/ukb-d-30680_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 Calcium (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-30680_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.0065371853 0.0002032958 #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
100.05597 12.19469 #report sample size
print(sample_size)[1] 315153#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.02262449 0.06841696 #compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))[1] 0.2292675 1.0630698Genes with highest PIPs
#distribution of PIPs
hist(ctwas_gene_res$susie_pip, xlim=c(0,1), main="Distribution of Gene PIPs")
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
#genes with PIP>0.8 or 20 highest PIPs
head(ctwas_gene_res[order(-ctwas_gene_res$susie_pip),report_cols], max(sum(ctwas_gene_res$susie_pip>0.8), 20)) genename region_tag susie_pip mu2 PVE z
1144 ASAP3 1_16 1.000 72.22 2.3e-04 8.43
2914 WDR75 2_113 1.000 1984.70 6.3e-03 0.15
5389 RPS11 19_34 1.000 66703.87 2.1e-01 9.19
2173 TMEM176B 7_93 0.993 80.23 2.5e-04 -8.80
11948 RP11-254F7.1 2_6 0.992 103.72 3.3e-04 7.21
10384 NTRK1 1_77 0.990 60.64 1.9e-04 -7.56
9613 C14orf80 14_55 0.963 58.76 1.8e-04 7.36
10134 VKORC1L1 7_44 0.956 150.12 4.6e-04 12.26
8765 ZNF77 19_3 0.936 26.97 8.0e-05 4.91
2546 LTBR 12_7 0.914 27.94 8.1e-05 4.94
8531 TNKS 8_12 0.913 53.84 1.6e-04 9.62
5400 EPHA2 1_11 0.845 58.18 1.6e-04 -7.86
10763 NYNRIN 14_3 0.822 40.95 1.1e-04 6.15
12467 RP11-219B17.3 15_27 0.819 30.83 8.0e-05 -5.52
9855 PALM3 19_11 0.808 30.07 7.7e-05 -5.03
2474 CBL 11_71 0.794 46.46 1.2e-04 -7.01
12074 RP11-131K5.2 17_12 0.794 27.15 6.8e-05 -4.73
9588 IFITM2 11_1 0.786 24.46 6.1e-05 -4.06
8803 DLEU1 13_21 0.777 29.07 7.2e-05 5.01
11389 GATSL3 22_10 0.775 34.02 8.4e-05 -5.85Genes with largest effect sizes
#plot PIP vs effect size
plot(ctwas_gene_res$susie_pip, ctwas_gene_res$mu2, xlab="PIP", ylab="mu^2", main="Gene PIPs vs Effect Size")
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
#genes with 20 largest effect sizes
head(ctwas_gene_res[order(-ctwas_gene_res$mu2),report_cols],20) genename region_tag susie_pip mu2 PVE z
5389 RPS11 19_34 1 66703.87 2.1e-01 9.19
1227 FLT3LG 19_34 0 57617.36 0.0e+00 -7.94
5393 RCN3 19_34 0 21431.87 0.0e+00 -5.74
1931 FCGRT 19_34 0 19612.29 0.0e+00 -1.88
4556 TMEM60 7_49 0 17872.48 0.0e+00 8.56
3804 PRRG2 19_34 0 9641.17 0.0e+00 -7.38
3803 PRMT1 19_34 0 6530.09 0.0e+00 -2.46
3805 SCAF1 19_34 0 6514.17 0.0e+00 -4.51
3802 IRF3 19_34 0 6315.07 0.0e+00 -4.00
10674 CASC10 10_16 0 6066.44 0.0e+00 -4.77
1940 SLC17A7 19_34 0 4625.84 0.0e+00 -4.14
10903 APTR 7_49 0 3375.72 0.0e+00 2.11
3411 CSTA 3_76 0 2324.66 1.9e-15 -45.63
1932 PIH1D1 19_34 0 2020.63 0.0e+00 3.44
2914 WDR75 2_113 1 1984.70 6.3e-03 0.15
9811 RSBN1L 7_49 0 1893.42 0.0e+00 2.19
10291 CTC-301O7.4 19_34 0 1592.22 0.0e+00 1.40
92 PHTF2 7_49 0 1363.23 0.0e+00 -0.24
11199 LINC00271 6_89 0 1303.72 0.0e+00 -3.84
2794 KPNA1 3_76 0 988.90 1.4e-14 -29.60Genes with highest PVE
#genes with 20 highest pve
head(ctwas_gene_res[order(-ctwas_gene_res$PVE),report_cols],20) genename region_tag susie_pip mu2 PVE z
5389 RPS11 19_34 1.000 66703.87 2.1e-01 9.19
2914 WDR75 2_113 1.000 1984.70 6.3e-03 0.15
10134 VKORC1L1 7_44 0.956 150.12 4.6e-04 12.26
7656 CATSPER2 15_16 0.774 151.31 3.7e-04 12.17
11948 RP11-254F7.1 2_6 0.992 103.72 3.3e-04 7.21
2173 TMEM176B 7_93 0.993 80.23 2.5e-04 -8.80
1144 ASAP3 1_16 1.000 72.22 2.3e-04 8.43
10384 NTRK1 1_77 0.990 60.64 1.9e-04 -7.56
9613 C14orf80 14_55 0.963 58.76 1.8e-04 7.36
5400 EPHA2 1_11 0.845 58.18 1.6e-04 -7.86
8531 TNKS 8_12 0.913 53.84 1.6e-04 9.62
10734 NAP1L4 11_2 0.631 74.95 1.5e-04 -9.36
12520 RP11-471B22.3 14_45 0.754 53.09 1.3e-04 -6.90
2474 CBL 11_71 0.794 46.46 1.2e-04 -7.01
2536 SH2B3 12_67 0.591 58.86 1.1e-04 -3.77
10763 NYNRIN 14_3 0.822 40.95 1.1e-04 6.15
4862 GGH 8_48 0.667 47.65 1.0e-04 6.56
9478 KMT5A 12_75 0.535 54.89 9.3e-05 -5.43
7915 GLYCTK 3_36 0.735 37.41 8.7e-05 -5.42
879 DGKD 2_137 0.637 41.96 8.5e-05 1.05Genes 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
3411 CSTA 3_76 0.000 2324.66 1.9e-15 -45.63
4999 PARP9 3_76 0.000 932.61 3.6e-18 30.38
2794 KPNA1 3_76 0.000 988.90 1.4e-14 -29.60
10167 WDR5B 3_76 0.000 141.47 0.0e+00 -13.97
5674 AHSG 3_114 0.000 96.06 1.1e-07 -12.38
10134 VKORC1L1 7_44 0.956 150.12 4.6e-04 12.26
7656 CATSPER2 15_16 0.774 151.31 3.7e-04 12.17
2887 NRBP1 2_16 0.004 139.41 1.9e-06 11.82
6480 CD109 6_51 0.004 134.52 1.6e-06 -11.69
6734 CCDC58 3_76 0.000 207.60 0.0e+00 11.20
8040 THBS3 1_76 0.037 85.05 1.0e-05 10.27
11738 RP11-115J16.2 8_12 0.003 87.79 8.3e-07 10.18
8284 RBKS 2_16 0.014 102.98 4.6e-06 9.93
2891 SNX17 2_16 0.024 109.54 8.2e-06 -9.75
8531 TNKS 8_12 0.913 53.84 1.6e-04 9.62
1058 GCKR 2_16 0.082 102.29 2.7e-05 -9.51
10987 C2orf16 2_16 0.082 102.29 2.7e-05 -9.51
8041 SLC50A1 1_76 0.008 69.34 1.7e-06 -9.37
10734 NAP1L4 11_2 0.631 74.95 1.5e-04 -9.36
3731 MED1 17_23 0.044 88.58 1.2e-05 9.25Comparing 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.0177048#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
3411 CSTA 3_76 0.000 2324.66 1.9e-15 -45.63
4999 PARP9 3_76 0.000 932.61 3.6e-18 30.38
2794 KPNA1 3_76 0.000 988.90 1.4e-14 -29.60
10167 WDR5B 3_76 0.000 141.47 0.0e+00 -13.97
5674 AHSG 3_114 0.000 96.06 1.1e-07 -12.38
10134 VKORC1L1 7_44 0.956 150.12 4.6e-04 12.26
7656 CATSPER2 15_16 0.774 151.31 3.7e-04 12.17
2887 NRBP1 2_16 0.004 139.41 1.9e-06 11.82
6480 CD109 6_51 0.004 134.52 1.6e-06 -11.69
6734 CCDC58 3_76 0.000 207.60 0.0e+00 11.20
8040 THBS3 1_76 0.037 85.05 1.0e-05 10.27
11738 RP11-115J16.2 8_12 0.003 87.79 8.3e-07 10.18
8284 RBKS 2_16 0.014 102.98 4.6e-06 9.93
2891 SNX17 2_16 0.024 109.54 8.2e-06 -9.75
8531 TNKS 8_12 0.913 53.84 1.6e-04 9.62
1058 GCKR 2_16 0.082 102.29 2.7e-05 -9.51
10987 C2orf16 2_16 0.082 102.29 2.7e-05 -9.51
8041 SLC50A1 1_76 0.008 69.34 1.7e-06 -9.37
10734 NAP1L4 11_2 0.631 74.95 1.5e-04 -9.36
3731 MED1 17_23 0.044 88.58 1.2e-05 9.25Locus 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: 3_76"
genename region_tag susie_pip mu2 PVE z
3411 CSTA 3_76 0 2324.66 1.9e-15 -45.63
2792 FAM162A 3_76 0 30.22 0.0e+00 2.40
6734 CCDC58 3_76 0 207.60 0.0e+00 11.20
10167 WDR5B 3_76 0 141.47 0.0e+00 -13.97
2794 KPNA1 3_76 0 988.90 1.4e-14 -29.60
4999 PARP9 3_76 0 932.61 3.6e-18 30.38
8518 PARP15 3_76 0 82.00 3.2e-19 -6.65
8517 PARP14 3_76 0 13.66 0.0e+00 1.59
8023 HSPBAP1 3_76 0 15.41 0.0e+00 3.28
4996 DIRC2 3_76 0 18.86 0.0e+00 3.96
12407 LINC02035 3_76 0 23.84 0.0e+00 3.53
1008 SEMA5B 3_76 0 33.59 0.0e+00 -5.44
3410 SEC22A 3_76 0 8.89 0.0e+00 -1.07
10775 HACD2 3_76 0 27.98 0.0e+00 -2.24
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 3_114"
genename region_tag susie_pip mu2 PVE z
7194 LIPH 3_114 0.000 6.30 7.5e-09 0.55
7196 SENP2 3_114 0.000 5.08 5.0e-09 0.17
11529 ETV5 3_114 0.001 16.18 4.9e-08 1.51
506 DGKG 3_114 0.001 14.60 3.9e-08 -1.34
10806 CRYGS 3_114 0.001 12.36 2.0e-08 -2.11
2784 TBCCD1 3_114 0.001 10.45 2.3e-08 0.14
5674 AHSG 3_114 0.000 96.06 1.1e-07 -12.38
2786 HRG 3_114 0.001 20.25 4.7e-08 -3.19
6503 EIF4A2 3_114 0.000 7.11 7.8e-09 1.33
7199 RFC4 3_114 0.000 6.29 6.4e-09 1.07
11225 LINC02043 3_114 0.000 5.04 5.0e-09 0.43
800 ST6GAL1 3_114 0.000 8.25 1.2e-08 0.03
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 7_44"
genename region_tag susie_pip mu2 PVE z
11283 ZNF736 7_44 0.002 9.22 4.4e-08 -0.19
10043 ZNF107 7_44 0.043 35.52 4.8e-06 0.51
10170 ZNF138 7_44 0.002 12.52 9.4e-08 -0.37
10333 ZNF273 7_44 0.001 6.50 1.8e-08 1.25
6235 ZNF117 7_44 0.009 30.49 8.6e-07 -3.21
10829 ERV3-1 7_44 0.003 27.66 2.8e-07 -3.90
5818 ZNF92 7_44 0.002 12.94 9.9e-08 -0.17
12364 RP11-479O9.4 7_44 0.010 53.18 1.7e-06 -6.09
10134 VKORC1L1 7_44 0.956 150.12 4.6e-04 12.26
8123 GUSB 7_44 0.001 11.71 3.3e-08 -2.33
11438 CRCP 7_44 0.001 8.79 2.4e-08 -1.96
8118 TPST1 7_44 0.001 6.85 2.1e-08 -0.60
11355 GS1-124K5.4 7_44 0.001 63.88 2.5e-07 -7.80
11493 KCTD7 7_44 0.014 42.80 1.9e-06 4.51
6356 RABGEF1 7_44 0.001 9.25 2.6e-08 -2.00
12416 RP11-458F8.4 7_44 0.001 48.67 2.1e-07 6.67
2180 TMEM248 7_44 0.001 18.07 7.2e-08 -3.04
10490 TYW1 7_44 0.002 31.84 2.5e-07 -4.73
12444 RP11-166O4.6 7_44 0.001 6.67 2.3e-08 0.21
11321 LINC01372 7_44 0.001 9.18 2.6e-08 2.23
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 15_16"
genename region_tag susie_pip mu2 PVE z
1853 ZNF106 15_16 0.004 12.45 1.4e-07 -2.68
9202 LRRC57 15_16 0.008 12.48 3.3e-07 0.97
6691 STARD9 15_16 0.007 10.81 2.3e-07 -0.94
5189 CDAN1 15_16 0.004 5.18 5.9e-08 0.36
3962 TTBK2 15_16 0.004 6.85 8.0e-08 1.42
4903 TMEM62 15_16 0.005 7.28 1.1e-07 0.58
7984 ADAL 15_16 0.004 27.22 3.1e-07 -4.73
7985 LCMT2 15_16 0.004 81.12 9.6e-07 -8.76
4898 TUBGCP4 15_16 0.004 32.67 3.7e-07 5.27
5180 ZSCAN29 15_16 0.004 8.38 1.2e-07 -1.30
7702 MAP1A 15_16 0.005 60.81 8.8e-07 7.38
7656 CATSPER2 15_16 0.774 151.31 3.7e-04 12.17
7709 PDIA3 15_16 0.004 17.71 2.4e-07 3.44
5178 MFAP1 15_16 0.004 32.42 4.4e-07 5.10
1286 WDR76 15_16 0.004 29.94 3.8e-07 4.94
| 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.004 6.24 7.6e-08 1.32
11149 OST4 2_16 0.021 20.71 1.4e-06 2.06
4939 EMILIN1 2_16 0.004 10.00 1.2e-07 -2.73
4927 KHK 2_16 0.012 16.86 6.6e-07 2.25
4935 PREB 2_16 0.004 13.90 1.7e-07 3.36
4941 ATRAID 2_16 0.009 68.09 1.9e-06 -7.86
4936 SLC5A6 2_16 0.010 70.51 2.2e-06 7.98
1060 CAD 2_16 0.005 41.51 6.7e-07 6.05
2885 SLC30A3 2_16 0.004 28.81 3.7e-07 4.88
7169 UCN 2_16 0.005 9.59 1.5e-07 1.83
2891 SNX17 2_16 0.024 109.54 8.2e-06 -9.75
7170 ZNF513 2_16 0.004 29.10 3.7e-07 4.97
2887 NRBP1 2_16 0.004 139.41 1.9e-06 11.82
4925 IFT172 2_16 0.004 16.13 2.1e-07 3.43
1058 GCKR 2_16 0.082 102.29 2.7e-05 -9.51
10987 C2orf16 2_16 0.082 102.29 2.7e-05 -9.51
10407 GPN1 2_16 0.004 36.30 4.4e-07 5.52
8847 CCDC121 2_16 0.005 17.07 2.8e-07 -3.26
6575 BRE 2_16 0.005 16.58 2.6e-07 -3.34
8284 RBKS 2_16 0.014 102.98 4.6e-06 9.93
| 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
102 rs12135382 1_1 1.000 66.40 2.1e-04 8.54
13505 rs3176461 1_32 1.000 130.43 4.1e-04 12.04
55307 rs708774 1_116 1.000 62.97 2.0e-04 -6.93
71848 rs780093 2_16 1.000 279.28 8.9e-04 -18.36
72149 rs569546056 2_17 1.000 583.92 1.9e-03 2.86
92632 rs6730773 2_57 1.000 114.62 3.6e-04 -15.10
130537 rs142215640 2_136 1.000 229.93 7.3e-04 3.30
131202 rs7584554 2_137 1.000 206.82 6.6e-04 -16.92
167161 rs16853539 3_67 1.000 65.42 2.1e-04 9.87
171231 rs492301 3_74 1.000 34.37 1.1e-04 -5.36
171649 rs7638322 3_75 1.000 158.97 5.0e-04 16.77
171668 rs35829261 3_75 1.000 243.40 7.7e-04 -20.36
181906 rs9817452 3_97 1.000 42.79 1.4e-04 6.89
210540 rs16852326 4_33 1.000 91.35 2.9e-04 -10.22
229935 rs35518360 4_67 1.000 115.29 3.7e-04 -11.52
230001 rs13140033 4_68 1.000 59.51 1.9e-04 -7.99
313497 rs11746728 5_105 1.000 43.68 1.4e-04 6.63
362289 rs1763519 6_89 1.000 139.81 4.4e-04 -12.43
362737 rs199804242 6_89 1.000 4903.70 1.6e-02 2.91
362996 rs7750197 6_90 1.000 48.10 1.5e-04 -7.18
370704 rs60425481 6_104 1.000 73.86 2.3e-04 -4.98
399876 rs761767938 7_49 1.000 17323.74 5.5e-02 -8.61
399884 rs1544459 7_49 1.000 17453.28 5.5e-02 -8.52
488666 rs11144105 9_35 1.000 61.70 2.0e-04 -8.09
495141 rs7873928 9_48 1.000 55.91 1.8e-04 7.48
505324 rs62578106 9_65 1.000 74.42 2.4e-04 2.12
512826 rs17486892 10_9 1.000 121.90 3.9e-04 -13.30
512835 rs12219116 10_9 1.000 103.37 3.3e-04 7.30
517641 rs72806996 10_16 1.000 6877.51 2.2e-02 4.97
517644 rs72810731 10_16 1.000 6826.97 2.2e-02 5.32
517645 rs551737161 10_16 1.000 7084.21 2.2e-02 4.97
592812 rs4937122 11_77 1.000 143.82 4.6e-04 -12.78
597326 rs117213754 12_4 1.000 109.31 3.5e-04 10.16
599622 rs11616030 12_11 1.000 47.14 1.5e-04 -7.01
612650 rs7303812 12_33 1.000 35.66 1.1e-04 -5.93
633403 rs116736246 12_75 1.000 59.72 1.9e-04 -7.57
690982 rs7155616 14_46 1.000 57.99 1.8e-04 -7.73
691854 rs12893029 14_48 1.000 97.61 3.1e-04 -1.16
691870 rs1243165 14_48 1.000 60.87 1.9e-04 3.06
737989 rs12445932 16_39 1.000 50.31 1.6e-04 -7.60
741885 rs57960111 16_45 1.000 100.72 3.2e-04 -3.80
745610 rs7206699 16_53 1.000 135.85 4.3e-04 -12.56
747062 rs11078597 17_2 1.000 223.19 7.1e-04 16.05
748958 rs218676 17_6 1.000 35.87 1.1e-04 5.15
748963 rs189551802 17_6 1.000 35.94 1.1e-04 5.38
761238 rs2097730 17_33 1.000 60.48 1.9e-04 7.99
804510 rs4806075 19_24 1.000 150.10 4.8e-04 -3.94
804703 rs2251125 19_24 1.000 43.57 1.4e-04 -6.83
805059 rs149366150 19_25 1.000 40.11 1.3e-04 -6.51
805461 rs35496032 19_26 1.000 67.35 2.1e-04 -8.12
813920 rs74273659 20_5 1.000 48.22 1.5e-04 7.23
823270 rs1535025 20_24 1.000 51.57 1.6e-04 3.42
823272 rs4812458 20_24 1.000 62.17 2.0e-04 -9.18
823277 rs4812465 20_24 1.000 66.92 2.1e-04 -8.54
827437 rs6123359 20_32 1.000 52.03 1.7e-04 10.54
827439 rs209955 20_32 1.000 179.94 5.7e-04 14.80
827443 rs2585441 20_32 1.000 110.69 3.5e-04 11.34
874608 rs11589479 1_76 1.000 60.09 1.9e-04 10.04
912872 rs1042636 3_76 1.000 223.80 7.1e-04 -20.49
912923 rs73186030 3_76 1.000 2455.64 7.8e-03 52.15
985756 rs190424317 19_4 1.000 170.24 5.4e-04 -15.02
996023 rs113176985 19_34 1.000 58687.29 1.9e-01 -9.54
996026 rs374141296 19_34 1.000 59391.87 1.9e-01 -8.82
6100 rs1780323 1_15 0.999 86.64 2.7e-04 9.23
13566 rs41292511 1_33 0.999 39.38 1.2e-04 2.79
53143 rs11117836 1_110 0.999 64.88 2.1e-04 8.72
196500 rs3748034 4_4 0.999 55.69 1.8e-04 9.71
196515 rs36205397 4_4 0.999 65.50 2.1e-04 13.50
512840 rs72786681 10_9 0.999 289.46 9.2e-04 16.27
517663 rs113349503 10_16 0.999 6910.86 2.2e-02 5.04
534864 rs1749850 10_51 0.999 65.01 2.1e-04 8.34
540005 rs17109928 10_60 0.999 34.52 1.1e-04 3.40
764143 rs113408695 17_39 0.999 38.62 1.2e-04 -6.98
851063 rs9611850 22_18 0.999 36.77 1.2e-04 6.15
196501 rs3752442 4_4 0.998 47.15 1.5e-04 -11.12
732785 rs10083818 16_29 0.998 39.77 1.3e-04 -6.25
804509 rs1688031 19_24 0.998 134.63 4.3e-04 14.08
808890 rs10408866 19_35 0.998 35.39 1.1e-04 4.99
130536 rs12478406 2_136 0.997 99.50 3.1e-04 3.02
134183 rs12619647 2_144 0.996 45.58 1.4e-04 -6.94
827441 rs78767971 20_32 0.996 63.04 2.0e-04 -5.48
746146 rs60990680 16_54 0.995 41.24 1.3e-04 7.34
187004 rs234043 3_106 0.994 29.10 9.2e-05 5.32
191221 rs62294588 3_114 0.993 56.97 1.8e-04 -8.01
136555 rs11712103 3_4 0.992 41.61 1.3e-04 -6.49
630023 rs141105880 12_67 0.991 37.38 1.2e-04 -5.29
318516 rs6597256 6_7 0.989 27.58 8.7e-05 -5.09
330832 rs1187117 6_28 0.989 121.93 3.8e-04 -13.37
116817 rs7575118 2_110 0.988 27.50 8.6e-05 -5.03
435275 rs12547987 8_23 0.988 41.58 1.3e-04 -6.60
52071 rs79687284 1_108 0.986 26.18 8.2e-05 4.90
124683 rs150868267 2_125 0.986 28.47 8.9e-05 5.30
130572 rs1667308 2_136 0.985 96.02 3.0e-04 5.05
747199 rs3760230 17_3 0.985 53.92 1.7e-04 -7.61
191094 rs2070633 3_114 0.983 111.76 3.5e-04 13.85
381846 rs7811500 7_16 0.982 29.81 9.3e-05 -5.14
449108 rs62508470 8_54 0.982 33.58 1.0e-04 5.76
53074 rs1502341 1_110 0.980 35.03 1.1e-04 6.13
370785 rs56393506 6_104 0.979 39.53 1.2e-04 -2.79
209860 rs11096949 4_31 0.977 29.22 9.1e-05 -5.45
613352 rs7485577 12_36 0.976 34.72 1.1e-04 -6.18
301192 rs4958244 5_80 0.975 48.10 1.5e-04 -7.23
435217 rs118162691 8_23 0.974 29.49 9.1e-05 -5.43
597342 rs11062925 12_4 0.974 31.14 9.6e-05 -3.52
808855 rs73050880 19_35 0.974 37.92 1.2e-04 5.43
486890 rs11143715 9_30 0.969 43.02 1.3e-04 -7.62
645931 rs2167059 13_16 0.968 31.51 9.7e-05 -3.62
565432 rs2785172 11_23 0.966 26.19 8.0e-05 -5.01
669590 rs1952554 14_2 0.964 33.74 1.0e-04 5.84
762035 rs11650989 17_36 0.963 39.29 1.2e-04 -7.94
435767 rs13253859 8_24 0.960 42.21 1.3e-04 6.58
711391 rs56357772 15_36 0.960 44.03 1.3e-04 -6.65
83645 rs13012253 2_39 0.959 24.86 7.6e-05 -4.81
152295 rs759228865 3_39 0.958 25.30 7.7e-05 4.92
797372 rs113622516 19_10 0.958 31.67 9.6e-05 -5.06
512830 rs1144645 10_9 0.953 111.80 3.4e-04 -7.02
98753 rs10864869 2_70 0.951 30.13 9.1e-05 5.59
623081 rs10648843 12_54 0.947 28.73 8.6e-05 -4.60
684577 rs1952196 14_32 0.946 28.07 8.4e-05 5.23
985889 rs34944502 19_4 0.944 175.99 5.3e-04 15.20
667260 rs4408410 13_59 0.932 38.33 1.1e-04 -6.27
406429 rs3197597 7_61 0.925 28.32 8.3e-05 -4.63
363835 rs540973884 6_92 0.920 23.93 7.0e-05 -4.53
547189 rs10736272 10_74 0.920 27.64 8.1e-05 -5.07
46764 rs3043084 1_98 0.919 33.68 9.8e-05 5.84
812568 rs6111987 20_2 0.915 24.00 7.0e-05 -4.51
205390 rs28570016 4_22 0.912 34.00 9.8e-05 -5.59
302225 rs145769851 5_82 0.912 24.78 7.2e-05 -4.62
171110 rs114594872 3_74 0.909 30.76 8.9e-05 -5.04
81771 rs6545412 2_36 0.908 42.75 1.2e-04 6.59
6029 rs35204330 1_14 0.907 75.82 2.2e-04 -8.31
13593 rs78767748 1_33 0.905 63.49 1.8e-04 -4.53
224396 rs190282180 4_57 0.900 23.90 6.8e-05 -4.64
71369 rs13409702 2_15 0.898 25.94 7.4e-05 -5.08
305840 rs4958272 5_89 0.898 25.06 7.1e-05 4.72
738002 rs62053230 16_39 0.896 29.93 8.5e-05 6.80
602326 rs66720652 12_15 0.894 25.93 7.4e-05 -4.71
569051 rs7123635 11_28 0.892 31.72 9.0e-05 -5.49
741867 rs11639652 16_45 0.889 66.85 1.9e-04 7.44
759093 rs7221118 17_29 0.889 41.85 1.2e-04 -6.51
277409 rs458036 5_33 0.888 30.59 8.6e-05 4.87
72152 rs4580350 2_17 0.885 579.58 1.6e-03 -3.12
798976 rs368707654 19_14 0.884 30.46 8.5e-05 5.33
319799 rs183745515 6_10 0.879 25.74 7.2e-05 4.62
512311 rs263418 10_8 0.879 60.01 1.7e-04 -8.69
171626 rs150804648 3_75 0.873 93.90 2.6e-04 9.29
897548 rs1448208 2_66 0.873 42.88 1.2e-04 -6.47
420933 rs10224210 7_94 0.871 40.70 1.1e-04 6.37
552457 rs75184896 10_84 0.870 24.75 6.8e-05 -4.49
92627 rs189255944 2_57 0.858 45.08 1.2e-04 8.34
540165 rs149794046 10_61 0.853 24.72 6.7e-05 4.84
78449 rs2176348 2_28 0.852 25.44 6.9e-05 4.60
495415 rs10993381 9_48 0.842 128.99 3.4e-04 -12.07
130543 rs6727510 2_136 0.834 123.35 3.3e-04 3.90
399880 rs11972122 7_49 0.832 16214.63 4.3e-02 -8.59
138101 rs747130 3_6 0.824 24.41 6.4e-05 -4.55
129728 rs11900497 2_135 0.821 25.94 6.8e-05 -4.68
703802 rs12437803 15_20 0.811 42.59 1.1e-04 -6.99
513967 rs10906148 10_10 0.805 27.27 7.0e-05 -4.99
169134 rs11929640 3_70 0.804 26.95 6.9e-05 2.98
453851 rs1044730 8_62 0.803 24.23 6.2e-05 4.57
92542 rs746864936 2_57 0.802 40.01 1.0e-04 -7.84
760711 rs4617927 17_32 0.801 24.30 6.2e-05 -4.38SNPs 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
996026 rs374141296 19_34 1 59391.87 1.9e-01 -8.82
996014 rs61371437 19_34 0 58819.91 0.0e+00 -9.35
996004 rs739349 19_34 0 58702.65 0.0e+00 -9.29
996005 rs756628 19_34 0 58701.65 0.0e+00 -9.29
996023 rs113176985 19_34 1 58687.29 1.9e-01 -9.54
996016 rs35295508 19_34 0 58615.14 0.0e+00 -9.42
996001 rs739347 19_34 0 58608.84 0.0e+00 -9.28
996002 rs2073614 19_34 0 58535.96 0.0e+00 -9.21
995997 rs4802613 19_34 0 58432.17 0.0e+00 -9.15
996030 rs2946865 19_34 0 58394.67 3.1e-11 -9.51
996021 rs73056069 19_34 0 58341.41 0.0e+00 -9.31
996007 rs2077300 19_34 0 58339.91 0.0e+00 -9.07
996018 rs2878354 19_34 0 58246.49 0.0e+00 -9.18
996011 rs73056059 19_34 0 58226.34 0.0e+00 -9.15
995995 rs10403394 19_34 0 58121.18 0.0e+00 -9.19
995996 rs17555056 19_34 0 58062.31 0.0e+00 -9.15
996031 rs60815603 19_34 0 57537.48 0.0e+00 -9.43
996034 rs1316885 19_34 0 57282.27 0.0e+00 -9.48
996039 rs2946863 19_34 0 57167.27 0.0e+00 -9.54
996036 rs60746284 19_34 0 57160.10 0.0e+00 -9.31
996032 rs35443645 19_34 0 57124.41 0.0e+00 -9.48
996012 rs73056062 19_34 0 56690.16 0.0e+00 -8.82
996042 rs553431297 19_34 0 55673.03 0.0e+00 -8.90
996025 rs112283514 19_34 0 55549.48 0.0e+00 -8.85
996027 rs11270139 19_34 0 55196.40 0.0e+00 -8.83
995992 rs10421294 19_34 0 52518.68 0.0e+00 -8.45
995991 rs8108175 19_34 0 52512.33 0.0e+00 -8.45
995984 rs59192944 19_34 0 52418.15 0.0e+00 -8.43
995990 rs1858742 19_34 0 52407.38 0.0e+00 -8.41
995981 rs55991145 19_34 0 52379.68 0.0e+00 -8.46
995976 rs3786567 19_34 0 52361.22 0.0e+00 -8.46
995972 rs2271952 19_34 0 52342.02 0.0e+00 -8.46
995975 rs4801801 19_34 0 52338.44 0.0e+00 -8.46
995971 rs2271953 19_34 0 52280.25 0.0e+00 -8.42
995973 rs2271951 19_34 0 52277.52 0.0e+00 -8.42
995962 rs60365978 19_34 0 52239.48 0.0e+00 -8.44
995968 rs4802612 19_34 0 52019.38 0.0e+00 -8.25
995978 rs2517977 19_34 0 51885.19 0.0e+00 -8.39
995965 rs55893003 19_34 0 51822.92 0.0e+00 -8.15
995957 rs55992104 19_34 0 50667.79 0.0e+00 -8.17
995951 rs60403475 19_34 0 50660.52 0.0e+00 -8.17
995954 rs4352151 19_34 0 50650.85 0.0e+00 -8.17
995948 rs11878448 19_34 0 50617.00 0.0e+00 -8.17
995942 rs9653100 19_34 0 50604.82 0.0e+00 -8.17
995938 rs4802611 19_34 0 50572.33 0.0e+00 -8.16
995930 rs7251338 19_34 0 50503.01 0.0e+00 -8.19
995929 rs59269605 19_34 0 50497.09 0.0e+00 -8.17
995950 rs1042120 19_34 0 50348.35 0.0e+00 -7.96
995946 rs113220577 19_34 0 50305.65 0.0e+00 -7.95
995940 rs9653118 19_34 0 50228.37 0.0e+00 -7.94SNPs 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
996023 rs113176985 19_34 1.000 58687.29 0.19000 -9.54
996026 rs374141296 19_34 1.000 59391.87 0.19000 -8.82
399876 rs761767938 7_49 1.000 17323.74 0.05500 -8.61
399884 rs1544459 7_49 1.000 17453.28 0.05500 -8.52
399880 rs11972122 7_49 0.832 16214.63 0.04300 -8.59
517641 rs72806996 10_16 1.000 6877.51 0.02200 4.97
517644 rs72810731 10_16 1.000 6826.97 0.02200 5.32
517645 rs551737161 10_16 1.000 7084.21 0.02200 4.97
517663 rs113349503 10_16 0.999 6910.86 0.02200 5.04
362737 rs199804242 6_89 1.000 4903.70 0.01600 2.91
362753 rs6923513 6_89 0.786 4905.74 0.01200 3.32
362736 rs2327654 6_89 0.710 4905.57 0.01100 3.32
399881 rs11406602 7_49 0.168 16190.15 0.00860 -8.51
912923 rs73186030 3_76 1.000 2455.64 0.00780 52.15
72149 rs569546056 2_17 1.000 583.92 0.00190 2.86
72152 rs4580350 2_17 0.885 579.58 0.00160 -3.12
171660 rs34863476 3_75 0.780 468.66 0.00120 22.39
903712 rs2304701 2_113 0.190 1724.17 0.00100 3.76
912941 rs17199211 3_76 0.660 461.79 0.00097 -10.60
512840 rs72786681 10_9 0.999 289.46 0.00092 16.27
72148 rs7562170 2_17 0.493 576.56 0.00090 3.04
71848 rs780093 2_16 1.000 279.28 0.00089 -18.36
171668 rs35829261 3_75 1.000 243.40 0.00077 -20.36
130537 rs142215640 2_136 1.000 229.93 0.00073 3.30
747062 rs11078597 17_2 1.000 223.19 0.00071 16.05
912872 rs1042636 3_76 1.000 223.80 0.00071 -20.49
131202 rs7584554 2_137 1.000 206.82 0.00066 -16.92
827439 rs209955 20_32 1.000 179.94 0.00057 14.80
903779 rs62187106 2_113 0.102 1719.87 0.00056 3.76
985756 rs190424317 19_4 1.000 170.24 0.00054 -15.02
985889 rs34944502 19_4 0.944 175.99 0.00053 15.20
171649 rs7638322 3_75 1.000 158.97 0.00050 16.77
913009 rs4678179 3_76 0.340 460.39 0.00050 -10.43
804510 rs4806075 19_24 1.000 150.10 0.00048 -3.94
903590 rs17271036 2_113 0.087 1723.69 0.00048 3.72
592812 rs4937122 11_77 1.000 143.82 0.00046 -12.78
912766 rs76947531 3_76 0.577 247.61 0.00045 -22.79
362289 rs1763519 6_89 1.000 139.81 0.00044 -12.43
745610 rs7206699 16_53 1.000 135.85 0.00043 -12.56
804509 rs1688031 19_24 0.998 134.63 0.00043 14.08
343697 rs7745364 6_51 0.538 243.36 0.00042 -16.75
13505 rs3176461 1_32 1.000 130.43 0.00041 12.04
171678 rs34395935 3_75 0.751 171.44 0.00041 17.11
512826 rs17486892 10_9 1.000 121.90 0.00039 -13.30
330832 rs1187117 6_28 0.989 121.93 0.00038 -13.37
229935 rs35518360 4_67 1.000 115.29 0.00037 -11.52
804507 rs2445819 19_24 0.464 252.88 0.00037 14.99
903569 rs150157254 2_113 0.067 1715.53 0.00037 3.77
903677 rs2289225 2_113 0.068 1723.37 0.00037 3.71
92632 rs6730773 2_57 1.000 114.62 0.00036 -15.10SNPs with largest z scores
#SNPs with 50 largest z scores
head(ctwas_snp_res[order(-abs(ctwas_snp_res$z)),report_cols_snps],50) id region_tag susie_pip mu2 PVE z
912923 rs73186030 3_76 1 2455.64 7.8e-03 52.15
912871 rs1801725 3_76 0 2295.15 8.1e-19 50.89
912830 rs73186015 3_76 0 2290.82 8.1e-19 50.86
912837 rs17251221 3_76 0 2290.69 8.1e-19 50.86
912791 rs34408666 3_76 0 2291.31 8.1e-19 50.84
912895 rs145063534 3_76 0 2215.68 7.8e-19 50.58
912960 rs55716378 3_76 0 2064.69 2.9e-18 46.06
912996 rs56139876 3_76 0 2033.73 2.9e-18 45.71
913001 rs55838400 3_76 0 2034.31 2.9e-18 45.71
912983 rs73186044 3_76 0 2032.94 2.1e-18 45.70
912990 rs5008830 3_76 0 2032.98 2.9e-18 45.70
913003 rs2001548 3_76 0 2033.33 2.9e-18 45.70
913010 rs73186048 3_76 0 2030.94 2.9e-18 45.68
913025 rs73186050 3_76 0 2028.09 3.6e-18 45.64
913032 rs73186053 3_76 0 2027.54 3.6e-18 45.64
913069 rs17265703 3_76 0 2024.65 3.6e-18 45.61
913043 rs111928357 3_76 0 2023.08 2.9e-18 45.60
913048 rs73186057 3_76 0 2022.67 2.9e-18 45.60
913075 rs73186060 3_76 0 2022.46 3.6e-18 45.59
913080 rs112779893 3_76 0 2022.05 2.8e-18 45.59
913088 rs56085624 3_76 0 2022.08 3.6e-18 45.59
913057 rs112527366 3_76 0 2021.19 2.8e-18 45.58
913087 rs55784797 3_76 0 2021.95 3.6e-18 45.58
913072 rs73186059 3_76 0 2021.05 2.8e-18 45.57
913040 rs2877603 3_76 0 2020.32 2.8e-18 45.56
913041 rs2332182 3_76 0 1836.21 6.5e-19 43.48
912946 rs16832956 3_76 0 1508.08 0.0e+00 43.13
913107 rs56098325 3_76 0 1762.50 7.5e-18 42.22
913127 rs55636509 3_76 0 1731.11 7.9e-18 41.88
913157 rs4491840 3_76 0 1719.30 7.9e-18 41.74
913153 rs2332179 3_76 0 1718.63 7.9e-18 41.73
913207 rs16833078 3_76 0 1716.55 7.9e-18 41.70
913213 rs16833080 3_76 0 1715.86 7.9e-18 41.70
913170 rs111289055 3_76 0 1714.30 6.6e-18 41.68
913169 rs66508963 3_76 0 1713.97 6.6e-18 41.67
913177 rs55744726 3_76 0 1713.65 6.6e-18 41.67
913171 rs9812341 3_76 0 1712.57 6.6e-18 41.66
913186 rs73186095 3_76 0 1713.08 6.6e-18 41.66
913201 rs67446552 3_76 0 1713.07 7.2e-18 41.66
913202 rs67706840 3_76 0 1713.07 7.2e-18 41.66
913179 rs6438725 3_76 0 1711.66 6.6e-18 41.65
913184 rs9834317 3_76 0 1711.56 6.6e-18 41.65
913194 rs9881201 3_76 0 1711.16 6.6e-18 41.64
913203 rs56236741 3_76 0 1711.07 6.6e-18 41.64
913104 rs34173813 3_76 0 1717.73 2.4e-18 41.58
913228 rs28630628 3_76 0 1559.83 6.0e-17 41.18
913229 rs67931976 3_76 0 1494.05 2.8e-17 40.74
913237 rs55780469 3_76 0 1450.55 1.9e-15 39.70
913281 rs73188388 3_76 0 1451.34 2.2e-15 39.70
913284 rs371655862 3_76 0 1454.32 3.4e-15 39.70Gene 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] 15if (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 Overlap
1 positive regulation of transferase activity (GO:0051347) 3/148
2 positive regulation of GTPase activity (GO:0043547) 3/214
Adjusted.P.value Genes
1 0.03300417 NTRK1;TNKS;EPHA2
2 0.04873205 NTRK1;ASAP3;EPHA2
[1] "GO_Cellular_Component_2021"
Term Overlap Adjusted.P.value
1 focal adhesion (GO:0005925) 3/387 0.03183948
2 cell-substrate junction (GO:0030055) 3/394 0.03183948
Genes
1 ASAP3;RPS11;EPHA2
2 ASAP3;RPS11;EPHA2
[1] "GO_Molecular_Function_2021"
Term
1 transmembrane receptor protein kinase activity (GO:0019199)
2 transmembrane receptor protein tyrosine kinase activity (GO:0004714)
3 protein tyrosine kinase activity (GO:0004713)
4 nerve growth factor binding (GO:0048406)
5 neurotrophin binding (GO:0043121)
6 transmembrane-ephrin receptor activity (GO:0005005)
7 ephrin receptor activity (GO:0005003)
Overlap Adjusted.P.value Genes
1 2/60 0.01132767 NTRK1;EPHA2
2 2/60 0.01132767 NTRK1;EPHA2
3 2/108 0.02340436 NTRK1;EPHA2
4 1/5 0.02340436 NTRK1
5 1/8 0.02992621 NTRK1
6 1/17 0.04528111 EPHA2
7 1/17 0.04528111 EPHA2C14orf80 gene(s) from the input list not found in DisGeNET CURATEDRP11-219B17.3 gene(s) from the input list not found in DisGeNET CURATEDTNKS gene(s) from the input list not found in DisGeNET CURATEDPALM3 gene(s) from the input list not found in DisGeNET CURATEDRP11-254F7.1 gene(s) from the input list not found in DisGeNET CURATEDNYNRIN gene(s) from the input list not found in DisGeNET CURATEDZNF77 gene(s) from the input list not found in DisGeNET CURATEDWDR75 gene(s) from the input list not found in DisGeNET CURATEDTMEM176B gene(s) from the input list not found in DisGeNET CURATEDRPS11 gene(s) from the input list not found in DisGeNET CURATEDVKORC1L1 gene(s) from the input list not found in DisGeNET CURATED Description FDR Ratio
55 CATARACT, POSTERIOR POLAR, 1 0.01216244 1/4
57 Age-related cortical cataract 0.01216244 1/4
1 Congenital Pain Insensitivity 0.01338953 1/4
7 Hereditary Sensory Autonomic Neuropathy, Type 1 0.01338953 1/4
8 Hereditary Sensory Autonomic Neuropathy, Type 2 0.01338953 1/4
9 HSAN Type IV 0.01338953 1/4
10 Hereditary Sensory Autonomic Neuropathy, Type 5 0.01338953 1/4
14 Neuralgia 0.01338953 1/4
16 Hereditary Sensory and Autonomic Neuropathies 0.01338953 1/4
18 Psychosis, Brief Reactive 0.01338953 1/4
BgRatio
55 1/9703
57 1/9703
1 7/9703
7 8/9703
8 5/9703
9 4/9703
10 7/9703
14 16/9703
16 4/9703
18 14/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