C-reactive protein (quantile) - Whole_Blood
wesleycrouse
2021-09-09
Last updated: 2021-09-09
Checks: 6 1
Knit directory: ctwas_applied/
This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
The R Markdown file has unstaged changes. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20210726) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 59e5f4d. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Unstaged changes:
Modified: analysis/ukb-d-30500_irnt_Liver.Rmd
Modified: analysis/ukb-d-30500_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30600_irnt_Liver.Rmd
Modified: analysis/ukb-d-30600_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30610_irnt_Liver.Rmd
Modified: analysis/ukb-d-30610_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30620_irnt_Liver.Rmd
Modified: analysis/ukb-d-30620_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30630_irnt_Liver.Rmd
Modified: analysis/ukb-d-30630_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30640_irnt_Liver.Rmd
Modified: analysis/ukb-d-30640_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30650_irnt_Liver.Rmd
Modified: analysis/ukb-d-30650_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30660_irnt_Liver.Rmd
Modified: analysis/ukb-d-30660_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30670_irnt_Liver.Rmd
Modified: analysis/ukb-d-30670_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30680_irnt_Liver.Rmd
Modified: analysis/ukb-d-30690_irnt_Liver.Rmd
Modified: analysis/ukb-d-30690_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30700_irnt_Liver.Rmd
Modified: analysis/ukb-d-30700_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30710_irnt_Liver.Rmd
Modified: analysis/ukb-d-30710_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30720_irnt_Liver.Rmd
Modified: analysis/ukb-d-30730_irnt_Liver.Rmd
Modified: analysis/ukb-d-30740_irnt_Liver.Rmd
Modified: analysis/ukb-d-30750_irnt_Liver.Rmd
Modified: analysis/ukb-d-30760_irnt_Liver.Rmd
Modified: analysis/ukb-d-30770_irnt_Liver.Rmd
Modified: analysis/ukb-d-30780_irnt_Liver.Rmd
Modified: analysis/ukb-d-30790_irnt_Liver.Rmd
Modified: analysis/ukb-d-30800_irnt_Liver.Rmd
Modified: analysis/ukb-d-30810_irnt_Liver.Rmd
Modified: analysis/ukb-d-30820_irnt_Liver.Rmd
Modified: analysis/ukb-d-30830_irnt_Liver.Rmd
Modified: analysis/ukb-d-30840_irnt_Liver.Rmd
Modified: analysis/ukb-d-30850_irnt_Liver.Rmd
Modified: analysis/ukb-d-30860_irnt_Liver.Rmd
Modified: analysis/ukb-d-30870_irnt_Liver.Rmd
Modified: analysis/ukb-d-30880_irnt_Liver.Rmd
Modified: analysis/ukb-d-30890_irnt_Liver.Rmd
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were made to the R Markdown (analysis/ukb-d-30710_irnt_Whole_Blood.Rmd) and HTML (docs/ukb-d-30710_irnt_Whole_Blood.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | cbf7408 | wesleycrouse | 2021-09-08 | adding enrichment to reports |
| html | cbf7408 | wesleycrouse | 2021-09-08 | adding enrichment to reports |
| Rmd | 4970e3e | wesleycrouse | 2021-09-08 | updating reports |
| html | 4970e3e | wesleycrouse | 2021-09-08 | updating reports |
| Rmd | dfd2b5f | wesleycrouse | 2021-09-07 | regenerating reports |
| html | dfd2b5f | wesleycrouse | 2021-09-07 | regenerating reports |
| Rmd | 61b53b3 | wesleycrouse | 2021-09-06 | updated PVE calculation |
| html | 61b53b3 | wesleycrouse | 2021-09-06 | updated PVE calculation |
| Rmd | 837dd01 | wesleycrouse | 2021-09-01 | adding additional fixedsigma report |
| Rmd | d0a5417 | wesleycrouse | 2021-08-30 | adding new reports to the index |
| Rmd | 0922de7 | wesleycrouse | 2021-08-18 | updating all reports with locus plots |
| html | 0922de7 | wesleycrouse | 2021-08-18 | updating all reports with locus plots |
| html | 1c62980 | wesleycrouse | 2021-08-11 | Updating reports |
| Rmd | 5981e80 | wesleycrouse | 2021-08-11 | Adding more reports |
| html | 5981e80 | wesleycrouse | 2021-08-11 | Adding more reports |
| Rmd | da9f015 | wesleycrouse | 2021-08-07 | adding more reports |
| html | da9f015 | wesleycrouse | 2021-08-07 | adding more reports |
Overview
These are the results of a ctwas analysis of the UK Biobank trait C-reactive protein (quantile) using Whole_Blood gene weights.
The GWAS was conducted by the Neale Lab, and the biomarker traits we analyzed are discussed here. Summary statistics were obtained from IEU OpenGWAS using GWAS ID: ukb-d-30710_irnt. Results were obtained from from IEU rather than Neale Lab because they are in a standardard format (GWAS VCF). Note that 3 of the 34 biomarker traits were not available from IEU and were excluded from analysis.
The weights are mashr GTEx v8 models on Whole_Blood eQTL obtained from PredictDB. We performed a full harmonization of the variants, including recovering strand ambiguous variants. This procedure is discussed in a separate document. (TO-DO: add report that describes harmonization)
LD matrices were computed from a 10% subset of Neale lab subjects. Subjects were matched using the plate and well information from genotyping. We included only biallelic variants with MAF>0.01 in the original Neale Lab GWAS. (TO-DO: add more details [number of subjects, variants, etc])
Weight QC
TO-DO: add enhanced QC reporting (total number of weights, why each variant was missing for all genes)
qclist_all <- list()
qc_files <- paste0(results_dir, "/", list.files(results_dir, pattern="exprqc.Rd"))
for (i in 1:length(qc_files)){
load(qc_files[i])
chr <- unlist(strsplit(rev(unlist(strsplit(qc_files[i], "_")))[1], "[.]"))[1]
qclist_all[[chr]] <- cbind(do.call(rbind, lapply(qclist,unlist)), as.numeric(substring(chr,4)))
}
qclist_all <- data.frame(do.call(rbind, qclist_all))
colnames(qclist_all)[ncol(qclist_all)] <- "chr"
rm(qclist, wgtlist, z_gene_chr)
#number of imputed weights
nrow(qclist_all)[1] 11095#number of imputed weights by chromosome
table(qclist_all$chr)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
1129 747 624 400 479 621 560 383 404 430 682 652 192 362 331
16 17 18 19 20 21 22
551 725 159 911 313 130 310 #proportion of imputed weights without missing variants
mean(qclist_all$nmiss==0)[1] 0.762776Load ctwas results
#load ctwas results
ctwas_res <- data.table::fread(paste0(results_dir, "/", analysis_id, "_ctwas.susieIrss.txt"))
#make unique identifier for regions
ctwas_res$region_tag <- paste(ctwas_res$region_tag1, ctwas_res$region_tag2, sep="_")
#compute PVE for each gene/SNP
ctwas_res$PVE = ctwas_res$susie_pip*ctwas_res$mu2/sample_size #check PVE calculation
#separate gene and SNP results
ctwas_gene_res <- ctwas_res[ctwas_res$type == "gene", ]
ctwas_gene_res <- data.frame(ctwas_gene_res)
ctwas_snp_res <- ctwas_res[ctwas_res$type == "SNP", ]
ctwas_snp_res <- data.frame(ctwas_snp_res)
#add gene information to results
sqlite <- RSQLite::dbDriver("SQLite")
db = RSQLite::dbConnect(sqlite, paste0("/project2/compbio/predictdb/mashr_models/mashr_", weight, ".db"))
query <- function(...) RSQLite::dbGetQuery(db, ...)
gene_info <- query("select gene, genename from extra")
gene_info <- query("select gene, genename, gene_type from extra")
RSQLite::dbDisconnect(db)
ctwas_gene_res <- cbind(ctwas_gene_res, gene_info[sapply(ctwas_gene_res$id, match, gene_info$gene), c("genename", "gene_type")])
#add z score to results
load(paste0(results_dir, "/", analysis_id, "_expr_z_gene.Rd"))
ctwas_gene_res$z <- z_gene[ctwas_gene_res$id,]$z
#load(paste0(results_dir, "/", analysis_id, "_expr_z_snp.Rd")) #for new version, stored after harmonization
z_snp <- readRDS(paste0(results_dir, "/", trait_id, ".RDS")) #for old version, unharmonized
z_snp <- z_snp[z_snp$id %in% ctwas_snp_res$id,] #subset snps to those included in analysis, note some are duplicated, need to match which allele was used
ctwas_snp_res$z <- z_snp$z[match(ctwas_snp_res$id, z_snp$id)] #for duplicated snps, this arbitrarily uses the first allele
ctwas_snp_res$z_flag <- as.numeric(ctwas_snp_res$id %in% z_snp$id[duplicated(z_snp$id)]) #mark the unclear z scores, flag=1
#formatting and rounding for tables
ctwas_gene_res$z <- round(ctwas_gene_res$z,2)
ctwas_snp_res$z <- round(ctwas_snp_res$z,2)
ctwas_gene_res$susie_pip <- round(ctwas_gene_res$susie_pip,3)
ctwas_snp_res$susie_pip <- round(ctwas_snp_res$susie_pip,3)
ctwas_gene_res$mu2 <- round(ctwas_gene_res$mu2,2)
ctwas_snp_res$mu2 <- round(ctwas_snp_res$mu2,2)
ctwas_gene_res$PVE <- signif(ctwas_gene_res$PVE, 2)
ctwas_snp_res$PVE <- signif(ctwas_snp_res$PVE, 2)
#merge gene and snp results with added information
ctwas_gene_res$z_flag=NA
ctwas_snp_res$genename=NA
ctwas_snp_res$gene_type=NA
ctwas_res <- rbind(ctwas_gene_res,
ctwas_snp_res[,colnames(ctwas_gene_res)])
#store columns to report
report_cols <- colnames(ctwas_gene_res)[!(colnames(ctwas_gene_res) %in% c("type", "region_tag1", "region_tag2", "cs_index", "gene_type", "z_flag", "id", "chrom", "pos"))]
first_cols <- c("genename", "region_tag")
report_cols <- c(first_cols, report_cols[!(report_cols %in% first_cols)])
report_cols_snps <- c("id", report_cols[-1])
#get number of SNPs from s1 results; adjust for thin
ctwas_res_s1 <- data.table::fread(paste0(results_dir, "/", analysis_id, "_ctwas.s1.susieIrss.txt"))
n_snps <- sum(ctwas_res_s1$type=="SNP")/thin
rm(ctwas_res_s1)Check convergence of parameters
library(ggplot2)
library(cowplot)
********************************************************Note: As of version 1.0.0, cowplot does not change the default ggplot2 theme anymore. To recover the previous behavior, execute:
theme_set(theme_cowplot())********************************************************load(paste0(results_dir, "/", analysis_id, "_ctwas.s2.susieIrssres.Rd"))
df <- data.frame(niter = rep(1:ncol(group_prior_rec), 2),
value = c(group_prior_rec[1,], group_prior_rec[2,]),
group = rep(c("Gene", "SNP"), each = ncol(group_prior_rec)))
df$group <- as.factor(df$group)
df$value[df$group=="SNP"] <- df$value[df$group=="SNP"]*thin #adjust parameter to account for thin argument
p_pi <- ggplot(df, aes(x=niter, y=value, group=group)) +
geom_line(aes(color=group)) +
geom_point(aes(color=group)) +
xlab("Iteration") + ylab(bquote(pi)) +
ggtitle("Prior mean") +
theme_cowplot()
df <- data.frame(niter = rep(1:ncol(group_prior_var_rec), 2),
value = c(group_prior_var_rec[1,], group_prior_var_rec[2,]),
group = rep(c("Gene", "SNP"), each = ncol(group_prior_var_rec)))
df$group <- as.factor(df$group)
p_sigma2 <- ggplot(df, aes(x=niter, y=value, group=group)) +
geom_line(aes(color=group)) +
geom_point(aes(color=group)) +
xlab("Iteration") + ylab(bquote(sigma^2)) +
ggtitle("Prior variance") +
theme_cowplot()
plot_grid(p_pi, p_sigma2)
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
#estimated group prior
estimated_group_prior <- group_prior_rec[,ncol(group_prior_rec)]
names(estimated_group_prior) <- c("gene", "snp")
estimated_group_prior["snp"] <- estimated_group_prior["snp"]*thin #adjust parameter to account for thin argument
print(estimated_group_prior) gene snp
0.0120339441 0.0002076959 #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
21.25937 13.74669 #report sample size
print(sample_size)[1] 343524#report group size
group_size <- c(nrow(ctwas_gene_res), n_snps)
print(group_size)[1] 11095 8697330#estimated group PVE
estimated_group_pve <- estimated_group_prior_var*estimated_group_prior*group_size/sample_size #check PVE calculation
names(estimated_group_pve) <- c("gene", "snp")
print(estimated_group_pve) gene snp
0.008262827 0.072286113 #compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))[1] 0.0316764 0.6720256Genes 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
2997 CHST10 2_58 1.000 337.73 9.8e-04 -5.33
2312 LIPA 10_57 0.994 106.00 3.1e-04 10.55
9863 LAMP1 13_62 0.990 63.32 1.8e-04 -7.99
11668 CEBPA 19_23 0.990 25.80 7.4e-05 -5.01
4610 ACP2 11_29 0.981 67.52 1.9e-04 10.30
10950 ZNF316 7_9 0.980 29.67 8.5e-05 -5.15
11023 SIPA1 11_36 0.973 29.79 8.4e-05 5.25
8924 CCDC85B 11_36 0.946 32.91 9.1e-05 5.16
11104 STARD10 11_41 0.946 55.95 1.5e-04 7.35
3323 NEK6 9_64 0.923 24.55 6.6e-05 4.73
5843 TNIP1 5_88 0.918 38.21 1.0e-04 -5.87
4731 CD164 6_73 0.912 25.58 6.8e-05 -4.83
713 SRBD1 2_28 0.899 24.03 6.3e-05 -4.52
12009 RP11-77H9.8 16_9 0.893 24.86 6.5e-05 4.68
11396 LINC01700 21_18 0.892 58.61 1.5e-04 10.43
6366 CMIP 16_46 0.868 23.73 6.0e-05 4.37
11291 LINC01806 2_97 0.861 21.62 5.4e-05 -4.22
9073 HIC1 17_3 0.860 22.65 5.7e-05 4.36
6590 NTAN1 16_15 0.858 22.57 5.6e-05 -4.35
4564 PSRC1 1_67 0.841 35.82 8.8e-05 6.13Genes 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
8464 WIPF2 17_23 0.000 5043.23 0.0e+00 -3.47
1346 CDC6 17_23 0.000 4534.82 0.0e+00 -2.73
10381 ZGPAT 20_38 0.000 4481.07 0.0e+00 4.97
1699 ARFRP1 20_38 0.700 4275.81 8.7e-03 -9.38
10436 STMN3 20_38 0.000 3386.34 0.0e+00 4.55
4733 AHI1 6_89 0.000 3098.39 0.0e+00 -2.43
1694 GMEB2 20_38 0.000 1889.62 0.0e+00 -4.08
11906 RTEL1 20_38 0.000 1306.08 0.0e+00 6.42
154 MED24 17_23 0.000 540.86 0.0e+00 4.73
6912 IL6R 1_75 0.000 522.42 1.7e-07 24.36
4687 TMEM60 7_49 0.029 420.40 3.5e-05 -3.95
4320 RARA 17_23 0.000 408.85 0.0e+00 -0.57
2997 CHST10 2_58 1.000 337.73 9.8e-04 -5.33
11612 TNFRSF6B 20_38 0.000 311.23 0.0e+00 -1.52
2953 NRBP1 2_16 0.007 303.74 6.1e-06 18.34
11036 LEPROT 1_41 0.000 300.37 3.9e-12 25.60
2956 SNX17 2_16 0.007 268.78 5.7e-06 16.87
8411 TRMT61B 2_19 0.655 235.79 4.5e-04 -4.06
11039 PPP1CB 2_19 0.021 226.71 1.4e-05 3.35
10746 LIME1 20_38 0.000 213.41 0.0e+00 -7.73Genes 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
1699 ARFRP1 20_38 0.700 4275.81 8.7e-03 -9.38
2997 CHST10 2_58 1.000 337.73 9.8e-04 -5.33
8411 TRMT61B 2_19 0.655 235.79 4.5e-04 -4.06
2312 LIPA 10_57 0.994 106.00 3.1e-04 10.55
4610 ACP2 11_29 0.981 67.52 1.9e-04 10.30
9863 LAMP1 13_62 0.990 63.32 1.8e-04 -7.99
8188 SHE 1_75 0.340 159.64 1.6e-04 -16.46
11104 STARD10 11_41 0.946 55.95 1.5e-04 7.35
11396 LINC01700 21_18 0.892 58.61 1.5e-04 10.43
9322 F2 11_28 0.796 54.65 1.3e-04 8.02
2410 MLX 17_25 0.740 61.73 1.3e-04 7.06
12263 RP11-574K11.29 10_49 0.756 49.90 1.1e-04 6.47
10038 CARD9 9_73 0.773 49.22 1.1e-04 7.50
929 IL4R 16_22 0.766 45.66 1.0e-04 -6.63
5843 TNIP1 5_88 0.918 38.21 1.0e-04 -5.87
5074 EMILIN1 2_16 0.552 59.97 9.6e-05 8.67
8924 CCDC85B 11_36 0.946 32.91 9.1e-05 5.16
208 PPP5C 19_32 0.685 45.52 9.1e-05 -6.95
4564 PSRC1 1_67 0.841 35.82 8.8e-05 6.13
10950 ZNF316 7_9 0.980 29.67 8.5e-05 -5.15Genes 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
11036 LEPROT 1_41 0.000 300.37 3.9e-12 25.60
6912 IL6R 1_75 0.000 522.42 1.7e-07 24.36
2953 NRBP1 2_16 0.007 303.74 6.1e-06 18.34
2956 SNX17 2_16 0.007 268.78 5.7e-06 16.87
8188 SHE 1_75 0.340 159.64 1.6e-04 -16.46
4677 OASL 12_74 0.000 124.91 1.7e-07 13.33
2578 MLEC 12_74 0.013 148.22 5.8e-06 12.37
10733 CENPW 6_84 0.016 113.00 5.4e-06 11.93
3596 ACADS 12_74 0.000 121.58 1.6e-07 11.68
7304 ZNF513 2_16 0.018 128.65 6.8e-06 11.23
12116 RP11-806H10.4 17_44 0.008 94.24 2.1e-06 11.21
11464 AF064858.8 21_18 0.058 66.57 1.1e-05 11.08
6645 SPPL3 12_74 0.000 52.04 6.6e-08 -10.98
4680 P2RX4 12_74 0.002 46.42 3.2e-07 -10.85
5076 ATRAID 2_16 0.007 96.26 1.9e-06 -10.82
10807 SLC44A4 6_26 0.001 139.38 3.1e-07 -10.67
2312 LIPA 10_57 0.994 106.00 3.1e-04 10.55
11496 AF064858.11 21_18 0.049 59.60 8.6e-06 10.54
9277 FCER1A 1_79 0.000 170.07 0.0e+00 -10.43
11396 LINC01700 21_18 0.892 58.61 1.5e-04 10.43Comparing 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.02361424#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
11036 LEPROT 1_41 0.000 300.37 3.9e-12 25.60
6912 IL6R 1_75 0.000 522.42 1.7e-07 24.36
2953 NRBP1 2_16 0.007 303.74 6.1e-06 18.34
2956 SNX17 2_16 0.007 268.78 5.7e-06 16.87
8188 SHE 1_75 0.340 159.64 1.6e-04 -16.46
4677 OASL 12_74 0.000 124.91 1.7e-07 13.33
2578 MLEC 12_74 0.013 148.22 5.8e-06 12.37
10733 CENPW 6_84 0.016 113.00 5.4e-06 11.93
3596 ACADS 12_74 0.000 121.58 1.6e-07 11.68
7304 ZNF513 2_16 0.018 128.65 6.8e-06 11.23
12116 RP11-806H10.4 17_44 0.008 94.24 2.1e-06 11.21
11464 AF064858.8 21_18 0.058 66.57 1.1e-05 11.08
6645 SPPL3 12_74 0.000 52.04 6.6e-08 -10.98
4680 P2RX4 12_74 0.002 46.42 3.2e-07 -10.85
5076 ATRAID 2_16 0.007 96.26 1.9e-06 -10.82
10807 SLC44A4 6_26 0.001 139.38 3.1e-07 -10.67
2312 LIPA 10_57 0.994 106.00 3.1e-04 10.55
11496 AF064858.11 21_18 0.049 59.60 8.6e-06 10.54
9277 FCER1A 1_79 0.000 170.07 0.0e+00 -10.43
11396 LINC01700 21_18 0.892 58.61 1.5e-04 10.43Locus 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: 1_41"
genename region_tag susie_pip mu2 PVE z
7060 RAVER2 1_41 0.001 187.91 3.4e-07 -7.05
12303 RP4-535B20.4 1_41 0.000 32.24 3.6e-14 1.30
7059 JAK1 1_41 0.000 13.83 1.8e-15 -2.42
7058 AK4 1_41 0.000 12.10 8.8e-16 -1.12
3108 LEPR 1_41 0.000 41.27 7.8e-15 -8.87
11036 LEPROT 1_41 0.000 300.37 3.9e-12 25.60
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 1_75"
genename region_tag susie_pip mu2 PVE z
10736 LOR 1_75 0.000 6.51 2.3e-09 0.26
5636 S100A8 1_75 0.000 8.18 3.5e-09 -0.01
7190 S100A9 1_75 0.000 10.86 6.1e-09 -0.40
7191 S100A12 1_75 0.000 10.86 6.1e-09 -0.40
10519 S100A6 1_75 0.000 15.44 6.6e-09 2.70
10267 S100A5 1_75 0.000 16.95 1.2e-08 -2.36
10060 S100A3 1_75 0.000 5.11 1.5e-09 0.66
10217 S100A4 1_75 0.000 5.82 1.8e-09 -0.75
10193 S100A14 1_75 0.000 19.61 1.5e-08 -3.99
6906 CHTOP 1_75 0.000 31.02 4.2e-08 -5.12
10176 S100A13 1_75 0.000 10.11 3.5e-09 1.33
6905 S100A1 1_75 0.000 14.67 6.8e-09 3.55
5638 SNAPIN 1_75 0.000 14.70 1.0e-08 -2.77
8201 NPR1 1_75 0.000 18.20 1.1e-08 4.22
5639 SLC27A3 1_75 0.000 12.15 6.0e-09 -1.03
5646 GATAD2B 1_75 0.000 6.53 1.9e-09 -2.14
10678 DENND4B 1_75 0.003 73.51 6.0e-07 -7.56
5643 CREB3L4 1_75 0.000 13.72 6.4e-09 -3.77
5641 SLC39A1 1_75 0.000 5.86 1.7e-09 -1.77
5635 RAB13 1_75 0.000 27.43 2.3e-08 3.60
5637 TPM3 1_75 0.000 11.30 3.6e-09 3.01
5640 UBAP2L 1_75 0.000 12.60 4.4e-09 1.83
5642 HAX1 1_75 0.000 11.58 7.0e-09 -1.57
5645 AQP10 1_75 0.000 44.68 2.1e-08 -5.82
5631 ATP8B2 1_75 0.000 27.13 2.5e-08 0.59
6912 IL6R 1_75 0.000 522.42 1.7e-07 24.36
8188 SHE 1_75 0.340 159.64 1.6e-04 -16.46
6913 UBE2Q1 1_75 0.000 30.75 1.0e-08 -8.91
6914 CHRNB2 1_75 0.000 21.67 1.0e-08 0.19
6911 ADAR 1_75 0.023 40.08 2.7e-06 6.77
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 2_16"
genename region_tag susie_pip mu2 PVE z
11045 SLC35F6 2_16 0.041 42.26 5.0e-06 -6.13
3366 TMEM214 2_16 0.015 10.23 4.3e-07 2.73
5074 EMILIN1 2_16 0.552 59.97 9.6e-05 8.67
5061 KHK 2_16 0.010 8.02 2.4e-07 -2.52
5059 CGREF1 2_16 0.008 6.84 1.7e-07 -0.98
5070 PREB 2_16 0.011 16.78 5.4e-07 -2.42
5076 ATRAID 2_16 0.007 96.26 1.9e-06 -10.82
1090 CAD 2_16 0.117 32.11 1.1e-05 -4.33
5071 SLC5A6 2_16 0.054 28.56 4.5e-06 3.29
7303 UCN 2_16 0.024 36.43 2.6e-06 -6.54
2952 GTF3C2 2_16 0.015 30.42 1.4e-06 6.21
2956 SNX17 2_16 0.007 268.78 5.7e-06 16.87
7304 ZNF513 2_16 0.018 128.65 6.8e-06 11.23
2953 NRBP1 2_16 0.007 303.74 6.1e-06 18.34
5057 IFT172 2_16 0.072 57.81 1.2e-05 -7.85
1087 GCKR 2_16 0.056 54.88 8.9e-06 7.75
10613 GPN1 2_16 0.007 37.22 7.5e-07 5.31
9018 CCDC121 2_16 0.007 9.64 2.0e-07 -2.46
6660 BRE 2_16 0.155 55.13 2.5e-05 6.90
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 12_74"
genename region_tag susie_pip mu2 PVE z
2658 PRKAB1 12_74 0.000 5.43 6.9e-09 -0.88
4681 BICDL1 12_74 0.000 9.85 1.4e-08 -2.58
2662 RAB35 12_74 0.001 7.25 1.3e-08 0.16
1220 GCN1 12_74 0.003 51.06 4.3e-07 8.09
1221 RPLP0 12_74 0.001 8.21 1.3e-08 -3.46
1222 PXN 12_74 0.001 46.65 8.5e-08 8.83
1223 SIRT4 12_74 0.004 65.45 7.8e-07 -8.84
2664 COX6A1 12_74 0.001 20.05 3.5e-08 -0.99
11883 GATC 12_74 0.000 12.33 1.6e-08 -4.21
8384 TRIAP1 12_74 0.001 15.48 4.0e-08 -0.95
1205 DYNLL1 12_74 0.001 7.28 1.6e-08 1.63
2572 COQ5 12_74 0.001 7.44 1.5e-08 2.50
278 RNF10 12_74 0.001 9.91 1.4e-08 3.33
7876 POP5 12_74 0.000 6.56 9.3e-09 -2.55
6639 CABP1 12_74 0.001 30.06 6.2e-08 5.55
2578 MLEC 12_74 0.013 148.22 5.8e-06 12.37
8953 UNC119B 12_74 0.001 30.10 6.5e-08 3.77
3596 ACADS 12_74 0.000 121.58 1.6e-07 11.68
6645 SPPL3 12_74 0.000 52.04 6.6e-08 -10.98
6651 C12orf43 12_74 0.003 30.95 2.5e-07 -0.60
4677 OASL 12_74 0.000 124.91 1.7e-07 13.33
1211 P2RX7 12_74 0.000 11.97 1.6e-08 -1.91
4680 P2RX4 12_74 0.002 46.42 3.2e-07 -10.85
12411 RP11-340F14.6 12_74 0.001 12.39 4.3e-08 -4.37
2581 CAMKK2 12_74 0.001 14.12 4.5e-08 -1.51
1214 ANAPC5 12_74 0.001 12.36 2.1e-08 -2.96
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 6_84"
genename region_tag susie_pip mu2 PVE z
2687 HDDC2 6_84 0.021 14.31 8.9e-07 1.65
2689 NCOA7 6_84 0.008 5.41 1.3e-07 1.17
2688 HINT3 6_84 0.008 6.26 1.4e-07 -0.84
10733 CENPW 6_84 0.016 113.00 5.4e-06 11.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
6872 rs1772722 1_15 1.000 38.64 1.1e-04 -5.08
8327 rs79598313 1_18 1.000 107.17 3.1e-04 -11.02
18254 rs11811946 1_41 1.000 514.78 1.5e-03 32.07
35862 rs12036827 1_79 1.000 319.09 9.3e-04 6.61
35900 rs7551731 1_79 1.000 3251.70 9.5e-03 -68.36
51304 rs78270913 1_106 1.000 797.15 2.3e-03 -2.05
51322 rs61820396 1_106 1.000 759.37 2.2e-03 -2.65
59770 rs4006577 1_122 1.000 34.64 1.0e-04 -5.80
64048 rs12239046 1_131 1.000 172.03 5.0e-04 16.09
73774 rs780093 2_16 1.000 743.09 2.2e-03 -29.34
74075 rs569546056 2_19 1.000 369.92 1.1e-03 -2.15
88507 rs6737957 2_47 1.000 538.00 1.6e-03 -4.51
88515 rs550064763 2_47 1.000 776.80 2.3e-03 -1.38
113041 rs1862069 2_102 1.000 70.25 2.0e-04 8.97
131482 rs142215640 2_136 1.000 216.81 6.3e-04 3.76
313657 rs10866657 5_103 1.000 39.19 1.1e-04 -6.89
361993 rs199804242 6_89 1.000 13408.77 3.9e-02 4.09
362488 rs62432712 6_92 1.000 48.80 1.4e-04 7.23
382948 rs11983782 7_20 1.000 78.36 2.3e-04 -8.97
400030 rs761767938 7_49 1.000 839.00 2.4e-03 1.70
400038 rs1544459 7_49 1.000 819.87 2.4e-03 2.10
429615 rs1703982 8_11 1.000 309.99 9.0e-04 -7.50
429641 rs758184196 8_11 1.000 346.66 1.0e-03 -2.22
429895 rs7012814 8_12 1.000 206.94 6.0e-04 19.71
429906 rs13265179 8_12 1.000 179.19 5.2e-04 -14.17
556604 rs34623292 11_10 1.000 62.54 1.8e-04 8.73
621854 rs55692966 12_56 1.000 94.00 2.7e-04 -10.80
629775 rs1169286 12_74 1.000 887.88 2.6e-03 -46.10
690945 rs73349296 14_50 1.000 45.86 1.3e-04 -7.23
703134 rs72743115 15_22 1.000 54.84 1.6e-04 -7.51
705743 rs340029 15_27 1.000 138.89 4.0e-04 12.13
729620 rs11641493 16_27 1.000 84.67 2.5e-04 -6.81
729621 rs112848702 16_27 1.000 100.19 2.9e-04 -8.99
729645 rs62039688 16_27 1.000 259.03 7.5e-04 -15.21
729727 rs17616063 16_27 1.000 764.62 2.2e-03 -27.92
752242 rs3110630 17_22 1.000 47.15 1.4e-04 -2.60
752243 rs12938438 17_22 1.000 85.54 2.5e-04 -4.93
752244 rs2189303 17_22 1.000 57.83 1.7e-04 2.36
758665 rs1801689 17_38 1.000 68.63 2.0e-04 8.52
780406 rs1217565 18_30 1.000 57.31 1.7e-04 -8.10
800797 rs140965448 19_26 1.000 34.79 1.0e-04 -5.88
802002 rs116881820 19_31 1.000 545.02 1.6e-03 -30.35
802007 rs405509 19_31 1.000 2153.42 6.3e-03 19.05
802008 rs440446 19_31 1.000 2257.15 6.6e-03 -28.59
802011 rs814573 19_31 1.000 2316.31 6.7e-03 -73.69
802865 rs601338 19_33 1.000 97.49 2.8e-04 10.26
861702 rs12083537 1_75 1.000 295.13 8.6e-04 -22.86
884509 rs777700134 2_58 1.000 1653.26 4.8e-03 -2.14
949427 rs2266008 10_57 1.000 73.65 2.1e-04 8.84
1016038 rs138395719 17_23 1.000 5646.52 1.6e-02 3.84
1021409 rs11658269 17_44 1.000 142.72 4.2e-04 10.47
1045521 rs202143810 20_38 1.000 5945.28 1.7e-02 -5.38
32016 rs2618039 1_69 0.999 35.43 1.0e-04 5.92
333004 rs9471632 6_32 0.999 35.09 1.0e-04 5.93
596364 rs12824533 12_11 0.999 31.44 9.1e-05 -5.50
629790 rs2258043 12_74 0.999 683.28 2.0e-03 -40.50
683118 rs2074648 14_34 0.999 43.53 1.3e-04 6.02
703128 rs8040040 15_22 0.999 61.76 1.8e-04 -6.81
18285 rs35228056 1_41 0.998 130.19 3.8e-04 -3.20
330591 rs41258084 6_27 0.998 35.11 1.0e-04 6.46
710223 rs62009152 15_36 0.998 30.68 8.9e-05 -5.07
716272 rs6496132 15_46 0.998 29.92 8.7e-05 5.27
7630 rs10794644 1_17 0.997 35.78 1.0e-04 -6.01
18612 rs4655616 1_42 0.997 45.38 1.3e-04 6.30
35898 rs77289344 1_79 0.997 710.94 2.1e-03 42.02
52752 rs1223802 1_108 0.997 34.80 1.0e-04 -6.06
430181 rs506276 8_13 0.995 56.53 1.6e-04 9.82
629778 rs2393775 12_74 0.995 1829.69 5.3e-03 59.51
683160 rs61986270 14_34 0.993 32.29 9.3e-05 4.81
51346 rs113267426 1_106 0.992 46.96 1.4e-04 -0.88
135128 rs12619647 2_144 0.992 36.97 1.1e-04 6.21
744165 rs71368518 17_5 0.992 28.46 8.2e-05 -4.87
414627 rs3757387 7_79 0.990 33.31 9.6e-05 5.71
54984 rs7523735 1_112 0.988 58.36 1.7e-04 -7.82
6886 rs72660319 1_15 0.987 38.85 1.1e-04 5.22
383054 rs11770879 7_20 0.987 37.11 1.1e-04 -6.65
319351 rs10458103 6_7 0.985 39.08 1.1e-04 6.37
84433 rs4672266 2_39 0.983 26.90 7.7e-05 -4.67
506976 rs115478735 9_70 0.983 112.59 3.2e-04 10.52
609581 rs34071274 12_34 0.983 27.63 7.9e-05 4.98
64036 rs12048215 1_131 0.982 42.47 1.2e-04 9.45
624216 rs4764939 12_62 0.981 53.56 1.5e-04 -7.45
64030 rs4508048 1_131 0.979 36.67 1.0e-04 -7.00
795825 rs3794991 19_15 0.979 68.47 2.0e-04 8.53
800126 rs1688031 19_24 0.979 45.78 1.3e-04 -6.75
825655 rs8121794 20_37 0.979 30.98 8.8e-05 5.53
18309 rs200241668 1_41 0.978 1630.29 4.6e-03 -50.44
465746 rs6987702 8_83 0.977 103.28 2.9e-04 6.94
667645 rs2332328 14_3 0.977 32.32 9.2e-05 -5.56
289098 rs56821385 5_55 0.976 26.08 7.4e-05 4.92
534220 rs1269867 10_50 0.976 28.68 8.2e-05 -5.25
175751 rs28433979 3_81 0.974 26.68 7.6e-05 5.09
92677 rs4441469 2_54 0.968 28.17 7.9e-05 5.15
556034 rs10831676 11_9 0.968 38.02 1.1e-04 -6.14
600404 rs11047225 12_18 0.968 40.18 1.1e-04 6.41
662232 rs76620584 13_52 0.968 26.09 7.3e-05 -4.94
299777 rs6595550 5_76 0.967 27.03 7.6e-05 -5.06
73799 rs9679004 2_16 0.966 28.25 7.9e-05 2.54
97149 rs72831808 2_67 0.966 30.85 8.7e-05 -5.09
1021434 rs41522945 17_44 0.966 54.12 1.5e-04 -7.20
1021471 rs4969172 17_44 0.966 95.84 2.7e-04 -10.85
11370 rs3103771 1_25 0.965 27.55 7.7e-05 -4.84
719954 rs2601781 16_4 0.964 30.82 8.7e-05 5.41
277507 rs67715745 5_31 0.962 26.70 7.5e-05 4.97
825917 rs62206958 20_37 0.960 25.51 7.1e-05 4.91
779933 rs377181240 18_30 0.959 28.09 7.8e-05 5.17
506371 rs2966361 9_69 0.955 26.14 7.3e-05 4.96
821228 rs6020459 20_30 0.954 47.10 1.3e-04 6.99
603308 rs60798338 12_21 0.951 57.50 1.6e-04 -5.84
361992 rs2327654 6_89 0.948 13361.19 3.7e-02 4.19
820444 rs6066487 20_29 0.948 25.90 7.1e-05 -4.30
55268 rs12726661 1_113 0.944 35.64 9.8e-05 -6.07
64049 rs114193355 1_131 0.938 39.76 1.1e-04 -9.08
325887 rs75080831 6_19 0.937 31.56 8.6e-05 5.55
12390 rs943513 1_27 0.936 32.41 8.8e-05 -6.03
382841 rs6973240 7_20 0.936 37.30 1.0e-04 4.60
97205 rs17628215 2_67 0.934 29.48 8.0e-05 5.03
406417 rs10257273 7_61 0.931 99.57 2.7e-04 -8.05
330396 rs6916779 6_27 0.930 75.53 2.0e-04 6.96
684457 rs2270420 14_36 0.930 24.41 6.6e-05 4.49
351778 rs9496567 6_67 0.929 24.92 6.7e-05 4.73
465748 rs4604455 8_83 0.929 24.95 6.7e-05 -2.33
826007 rs6122476 20_37 0.926 24.01 6.5e-05 4.83
6834 rs192212396 1_15 0.925 25.33 6.8e-05 -4.52
499303 rs2777804 9_52 0.925 32.43 8.7e-05 5.56
226940 rs149027545 4_59 0.924 31.83 8.6e-05 -5.59
282730 rs200933157 5_42 0.924 29.14 7.8e-05 4.96
458986 rs35521407 8_69 0.923 41.71 1.1e-04 6.61
721998 rs6501005 16_8 0.921 24.33 6.5e-05 4.63
789046 rs351988 19_2 0.916 27.68 7.4e-05 -5.07
213343 rs12498157 4_35 0.914 25.15 6.7e-05 4.52
471409 rs1078141 8_92 0.912 30.62 8.1e-05 5.41
600403 rs11047224 12_18 0.911 54.76 1.5e-04 -8.45
166353 rs138332280 3_64 0.910 26.11 6.9e-05 5.31
125715 rs1441171 2_126 0.909 83.05 2.2e-04 -9.41
743596 rs35316912 17_2 0.908 31.75 8.4e-05 5.84
566366 rs35243581 11_27 0.903 42.59 1.1e-04 6.51
139485 rs2279440 3_7 0.895 39.42 1.0e-04 6.48
219045 rs116774130 4_45 0.894 24.33 6.3e-05 -4.41
18596 rs2166311 1_42 0.888 25.25 6.5e-05 -3.44
462813 rs2737251 8_78 0.884 67.78 1.7e-04 -9.54
621100 rs35000205 12_55 0.879 27.04 6.9e-05 4.92
144564 rs2362186 3_18 0.878 25.40 6.5e-05 -4.58
790087 rs111334206 19_4 0.875 26.38 6.7e-05 -4.71
800348 rs7254601 19_24 0.874 26.69 6.8e-05 -4.80
845303 rs112169312 22_17 0.874 39.89 1.0e-04 -7.05
73971 rs7606480 2_19 0.873 29.49 7.5e-05 -5.66
97137 rs45487895 2_67 0.872 25.49 6.5e-05 -3.12
607392 rs2429465 12_29 0.872 25.45 6.5e-05 -4.57
792953 rs113439314 19_9 0.869 36.87 9.3e-05 5.91
288083 rs7444298 5_52 0.867 28.68 7.2e-05 -5.06
289395 rs13170906 5_56 0.866 24.37 6.1e-05 -4.52
201045 rs150164330 4_11 0.858 23.89 6.0e-05 -4.55
801995 rs2972559 19_31 0.855 220.71 5.5e-04 -24.15
111525 rs9287838 2_99 0.853 24.92 6.2e-05 -4.60
742234 rs904787 16_54 0.852 55.40 1.4e-04 7.51
333537 rs3024991 6_33 0.846 25.68 6.3e-05 4.83
430972 rs3824116 8_15 0.845 41.68 1.0e-04 -8.30
197954 rs3752442 4_4 0.840 29.25 7.2e-05 5.35
626770 rs75622376 12_67 0.839 23.72 5.8e-05 -4.36
585516 rs68170644 11_69 0.837 27.71 6.8e-05 -4.69
206386 rs34811474 4_21 0.835 24.77 6.0e-05 -4.61
844720 rs4821799 22_16 0.835 155.25 3.8e-04 -12.41
1021395 rs57933884 17_44 0.830 47.81 1.2e-04 -2.96
277481 rs1499279 5_31 0.828 67.53 1.6e-04 -8.45
348575 rs854863 6_61 0.828 25.87 6.2e-05 -4.88
720850 rs1436387 16_6 0.828 26.18 6.3e-05 4.80
630924 rs4765621 12_76 0.825 26.26 6.3e-05 -4.88
532632 rs117206424 10_48 0.822 30.49 7.3e-05 -5.39
144944 rs11711833 3_18 0.821 27.76 6.6e-05 -4.95
382898 rs7782803 7_20 0.820 76.91 1.8e-04 -8.01
382868 rs35345753 7_20 0.808 47.99 1.1e-04 -5.81SNPs 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
361993 rs199804242 6_89 1.000 13408.77 3.9e-02 4.09
361992 rs2327654 6_89 0.948 13361.19 3.7e-02 4.19
362009 rs6923513 6_89 0.754 13360.46 2.9e-02 4.17
361996 rs113527452 6_89 0.000 13289.43 1.0e-05 4.07
362001 rs200796875 6_89 0.000 13214.09 3.1e-07 4.13
362014 rs7756915 6_89 0.000 13126.37 3.7e-12 4.04
362007 rs6570040 6_89 0.000 12608.54 0.0e+00 4.36
361994 rs6570031 6_89 0.000 12580.08 0.0e+00 4.40
361995 rs9389323 6_89 0.000 12574.10 0.0e+00 4.34
362011 rs9321531 6_89 0.000 11020.37 0.0e+00 3.54
361984 rs9321528 6_89 0.000 10890.02 0.0e+00 3.97
362012 rs9494389 6_89 0.000 10352.99 0.0e+00 3.56
362016 rs5880262 6_89 0.000 10327.55 0.0e+00 3.39
361990 rs2208574 6_89 0.000 9998.71 0.0e+00 3.79
361989 rs1033755 6_89 0.000 9991.56 0.0e+00 3.56
361987 rs2038551 6_89 0.000 9826.09 0.0e+00 4.35
361998 rs9494377 6_89 0.000 9811.91 0.0e+00 3.60
361985 rs2038550 6_89 0.000 9799.39 0.0e+00 4.31
361974 rs6570026 6_89 0.000 8117.34 0.0e+00 3.38
361970 rs6926161 6_89 0.000 8014.49 0.0e+00 3.40
361979 rs6930773 6_89 0.000 7881.99 0.0e+00 3.94
361966 rs6925959 6_89 0.000 6775.54 0.0e+00 2.94
361971 rs6917005 6_89 0.000 6551.03 0.0e+00 3.09
1045521 rs202143810 20_38 1.000 5945.28 1.7e-02 -5.38
1045518 rs6089961 20_38 0.460 5886.32 7.9e-03 -5.48
1045520 rs2738758 20_38 0.460 5886.32 7.9e-03 -5.48
1045501 rs2750483 20_38 0.320 5884.64 5.5e-03 -5.48
1045499 rs35201382 20_38 0.140 5884.11 2.4e-03 -5.46
1045500 rs67468102 20_38 0.234 5883.77 4.0e-03 -5.48
1045496 rs2315009 20_38 0.171 5882.69 2.9e-03 -5.47
1045505 rs2259858 20_38 0.006 5859.16 1.1e-04 -5.47
1045498 rs35046559 20_38 0.002 5855.24 3.8e-05 -5.46
1045497 rs547768273 20_38 0.006 5830.89 1.1e-04 -5.54
1045517 rs145835311 20_38 0.000 5829.38 3.6e-06 -5.36
1045533 rs201555224 20_38 0.000 5658.55 2.7e-16 5.50
1045548 rs34251386 20_38 0.000 5656.17 8.4e-17 -5.48
1045523 rs2263805 20_38 0.000 5655.55 4.0e-17 -5.45
1045553 rs2427529 20_38 0.000 5655.53 9.1e-17 -5.48
1045544 rs1758205 20_38 0.000 5654.74 4.4e-17 -5.46
1045545 rs2263806 20_38 0.000 5654.28 2.1e-16 -5.47
1045493 rs2750482 20_38 0.000 5652.14 1.8e-17 -5.46
1045492 rs2750481 20_38 0.000 5652.12 1.8e-17 -5.46
1045490 rs2872881 20_38 0.000 5652.09 1.8e-17 -5.46
1045488 rs2750480 20_38 0.000 5650.28 1.1e-17 -5.45
1045579 rs2427534 20_38 0.000 5650.10 1.8e-17 -5.45
1045577 rs1151622 20_38 0.000 5649.69 3.7e-17 -5.47
1045552 rs6122154 20_38 0.000 5649.38 5.7e-17 -5.49
1045566 rs71325464 20_38 0.000 5649.12 1.5e-17 -5.45
1045567 rs67616625 20_38 0.000 5648.62 1.3e-17 -5.44
1045557 rs2427530 20_38 0.000 5648.56 1.8e-17 -5.45SNPs 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
361993 rs199804242 6_89 1.000 13408.77 0.0390 4.09
361992 rs2327654 6_89 0.948 13361.19 0.0370 4.19
362009 rs6923513 6_89 0.754 13360.46 0.0290 4.17
1045521 rs202143810 20_38 1.000 5945.28 0.0170 -5.38
1016038 rs138395719 17_23 1.000 5646.52 0.0160 3.84
35900 rs7551731 1_79 1.000 3251.70 0.0095 -68.36
1045518 rs6089961 20_38 0.460 5886.32 0.0079 -5.48
1045520 rs2738758 20_38 0.460 5886.32 0.0079 -5.48
802011 rs814573 19_31 1.000 2316.31 0.0067 -73.69
802008 rs440446 19_31 1.000 2257.15 0.0066 -28.59
802007 rs405509 19_31 1.000 2153.42 0.0063 19.05
1045501 rs2750483 20_38 0.320 5884.64 0.0055 -5.48
629778 rs2393775 12_74 0.995 1829.69 0.0053 59.51
884509 rs777700134 2_58 1.000 1653.26 0.0048 -2.14
1016037 rs9893520 17_23 0.296 5597.09 0.0048 3.88
18309 rs200241668 1_41 0.978 1630.29 0.0046 -50.44
35893 rs2808628 1_79 0.474 3225.29 0.0044 -68.11
1045500 rs67468102 20_38 0.234 5883.77 0.0040 -5.48
1016042 rs9897092 17_23 0.219 5596.82 0.0036 3.87
1016073 rs9914531 17_23 0.181 5587.04 0.0029 3.93
1045496 rs2315009 20_38 0.171 5882.69 0.0029 -5.47
629775 rs1169286 12_74 1.000 887.88 0.0026 -46.10
1015978 rs72836641 17_23 0.156 5592.05 0.0025 3.88
400030 rs761767938 7_49 1.000 839.00 0.0024 1.70
400038 rs1544459 7_49 1.000 819.87 0.0024 2.10
1016006 rs16965687 17_23 0.147 5594.64 0.0024 3.87
1016007 rs9916782 17_23 0.150 5594.71 0.0024 3.87
1016008 rs79880154 17_23 0.148 5596.08 0.0024 3.86
1045499 rs35201382 20_38 0.140 5884.11 0.0024 -5.46
35895 rs2794521 1_79 0.533 1469.86 0.0023 3.73
51304 rs78270913 1_106 1.000 797.15 0.0023 -2.05
88515 rs550064763 2_47 1.000 776.80 0.0023 -1.38
1015993 rs7359537 17_23 0.138 5595.80 0.0023 3.86
51322 rs61820396 1_106 1.000 759.37 0.0022 -2.65
73774 rs780093 2_16 1.000 743.09 0.0022 -29.34
729727 rs17616063 16_27 1.000 764.62 0.0022 -27.92
1015986 rs72836642 17_23 0.134 5594.83 0.0022 3.86
1015996 rs9892404 17_23 0.138 5595.79 0.0022 3.86
1016000 rs9906409 17_23 0.136 5595.78 0.0022 3.85
1016079 rs73983025 17_23 0.137 5588.21 0.0022 3.91
35898 rs77289344 1_79 0.997 710.94 0.0021 42.02
35901 rs3116652 1_79 0.467 1470.27 0.0020 3.70
629790 rs2258043 12_74 0.999 683.28 0.0020 -40.50
1016098 rs72836660 17_23 0.121 5586.57 0.0020 3.92
1016091 rs58822685 17_23 0.117 5587.87 0.0019 3.91
1016094 rs9916460 17_23 0.115 5587.82 0.0019 3.91
88527 rs630241 2_47 0.758 765.24 0.0017 -2.19
1016087 rs72836656 17_23 0.105 5587.14 0.0017 3.91
88507 rs6737957 2_47 1.000 538.00 0.0016 -4.51
802002 rs116881820 19_31 1.000 545.02 0.0016 -30.35SNPs 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
802011 rs814573 19_31 1.000 2316.31 6.7e-03 -73.69
35900 rs7551731 1_79 1.000 3251.70 9.5e-03 -68.36
35893 rs2808628 1_79 0.474 3225.29 4.4e-03 -68.11
629778 rs2393775 12_74 0.995 1829.69 5.3e-03 59.51
629777 rs7979478 12_74 0.035 1833.87 1.9e-04 59.29
629788 rs1169311 12_74 0.001 1371.92 3.5e-06 -53.81
629765 rs2701194 12_74 0.037 1249.48 1.3e-04 50.78
629784 rs1169300 12_74 0.000 1120.78 1.1e-06 -50.74
18309 rs200241668 1_41 0.978 1630.29 4.6e-03 -50.44
18319 rs13375019 1_41 0.017 1622.76 8.1e-05 -50.37
18324 rs10789193 1_41 0.005 1619.85 2.2e-05 -50.34
18305 rs11208691 1_41 0.000 1606.76 4.8e-08 -50.15
35896 rs3116635 1_79 0.000 1044.27 0.0e+00 49.99
35902 rs12727021 1_79 0.000 1040.68 0.0e+00 49.98
18301 rs4655743 1_41 0.000 1559.65 5.5e-14 -49.67
18310 rs57274629 1_41 0.000 1605.20 7.9e-10 -49.50
18335 rs4655582 1_41 0.000 1517.13 6.6e-14 -49.01
18312 rs6697515 1_41 0.000 1503.56 9.6e-14 -48.34
18346 rs7515766 1_41 0.000 1438.18 4.3e-14 48.22
18348 rs7541434 1_41 0.000 1433.53 4.1e-14 48.18
629792 rs2258287 12_74 0.000 927.11 3.8e-07 -47.17
629775 rs1169286 12_74 1.000 887.88 2.6e-03 -46.10
35898 rs77289344 1_79 0.997 710.94 2.1e-03 42.02
35894 rs16842559 1_79 0.002 704.99 4.5e-06 41.86
35891 rs79162334 1_79 0.001 704.32 1.5e-06 41.69
35889 rs16842502 1_79 0.000 697.77 2.0e-07 41.60
629790 rs2258043 12_74 0.999 683.28 2.0e-03 -40.50
35904 rs12760041 1_79 0.000 454.78 0.0e+00 39.55
35903 rs12077166 1_79 0.000 605.43 0.0e+00 38.88
861830 rs12133641 1_75 0.479 1097.75 1.5e-03 -38.83
861824 rs12126142 1_75 0.266 1096.13 8.5e-04 -38.80
861808 rs61812598 1_75 0.151 1095.10 4.8e-04 -38.78
861827 rs4129267 1_75 0.084 1093.20 2.7e-04 -38.77
861806 rs12730935 1_75 0.018 1092.51 5.8e-05 -38.74
861829 rs2228145 1_75 0.009 1088.01 2.8e-05 -38.69
35905 rs11265266 1_79 0.000 482.22 0.0e+00 38.50
629757 rs7969196 12_74 0.001 697.29 1.8e-06 -38.32
861810 rs7512646 1_75 0.002 1063.10 5.0e-06 -38.25
861811 rs7529229 1_75 0.002 1063.12 5.0e-06 -38.25
861814 rs12404936 1_75 0.001 1058.26 4.2e-06 -38.18
629801 rs2859398 12_74 0.000 651.35 6.9e-07 -38.17
861813 rs12403159 1_75 0.001 1058.19 4.2e-06 -38.17
861800 rs4576655 1_75 0.001 1050.31 3.5e-06 -38.05
861801 rs4537545 1_75 0.001 1043.68 2.9e-06 -37.97
861799 rs11265613 1_75 0.001 1026.31 2.8e-06 -37.68
18331 rs55737039 1_41 0.000 609.01 4.1e-14 -37.67
18330 rs34618644 1_41 0.000 606.73 4.5e-14 -37.65
18339 rs6699885 1_41 0.000 603.95 4.6e-14 37.58
861790 rs12753254 1_75 0.001 1015.56 3.4e-06 -37.48
861791 rs12730036 1_75 0.001 1015.29 3.3e-06 -37.48Gene 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] 20if (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 negative regulation of cellular process (GO:0048523) 5/566
2 negative regulation of cell growth (GO:0030308) 3/126
3 negative regulation of growth (GO:0045926) 3/126
Adjusted.P.value Genes
1 0.0144967 CD164;CEBPA;PSRC1;CCDC85B;SIPA1
2 0.0144967 PSRC1;CCDC85B;SIPA1
3 0.0144967 PSRC1;CCDC85B;SIPA1
[1] "GO_Cellular_Component_2021"
Term Overlap Adjusted.P.value
1 lytic vacuole (GO:0000323) 4/219 0.001889577
2 lysosome (GO:0005764) 4/477 0.018282578
Genes
1 CD164;LAMP1;ACP2;LIPA
2 CD164;LAMP1;ACP2;LIPA
[1] "GO_Molecular_Function_2021"
[1] Term Overlap Adjusted.P.value Genes
<0 rows> (or 0-length row.names)LINC01806 gene(s) from the input list not found in DisGeNET CURATEDLAMP1 gene(s) from the input list not found in DisGeNET CURATEDCHST10 gene(s) from the input list not found in DisGeNET CURATEDPSRC1 gene(s) from the input list not found in DisGeNET CURATEDZNF316 gene(s) from the input list not found in DisGeNET CURATEDRP11-77H9.8 gene(s) from the input list not found in DisGeNET CURATEDSRBD1 gene(s) from the input list not found in DisGeNET CURATEDNTAN1 gene(s) from the input list not found in DisGeNET CURATEDCCDC85B gene(s) from the input list not found in DisGeNET CURATEDLINC01700 gene(s) from the input list not found in DisGeNET CURATED Description FDR Ratio
8 Cholesterol Ester Storage Disease 0.009161914 1/10
30 polyps 0.009161914 1/10
38 Wolman Disease 0.009161914 1/10
41 Adult Acute Myeloblastic Leukemia 0.009161914 1/10
42 Childhood Acute Myeloid Leukemia 0.009161914 1/10
48 Acid Phosphatase Deficiency 0.009161914 1/10
70 Acid cholesteryl ester hydrolase deficiency, type 2 0.009161914 1/10
75 DEAFNESS, AUTOSOMAL DOMINANT 66 0.009161914 1/10
78 Acute myeloid leukemia with CEBPA somatic mutations 0.009161914 1/10
27 Neoplasms, Radiation-Induced 0.011774139 1/10
BgRatio
8 1/9703
30 1/9703
38 1/9703
41 1/9703
42 1/9703
48 1/9703
70 1/9703
75 1/9703
78 1/9703
27 2/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