Apoliprotein B (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
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-30640_irnt_Liver.Rmd) and HTML (docs/ukb-d-30640_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 Apoliprotein B (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-30640_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.0121318812 0.0001595173 #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
33.18673 11.66628 #report sample size
print(sample_size)[1] 342590#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.01281104 0.04724453 #compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))[1] 0.04478502 0.44057500Genes 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
4435 PSRC1 1_67 1.000 2883.70 8.4e-03 -55.48
6391 TTC39B 9_13 1.000 55.22 1.6e-04 -7.30
5991 FADS1 11_34 1.000 303.61 8.9e-04 17.90
12008 HPR 16_38 1.000 234.79 6.9e-04 -17.11
5389 RPS11 19_34 1.000 2949.98 8.6e-03 -3.01
1597 PLTP 20_28 1.000 239.88 7.0e-04 -13.80
12687 RP4-781K5.7 1_121 0.999 200.54 5.8e-04 -15.25
9390 GAS6 13_62 0.995 81.45 2.4e-04 -9.23
5563 ABCG8 2_27 0.992 240.72 7.0e-04 -18.17
1999 PRKD2 19_33 0.990 28.92 8.4e-05 4.54
4704 DDX56 7_32 0.989 41.04 1.2e-04 7.99
6220 PELO 5_31 0.987 124.86 3.6e-04 10.82
6093 CSNK1G3 5_75 0.979 91.08 2.6e-04 9.47
10657 TRIM39 6_24 0.978 59.66 1.7e-04 7.87
3562 ACVR1C 2_94 0.977 29.63 8.4e-05 -5.20
5415 SYTL1 1_19 0.975 42.03 1.2e-04 -6.19
8531 TNKS 8_12 0.965 41.13 1.2e-04 8.46
11790 CYP2A6 19_28 0.960 29.27 8.2e-05 5.26
2092 SP4 7_19 0.957 105.39 2.9e-04 11.06
7040 INHBB 2_70 0.955 103.74 2.9e-04 -10.26
3212 CCND2 12_4 0.954 31.37 8.7e-05 -5.26
4608 REPS1 6_92 0.927 23.01 6.2e-05 4.14
8865 FUT2 19_33 0.927 77.86 2.1e-04 -11.86
3721 INSIG2 2_69 0.921 162.67 4.4e-04 -8.57
6778 PKN3 9_66 0.921 35.08 9.4e-05 -5.52
8579 STAT5B 17_25 0.913 25.25 6.7e-05 5.11
9072 SPTY2D1 11_13 0.892 51.45 1.3e-04 -7.12
3300 C10orf88 10_77 0.884 28.84 7.4e-05 -5.97
10621 HSPA1L 6_26 0.882 25.73 6.6e-05 4.69
1268 TMEM101 17_26 0.824 23.08 5.5e-05 4.42
3072 ADGRL2 1_50 0.814 23.05 5.5e-05 -4.25
4925 IFT172 2_16 0.810 42.77 1.0e-04 7.67
10519 SGMS1 10_33 0.809 25.60 6.0e-05 5.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
5389 RPS11 19_34 1.000 2949.98 8.6e-03 -3.01
4435 PSRC1 1_67 1.000 2883.70 8.4e-03 -55.48
1227 FLT3LG 19_34 0.000 2549.17 1.2e-07 2.72
5436 PSMA5 1_67 0.002 2087.57 1.2e-05 -47.11
4562 SRPK2 7_65 0.000 991.16 0.0e+00 -1.32
5393 RCN3 19_34 0.001 975.54 2.0e-06 3.79
1931 FCGRT 19_34 0.001 892.63 1.9e-06 3.11
10830 SYNJ2BP 14_32 0.742 761.84 1.6e-03 7.38
6970 ATXN7L2 1_67 0.002 625.61 3.9e-06 -25.64
11364 RP11-325F22.2 7_65 0.000 584.87 0.0e+00 1.95
3804 PRRG2 19_34 0.000 435.55 1.1e-07 2.69
781 PVR 19_32 0.000 434.22 0.0e+00 -12.58
11684 RP11-136O12.2 8_83 0.003 358.64 3.0e-06 18.57
5377 GEMIN7 19_32 0.000 341.98 0.0e+00 14.76
5431 SYPL2 1_67 0.007 334.16 6.4e-06 -18.18
3803 PRMT1 19_34 0.000 309.64 3.4e-08 2.88
5991 FADS1 11_34 1.000 303.61 8.9e-04 17.90
3805 SCAF1 19_34 0.000 303.59 1.6e-07 2.34
3802 IRF3 19_34 0.001 303.11 4.8e-07 2.72
4507 FADS2 11_34 0.011 267.14 8.2e-06 16.68Genes 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 2949.98 0.00860 -3.01
4435 PSRC1 1_67 1.000 2883.70 0.00840 -55.48
10830 SYNJ2BP 14_32 0.742 761.84 0.00160 7.38
5991 FADS1 11_34 1.000 303.61 0.00089 17.90
5563 ABCG8 2_27 0.992 240.72 0.00070 -18.17
1597 PLTP 20_28 1.000 239.88 0.00070 -13.80
12008 HPR 16_38 1.000 234.79 0.00069 -17.11
12687 RP4-781K5.7 1_121 0.999 200.54 0.00058 -15.25
3721 INSIG2 2_69 0.921 162.67 0.00044 -8.57
6957 USP1 1_39 0.799 155.48 0.00036 12.72
6220 PELO 5_31 0.987 124.86 0.00036 10.82
7040 INHBB 2_70 0.955 103.74 0.00029 -10.26
2092 SP4 7_19 0.957 105.39 0.00029 11.06
6093 CSNK1G3 5_75 0.979 91.08 0.00026 9.47
9390 GAS6 13_62 0.995 81.45 0.00024 -9.23
8865 FUT2 19_33 0.927 77.86 0.00021 -11.86
1932 PIH1D1 19_34 0.535 115.42 0.00018 -3.62
10657 TRIM39 6_24 0.978 59.66 0.00017 7.87
6391 TTC39B 9_13 1.000 55.22 0.00016 -7.30
10578 TMEM57 1_18 0.402 112.18 0.00013 -11.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
4435 PSRC1 1_67 1.000 2883.70 8.4e-03 -55.48
5436 PSMA5 1_67 0.002 2087.57 1.2e-05 -47.11
6970 ATXN7L2 1_67 0.002 625.61 3.9e-06 -25.64
11684 RP11-136O12.2 8_83 0.003 358.64 3.0e-06 18.57
5431 SYPL2 1_67 0.007 334.16 6.4e-06 -18.18
5563 ABCG8 2_27 0.992 240.72 7.0e-04 -18.17
5991 FADS1 11_34 1.000 303.61 8.9e-04 17.90
5240 NLRC5 16_31 0.015 241.23 1.1e-05 17.13
12008 HPR 16_38 1.000 234.79 6.9e-04 -17.11
4507 FADS2 11_34 0.011 267.14 8.2e-06 16.68
7955 FEN1 11_34 0.011 267.14 8.2e-06 16.68
1930 PPP1R37 19_32 0.000 164.19 0.0e+00 -16.41
11300 APOC2 19_32 0.000 200.78 0.0e+00 -16.02
12687 RP4-781K5.7 1_121 0.999 200.54 5.8e-04 -15.25
3441 POLK 5_45 0.003 158.20 1.2e-06 15.02
5377 GEMIN7 19_32 0.000 341.98 0.0e+00 14.76
1120 CETP 16_31 0.002 145.47 7.5e-07 14.43
12637 ZNF229 19_32 0.000 132.11 0.0e+00 14.21
2465 APOA5 11_70 0.019 209.39 1.1e-05 -14.07
1597 PLTP 20_28 1.000 239.88 7.0e-04 -13.80Comparing 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.02146592#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
4435 PSRC1 1_67 1.000 2883.70 8.4e-03 -55.48
5436 PSMA5 1_67 0.002 2087.57 1.2e-05 -47.11
6970 ATXN7L2 1_67 0.002 625.61 3.9e-06 -25.64
11684 RP11-136O12.2 8_83 0.003 358.64 3.0e-06 18.57
5431 SYPL2 1_67 0.007 334.16 6.4e-06 -18.18
5563 ABCG8 2_27 0.992 240.72 7.0e-04 -18.17
5991 FADS1 11_34 1.000 303.61 8.9e-04 17.90
5240 NLRC5 16_31 0.015 241.23 1.1e-05 17.13
12008 HPR 16_38 1.000 234.79 6.9e-04 -17.11
4507 FADS2 11_34 0.011 267.14 8.2e-06 16.68
7955 FEN1 11_34 0.011 267.14 8.2e-06 16.68
1930 PPP1R37 19_32 0.000 164.19 0.0e+00 -16.41
11300 APOC2 19_32 0.000 200.78 0.0e+00 -16.02
12687 RP4-781K5.7 1_121 0.999 200.54 5.8e-04 -15.25
3441 POLK 5_45 0.003 158.20 1.2e-06 15.02
5377 GEMIN7 19_32 0.000 341.98 0.0e+00 14.76
1120 CETP 16_31 0.002 145.47 7.5e-07 14.43
12637 ZNF229 19_32 0.000 132.11 0.0e+00 14.21
2465 APOA5 11_70 0.019 209.39 1.1e-05 -14.07
1597 PLTP 20_28 1.000 239.88 7.0e-04 -13.80Locus 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_67"
genename region_tag susie_pip mu2 PVE z
4434 VAV3 1_67 0.033 28.04 2.7e-06 -2.09
1073 SLC25A24 1_67 0.002 7.21 4.0e-08 1.21
6966 FAM102B 1_67 0.002 10.92 6.3e-08 -2.14
3009 STXBP3 1_67 0.006 28.45 4.9e-07 3.77
3438 GPSM2 1_67 0.002 10.99 5.1e-08 -2.49
3437 CLCC1 1_67 0.002 18.07 8.2e-08 3.63
10286 TAF13 1_67 0.002 7.07 3.4e-08 -1.18
10955 TMEM167B 1_67 0.002 7.13 3.6e-08 -1.01
315 SARS 1_67 0.006 165.04 3.0e-06 12.91
5431 SYPL2 1_67 0.007 334.16 6.4e-06 -18.18
5436 PSMA5 1_67 0.002 2087.57 1.2e-05 -47.11
6970 ATXN7L2 1_67 0.002 625.61 3.9e-06 -25.64
4435 PSRC1 1_67 1.000 2883.70 8.4e-03 -55.48
8615 CYB561D1 1_67 0.037 197.03 2.1e-05 12.67
9259 AMIGO1 1_67 0.006 43.41 8.2e-07 -5.47
6445 GPR61 1_67 0.002 40.59 1.9e-07 4.35
587 GNAI3 1_67 0.022 43.70 2.8e-06 -4.84
7977 GSTM4 1_67 0.002 41.08 2.4e-07 7.47
10821 GSTM2 1_67 0.002 19.59 9.1e-08 4.90
4430 GSTM1 1_67 0.003 35.88 2.6e-07 6.33
4432 GSTM5 1_67 0.002 20.04 1.1e-07 4.92
4433 GSTM3 1_67 0.003 29.40 2.2e-07 -4.29
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 8_83"
genename region_tag susie_pip mu2 PVE z
11684 RP11-136O12.2 8_83 0.003 358.64 3.0e-06 18.57
7970 FAM84B 8_83 0.002 5.31 3.3e-08 0.12
11702 PCAT1 8_83 0.004 11.51 1.4e-07 1.14
11735 CASC19 8_83 0.002 6.15 4.1e-08 -0.57
10794 POU5F1B 8_83 0.003 7.31 5.5e-08 0.70
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 2_27"
genename region_tag susie_pip mu2 PVE z
12661 LINC01126 2_27 0.000 10.07 5.0e-11 0.38
2977 THADA 2_27 0.000 13.78 1.3e-10 -3.06
6208 PLEKHH2 2_27 0.000 16.54 2.0e-10 -3.02
11022 C1GALT1C1L 2_27 0.000 27.21 3.3e-10 3.27
4930 DYNC2LI1 2_27 0.000 6.37 1.3e-11 0.15
5563 ABCG8 2_27 0.992 240.72 7.0e-04 -18.17
4943 LRPPRC 2_27 0.000 12.97 8.5e-11 -0.54
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 11_34"
genename region_tag susie_pip mu2 PVE z
9982 FAM111B 11_34 0.006 5.41 9.7e-08 -0.28
7662 FAM111A 11_34 0.006 5.13 8.9e-08 0.12
2444 DTX4 11_34 0.006 5.45 9.7e-08 -0.28
10267 MPEG1 11_34 0.006 6.04 1.1e-07 0.56
7684 PATL1 11_34 0.056 27.04 4.4e-06 2.89
7687 STX3 11_34 0.006 4.88 8.2e-08 0.24
7688 MRPL16 11_34 0.010 10.17 3.1e-07 1.16
5997 MS4A2 11_34 0.043 23.15 2.9e-06 -2.53
2453 MS4A6A 11_34 0.018 15.26 8.2e-07 2.09
10924 MS4A4E 11_34 0.013 12.95 4.9e-07 1.56
7698 MS4A14 11_34 0.030 20.10 1.7e-06 -1.92
7697 MS4A7 11_34 0.006 5.17 8.8e-08 -0.39
2455 CCDC86 11_34 0.007 7.11 1.5e-07 -0.60
2456 PRPF19 11_34 0.020 17.31 9.9e-07 1.88
2457 TMEM109 11_34 0.019 16.56 9.2e-07 1.70
2480 SLC15A3 11_34 0.008 8.93 2.0e-07 1.37
2481 CD5 11_34 0.006 5.51 9.5e-08 0.57
7874 VPS37C 11_34 0.006 5.46 9.4e-08 0.62
7875 VWCE 11_34 0.006 5.97 1.1e-07 -0.74
5990 TMEM138 11_34 0.007 10.22 2.1e-07 -2.01
6902 CYB561A3 11_34 0.007 10.22 2.1e-07 -2.01
9789 TMEM216 11_34 0.007 7.16 1.6e-07 0.46
11817 RP11-286N22.8 11_34 0.006 5.94 1.1e-07 0.55
5996 CPSF7 11_34 0.006 10.32 1.8e-07 -2.36
6903 PPP1R32 11_34 0.006 6.26 1.1e-07 -1.10
11812 RP11-794G24.1 11_34 0.007 7.15 1.5e-07 -0.61
4508 TMEM258 11_34 0.048 92.98 1.3e-05 -8.87
4507 FADS2 11_34 0.011 267.14 8.2e-06 16.68
7955 FEN1 11_34 0.011 267.14 8.2e-06 16.68
5991 FADS1 11_34 1.000 303.61 8.9e-04 17.90
1196 GANAB 11_34 0.011 123.98 4.2e-06 -11.12
11004 FADS3 11_34 0.018 30.99 1.6e-06 4.31
7876 BEST1 11_34 0.007 32.35 6.4e-07 -5.28
3676 DKFZP434K028 11_34 0.008 7.93 1.8e-07 0.76
5994 INCENP 11_34 0.006 6.42 1.1e-07 -1.11
6904 ASRGL1 11_34 0.007 6.86 1.5e-07 0.05
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 16_31"
genename region_tag susie_pip mu2 PVE z
6688 CES5A 16_31 0.001 5.22 1.4e-08 -0.54
1124 GNAO1 16_31 0.001 7.04 2.4e-08 -0.60
11561 RP11-461O7.1 16_31 0.001 8.75 3.8e-08 0.36
6695 AMFR 16_31 0.001 6.71 2.1e-08 0.01
7710 NUDT21 16_31 0.001 7.12 2.2e-08 -0.95
3681 BBS2 16_31 0.006 25.32 4.6e-07 -2.37
1122 MT3 16_31 0.001 5.18 1.4e-08 0.74
8094 MT1E 16_31 0.001 5.45 1.4e-08 0.65
10727 MT1M 16_31 0.004 25.96 3.1e-07 3.25
10725 MT1A 16_31 0.002 14.39 7.6e-08 2.18
10386 MT1F 16_31 0.004 26.92 3.4e-07 -3.20
9805 MT1X 16_31 0.001 5.43 1.5e-08 -0.05
1740 NUP93 16_31 0.006 26.70 4.7e-07 2.96
438 HERPUD1 16_31 0.002 27.93 1.3e-07 3.84
1120 CETP 16_31 0.002 145.47 7.5e-07 14.43
5240 NLRC5 16_31 0.015 241.23 1.1e-05 17.13
5239 CPNE2 16_31 0.001 7.42 2.5e-08 0.66
8472 FAM192A 16_31 0.001 6.70 2.0e-08 -0.86
6698 RSPRY1 16_31 0.001 13.16 5.7e-08 -2.21
1745 PLLP 16_31 0.002 15.89 8.4e-08 -2.79
81 CX3CL1 16_31 0.001 8.53 3.0e-08 -1.13
1747 CCL17 16_31 0.002 12.88 8.4e-08 0.99
52 CIAPIN1 16_31 0.005 22.14 3.4e-07 -1.99
1154 COQ9 16_31 0.002 10.56 5.1e-08 -1.02
3685 DOK4 16_31 0.001 6.83 2.1e-08 -0.89
4628 CCDC102A 16_31 0.001 5.13 1.4e-08 0.40
10722 ADGRG1 16_31 0.009 26.75 6.8e-07 -2.22
9366 ADGRG3 16_31 0.001 5.22 1.4e-08 -0.28
5241 KATNB1 16_31 0.002 14.05 9.4e-08 -1.39
5242 KIFC3 16_31 0.009 26.44 6.6e-07 -2.16
1754 USB1 16_31 0.001 5.41 1.5e-08 0.35
7571 ZNF319 16_31 0.001 5.05 1.3e-08 0.21
1753 MMP15 16_31 0.003 17.63 1.7e-07 -1.66
729 CFAP20 16_31 0.001 7.12 2.3e-08 0.69
730 CSNK2A2 16_31 0.001 6.36 1.9e-08 -0.56
9278 GINS3 16_31 0.001 8.83 3.4e-08 -0.91
1757 NDRG4 16_31 0.001 5.60 1.6e-08 0.39
3680 CNOT1 16_31 0.013 27.37 1.1e-06 -2.38
1759 SLC38A7 16_31 0.002 11.27 6.2e-08 1.25
3684 GOT2 16_31 0.004 18.20 2.1e-07 1.77
| 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
7595 rs79598313 1_18 1.000 114.37 3.3e-04 11.30
14109 rs11580527 1_34 1.000 96.07 2.8e-04 -11.12
14142 rs57498787 1_34 1.000 121.93 3.6e-04 -10.04
14150 rs10888896 1_34 1.000 132.40 3.9e-04 11.62
69373 rs11679386 2_12 1.000 237.04 6.9e-04 15.27
69422 rs1042034 2_13 1.000 682.94 2.0e-03 25.92
69428 rs934197 2_13 1.000 366.46 1.1e-03 36.28
69431 rs548145 2_13 1.000 1268.60 3.7e-03 43.00
69508 rs1848922 2_13 1.000 442.39 1.3e-03 33.23
76482 rs72800939 2_28 1.000 46.65 1.4e-04 -7.00
318407 rs11376017 6_13 1.000 64.43 1.9e-04 -8.44
322493 rs454182 6_22 1.000 86.78 2.5e-04 4.48
346338 rs9496567 6_67 1.000 47.39 1.4e-04 -7.07
364566 rs12208357 6_103 1.000 285.29 8.3e-04 13.53
364714 rs117733303 6_104 1.000 115.46 3.4e-04 10.35
364750 rs56393506 6_104 1.000 104.82 3.1e-04 14.79
401612 rs763798411 7_65 1.000 6915.66 2.0e-02 -3.84
428497 rs75835816 8_21 1.000 44.32 1.3e-04 6.88
438023 rs140753685 8_42 1.000 50.19 1.5e-04 7.39
439419 rs4738679 8_45 1.000 89.04 2.6e-04 -10.43
459082 rs13252684 8_83 1.000 264.31 7.7e-04 12.25
500227 rs115478735 9_70 1.000 160.03 4.7e-04 13.95
537966 rs17875416 10_71 1.000 41.71 1.2e-04 -6.65
582633 rs4937122 11_77 1.000 79.51 2.3e-04 12.93
603061 rs7397189 12_36 1.000 51.90 1.5e-04 -7.60
623737 rs11057830 12_76 1.000 35.73 1.0e-04 6.05
750198 rs1801689 17_38 1.000 92.13 2.7e-04 10.00
751114 rs113408695 17_39 1.000 123.90 3.6e-04 11.55
751140 rs8070232 17_39 1.000 142.57 4.2e-04 -7.51
784977 rs73013176 19_9 1.000 224.71 6.6e-04 -15.86
785015 rs137992968 19_9 1.000 228.66 6.7e-04 -10.90
785020 rs775693326 19_9 1.000 597.39 1.7e-03 -28.47
787812 rs3794991 19_15 1.000 338.22 9.9e-04 -18.65
794834 rs73036721 19_30 1.000 122.61 3.6e-04 -10.69
794879 rs62115478 19_30 1.000 310.89 9.1e-04 -17.85
795120 rs62117204 19_32 1.000 1463.85 4.3e-03 -59.82
795138 rs111794050 19_32 1.000 1408.27 4.1e-03 -45.67
795171 rs814573 19_32 1.000 4067.72 1.2e-02 74.70
795173 rs113345881 19_32 1.000 1607.93 4.7e-03 -48.87
795176 rs12721109 19_32 1.000 2346.00 6.8e-03 -60.87
805449 rs6075251 20_13 1.000 88.40 2.6e-04 -3.56
805450 rs34507316 20_13 1.000 116.57 3.4e-04 -7.80
874519 rs34287152 1_67 1.000 155.33 4.5e-04 -1.67
895911 rs10677712 2_69 1.000 6557.67 1.9e-02 5.16
999845 rs41274050 10_33 1.000 55.27 1.6e-04 7.99
1027259 rs964184 11_70 1.000 414.80 1.2e-03 -21.68
1057956 rs559100757 14_32 1.000 1646.63 4.8e-03 -1.67
1111714 rs374141296 19_34 1.000 2693.85 7.9e-03 3.02
14157 rs471705 1_34 0.999 142.99 4.2e-04 15.46
401623 rs4997569 7_65 0.999 6951.29 2.0e-02 -3.54
876109 rs140584594 1_67 0.999 63.41 1.8e-04 -11.12
69425 rs78610189 2_13 0.998 111.02 3.2e-04 -10.48
785060 rs2043302 19_10 0.998 47.31 1.4e-04 5.05
1057941 rs35013494 14_32 0.998 1673.13 4.9e-03 -3.54
785095 rs201868221 19_10 0.997 51.27 1.5e-04 7.69
787843 rs113619686 19_15 0.997 46.17 1.3e-04 0.73
794824 rs769162207 19_30 0.997 38.50 1.1e-04 0.61
279809 rs7701166 5_45 0.996 32.39 9.4e-05 -2.82
785008 rs2738447 19_9 0.996 381.28 1.1e-03 18.67
1118156 rs1800961 20_28 0.996 34.05 9.9e-05 -5.94
607427 rs148481241 12_44 0.995 28.49 8.3e-05 5.22
784998 rs147985405 19_9 0.993 1633.67 4.7e-03 -47.76
28857 rs1109112 1_69 0.991 28.64 8.3e-05 -5.03
427750 rs1495743 8_20 0.990 33.84 9.8e-05 -5.88
785041 rs4804149 19_10 0.988 35.24 1.0e-04 7.46
890286 rs1260326 2_16 0.988 274.24 7.9e-04 -20.45
801095 rs74273659 20_5 0.987 33.78 9.7e-05 6.07
810407 rs73124945 20_24 0.987 32.03 9.2e-05 -7.68
53576 rs2807848 1_112 0.986 30.67 8.8e-05 -6.39
28859 rs78221564 1_69 0.985 27.69 8.0e-05 -4.76
674821 rs13379043 14_34 0.984 28.88 8.3e-05 -5.43
144638 rs9834932 3_24 0.982 65.59 1.9e-04 -8.42
546829 rs10838525 11_4 0.982 41.27 1.2e-04 -5.31
737199 rs144129583 17_7 0.982 29.04 8.3e-05 -5.94
632624 rs1012130 13_10 0.981 71.91 2.1e-04 -3.82
279750 rs10062361 5_45 0.980 162.83 4.7e-04 17.92
924546 rs535137438 5_31 0.979 31.72 9.1e-05 -5.54
787452 rs2302209 19_14 0.978 34.61 9.9e-05 5.78
787861 rs12984303 19_15 0.977 26.54 7.6e-05 5.24
726521 rs4396539 16_37 0.976 27.19 7.7e-05 -5.06
810354 rs6029132 20_24 0.976 40.89 1.2e-04 -6.92
324198 rs112357706 6_27 0.964 25.11 7.1e-05 5.21
632629 rs206326 13_10 0.962 49.28 1.4e-04 -5.04
1066387 rs5880 16_31 0.962 56.00 1.6e-04 10.84
619647 rs653178 12_67 0.959 56.36 1.6e-04 8.65
576350 rs201912654 11_59 0.957 35.34 9.9e-05 -5.86
195027 rs3748034 4_4 0.956 51.27 1.4e-04 6.91
349074 rs12199109 6_73 0.955 25.58 7.1e-05 5.19
459071 rs79658059 8_83 0.954 395.14 1.1e-03 -20.19
322901 rs28986304 6_23 0.952 37.25 1.0e-04 6.87
537676 rs12244851 10_70 0.952 30.00 8.3e-05 -4.33
805430 rs78348000 20_13 0.938 33.38 9.1e-05 5.58
221631 rs1458038 4_54 0.937 38.77 1.1e-04 -6.28
322930 rs3130253 6_23 0.926 26.09 7.1e-05 5.17
1066353 rs5883 16_31 0.926 35.27 9.5e-05 -3.81
400146 rs3197597 7_61 0.925 25.20 6.8e-05 -4.71
582620 rs1614592 11_77 0.923 27.34 7.4e-05 -6.08
378698 rs10268799 7_23 0.918 27.97 7.5e-05 5.33
810372 rs6102034 20_24 0.913 76.26 2.0e-04 -10.11
364560 rs9456502 6_103 0.912 31.87 8.5e-05 6.00
810403 rs76981217 20_24 0.912 29.35 7.8e-05 6.60
272691 rs55681913 5_28 0.909 24.97 6.6e-05 -4.97
137557 rs307572 3_9 0.907 26.28 7.0e-05 -5.25
632616 rs1799955 13_10 0.907 131.87 3.5e-04 -9.26
14164 rs10493176 1_34 0.901 94.27 2.5e-04 -12.21
439387 rs56386732 8_45 0.899 27.96 7.3e-05 -5.92
1111736 rs142385484 19_34 0.897 51.54 1.3e-04 -5.32
895907 rs28850371 2_69 0.896 6552.20 1.7e-02 4.81
732934 rs77277579 16_52 0.894 25.75 6.7e-05 -4.83
531855 rs10882161 10_59 0.877 40.98 1.0e-04 -6.50
582636 rs74612335 11_77 0.877 89.26 2.3e-04 12.89
645285 rs9564985 13_36 0.865 25.57 6.5e-05 -4.75
985091 rs4841181 8_12 0.857 34.35 8.6e-05 -5.75
798905 rs34003091 19_39 0.856 97.81 2.4e-04 -10.07
279773 rs3843482 5_45 0.849 310.77 7.7e-04 21.88
837739 rs145678077 22_17 0.847 24.92 6.2e-05 -4.78
737090 rs402614 17_6 0.837 24.28 5.9e-05 4.09
660463 rs2332328 14_3 0.836 55.58 1.4e-04 7.75
751125 rs9303012 17_39 0.835 137.52 3.4e-04 2.40
745703 rs7212325 17_28 0.834 37.98 9.2e-05 -7.37
622602 rs1169300 12_74 0.828 56.34 1.4e-04 7.69
809148 rs11167269 20_21 0.826 51.58 1.2e-04 -7.30
817678 rs816955 20_38 0.822 24.67 5.9e-05 -3.86
794915 rs79280038 19_30 0.821 42.64 1.0e-04 -6.44
364580 rs3818678 6_103 0.819 210.81 5.0e-04 -9.82
56220 rs11122453 1_117 0.817 71.07 1.7e-04 -8.77
119725 rs7569317 2_120 0.811 38.44 9.1e-05 7.72
893195 rs12990177 2_27 0.804 41.92 9.8e-05 5.94SNPs 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
401623 rs4997569 7_65 0.999 6951.29 2.0e-02 -3.54
401615 rs10274607 7_65 0.027 6932.46 5.4e-04 -3.41
401630 rs6952534 7_65 0.001 6924.86 1.6e-05 -3.46
401618 rs13230660 7_65 0.046 6924.10 9.2e-04 -3.48
401629 rs4730069 7_65 0.000 6918.38 2.0e-06 -3.44
401612 rs763798411 7_65 1.000 6915.66 2.0e-02 -3.84
401622 rs10242713 7_65 0.000 6891.47 3.2e-09 -3.36
401625 rs10249965 7_65 0.000 6834.82 8.3e-15 -3.33
895911 rs10677712 2_69 1.000 6557.67 1.9e-02 5.16
895907 rs28850371 2_69 0.896 6552.20 1.7e-02 4.81
895868 rs6542385 2_69 0.445 6546.12 8.5e-03 4.84
401637 rs1013016 7_65 0.000 6541.67 0.0e+00 2.82
895915 rs1546621 2_69 0.084 6538.35 1.6e-03 4.81
895920 rs1516006 2_69 0.047 6537.66 8.9e-04 4.79
401662 rs8180737 7_65 0.000 6224.65 0.0e+00 -3.26
401655 rs17778396 7_65 0.000 6221.96 0.0e+00 -3.22
401656 rs2237621 7_65 0.000 6219.36 0.0e+00 -3.23
401627 rs71562637 7_65 0.000 6211.55 0.0e+00 -3.13
401689 rs10224564 7_65 0.000 6208.75 0.0e+00 -3.24
401674 rs10255779 7_65 0.000 6206.68 0.0e+00 -3.28
401691 rs78132606 7_65 0.000 6175.99 0.0e+00 -3.24
401694 rs4610671 7_65 0.000 6166.63 0.0e+00 -3.18
401696 rs12669532 7_65 0.000 5913.37 0.0e+00 -3.21
401653 rs2237618 7_65 0.000 5808.57 0.0e+00 -2.88
401698 rs118089279 7_65 0.000 5760.28 0.0e+00 -3.19
401685 rs73188303 7_65 0.000 5746.51 0.0e+00 -2.85
895896 rs151256241 2_69 0.000 5276.45 1.0e-12 5.64
895916 rs2551674 2_69 0.000 5270.56 8.6e-13 5.61
895928 rs2245241 2_69 0.000 5267.21 7.6e-13 5.59
895931 rs2674843 2_69 0.000 5264.64 7.2e-13 5.58
895979 rs2255565 2_69 0.000 5249.32 6.7e-13 5.56
895905 rs36144266 2_69 0.000 5159.69 1.0e-10 5.20
895903 rs201120706 2_69 0.000 5158.23 9.8e-11 5.19
895913 rs7575129 2_69 0.000 5153.20 1.2e-10 5.23
895922 rs965193 2_69 0.000 5152.04 8.4e-11 5.18
895859 rs4117160 2_69 0.000 5148.14 1.1e-10 5.22
895777 rs1516009 2_69 0.000 5126.69 1.6e-10 5.29
895857 rs2551687 2_69 0.000 5101.40 4.9e-11 6.14
895851 rs2551685 2_69 0.000 5100.69 7.4e-11 6.22
895787 rs2551684 2_69 0.000 5083.71 8.9e-11 6.25
895795 rs11695801 2_69 0.000 4789.42 3.0e-11 5.08
896005 rs2551679 2_69 0.000 4665.12 6.3e-12 4.87
895948 rs2263642 2_69 0.000 4592.07 1.4e-07 5.88
895938 rs2216845 2_69 0.000 4591.36 1.1e-07 5.88
895986 rs2161843 2_69 0.000 4591.13 1.2e-07 5.89
895987 rs2194495 2_69 0.000 4590.75 1.2e-07 5.88
895956 rs2421964 2_69 0.000 4590.72 1.0e-07 5.86
895960 rs2244678 2_69 0.000 4590.68 1.2e-07 5.88
895957 rs2421965 2_69 0.000 4590.62 1.0e-07 5.86
895976 rs2255586 2_69 0.000 4590.61 1.2e-07 5.88SNPs 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
401612 rs763798411 7_65 1.000 6915.66 0.02000 -3.84
401623 rs4997569 7_65 0.999 6951.29 0.02000 -3.54
895911 rs10677712 2_69 1.000 6557.67 0.01900 5.16
895907 rs28850371 2_69 0.896 6552.20 0.01700 4.81
795171 rs814573 19_32 1.000 4067.72 0.01200 74.70
895868 rs6542385 2_69 0.445 6546.12 0.00850 4.84
1111714 rs374141296 19_34 1.000 2693.85 0.00790 3.02
795176 rs12721109 19_32 1.000 2346.00 0.00680 -60.87
1057941 rs35013494 14_32 0.998 1673.13 0.00490 -3.54
1057956 rs559100757 14_32 1.000 1646.63 0.00480 -1.67
784998 rs147985405 19_9 0.993 1633.67 0.00470 -47.76
795173 rs113345881 19_32 1.000 1607.93 0.00470 -48.87
795120 rs62117204 19_32 1.000 1463.85 0.00430 -59.82
795138 rs111794050 19_32 1.000 1408.27 0.00410 -45.67
69431 rs548145 2_13 1.000 1268.60 0.00370 43.00
69422 rs1042034 2_13 1.000 682.94 0.00200 25.92
785020 rs775693326 19_9 1.000 597.39 0.00170 -28.47
895915 rs1546621 2_69 0.084 6538.35 0.00160 4.81
69508 rs1848922 2_13 1.000 442.39 0.00130 33.23
1027259 rs964184 11_70 1.000 414.80 0.00120 -21.68
69428 rs934197 2_13 1.000 366.46 0.00110 36.28
459071 rs79658059 8_83 0.954 395.14 0.00110 -20.19
785008 rs2738447 19_9 0.996 381.28 0.00110 18.67
787812 rs3794991 19_15 1.000 338.22 0.00099 -18.65
401618 rs13230660 7_65 0.046 6924.10 0.00092 -3.48
794879 rs62115478 19_30 1.000 310.89 0.00091 -17.85
895920 rs1516006 2_69 0.047 6537.66 0.00089 4.79
364566 rs12208357 6_103 1.000 285.29 0.00083 13.53
890286 rs1260326 2_16 0.988 274.24 0.00079 -20.45
279773 rs3843482 5_45 0.849 310.77 0.00077 21.88
459082 rs13252684 8_83 1.000 264.31 0.00077 12.25
1057123 rs8021180 14_32 0.333 738.39 0.00072 7.90
69373 rs11679386 2_12 1.000 237.04 0.00069 15.27
785015 rs137992968 19_9 1.000 228.66 0.00067 -10.90
784977 rs73013176 19_9 1.000 224.71 0.00066 -15.86
1057298 rs1044527 14_32 0.286 740.08 0.00062 7.88
401615 rs10274607 7_65 0.027 6932.46 0.00054 -3.41
364580 rs3818678 6_103 0.819 210.81 0.00050 -9.82
279750 rs10062361 5_45 0.980 162.83 0.00047 17.92
500227 rs115478735 9_70 1.000 160.03 0.00047 13.95
874519 rs34287152 1_67 1.000 155.33 0.00045 -1.67
14157 rs471705 1_34 0.999 142.99 0.00042 15.46
751140 rs8070232 17_39 1.000 142.57 0.00042 -7.51
14150 rs10888896 1_34 1.000 132.40 0.00039 11.62
459070 rs2980875 8_83 0.493 269.45 0.00039 -26.62
1057062 rs12879801 14_32 0.175 731.45 0.00037 7.83
14142 rs57498787 1_34 1.000 121.93 0.00036 -10.04
751114 rs113408695 17_39 1.000 123.90 0.00036 11.55
794834 rs73036721 19_30 1.000 122.61 0.00036 -10.69
632616 rs1799955 13_10 0.907 131.87 0.00035 -9.26SNPs with largest z scores
#SNPs with 50 largest z scores
head(ctwas_snp_res[order(-abs(ctwas_snp_res$z)),report_cols_snps],50) id region_tag susie_pip mu2 PVE z
795171 rs814573 19_32 1.000 4067.72 1.2e-02 74.70
795176 rs12721109 19_32 1.000 2346.00 6.8e-03 -60.87
795120 rs62117204 19_32 1.000 1463.85 4.3e-03 -59.82
795107 rs1551891 19_32 0.000 830.33 0.0e+00 -56.10
875232 rs12740374 1_67 0.000 2612.94 1.5e-06 -55.62
875239 rs646776 1_67 0.000 2599.97 1.0e-06 55.51
875238 rs629301 1_67 0.000 2596.23 8.9e-07 55.48
875228 rs7528419 1_67 0.000 2588.41 1.1e-06 -55.38
875250 rs583104 1_67 0.000 2510.61 8.4e-07 54.58
875253 rs4970836 1_67 0.000 2509.60 8.6e-07 54.55
875255 rs1277930 1_67 0.000 2497.79 8.5e-07 54.42
875256 rs599839 1_67 0.000 2492.67 8.5e-07 54.37
875236 rs3832016 1_67 0.000 2447.99 5.0e-07 53.86
875233 rs660240 1_67 0.000 2435.30 5.0e-07 53.72
875251 rs602633 1_67 0.000 2384.92 5.5e-07 53.17
795167 rs405509 19_32 0.000 2307.45 0.0e+00 -49.85
795173 rs113345881 19_32 1.000 1607.93 4.7e-03 -48.87
784998 rs147985405 19_9 0.993 1633.67 4.7e-03 -47.76
784991 rs138175288 19_9 0.002 1620.88 7.2e-06 -47.58
784993 rs73015020 19_9 0.004 1623.52 2.1e-05 -47.58
784992 rs138294113 19_9 0.000 1616.96 2.0e-06 -47.55
784995 rs112552009 19_9 0.000 1612.47 4.8e-07 -47.52
784994 rs77140532 19_9 0.000 1618.43 2.2e-06 -47.51
784996 rs10412048 19_9 0.000 1615.89 7.0e-07 -47.47
784990 rs55997232 19_9 0.000 1597.90 8.3e-10 -47.33
875219 rs4970834 1_67 0.000 1817.97 2.4e-06 -46.39
795138 rs111794050 19_32 1.000 1408.27 4.1e-03 -45.67
795091 rs62120566 19_32 0.000 2314.50 0.0e+00 -45.03
795085 rs188099946 19_32 0.000 2171.00 0.0e+00 -43.85
795144 rs4802238 19_32 0.000 1528.78 0.0e+00 43.14
69431 rs548145 2_13 1.000 1268.60 3.7e-03 43.00
875240 rs3902354 1_67 0.000 1535.18 3.3e-07 42.81
875229 rs11102967 1_67 0.000 1531.11 3.4e-07 42.74
875254 rs4970837 1_67 0.000 1516.12 3.9e-07 42.55
795079 rs201314191 19_32 0.000 1992.47 0.0e+00 -42.45
795146 rs56394238 19_32 0.000 1716.05 0.0e+00 42.20
795155 rs2972559 19_32 0.000 2087.44 0.0e+00 42.06
875224 rs611917 1_67 0.000 1444.10 2.9e-07 -41.46
795123 rs2965169 19_32 0.000 522.86 0.0e+00 -41.16
795147 rs3021439 19_32 0.000 1330.83 0.0e+00 40.36
795084 rs62119327 19_32 0.000 1762.68 0.0e+00 -40.32
69432 rs478588 2_13 0.000 1179.97 0.0e+00 39.96
795154 rs12162222 19_32 0.000 1798.87 0.0e+00 39.84
784999 rs17248769 19_9 0.000 1090.74 1.5e-11 -39.67
785000 rs2228671 19_9 0.000 1080.01 7.6e-12 -39.52
69459 rs532300 2_13 0.000 1083.48 0.0e+00 38.57
69460 rs558130 2_13 0.000 1083.47 0.0e+00 38.57
69461 rs533211 2_13 0.000 1083.47 0.0e+00 38.57
69457 rs312979 2_13 0.000 1082.55 0.0e+00 38.56
69462 rs528113 2_13 0.000 1083.02 0.0e+00 38.56Gene 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] 33if (length(genes)>0){
GO_enrichment <- enrichr(genes, dbs)
for (db in dbs){
print(db)
df <- GO_enrichment[[db]]
df <- df[df$Adjusted.P.value<0.05,c("Term", "Overlap", "Adjusted.P.value", "Genes")]
print(df)
}
#DisGeNET enrichment
# devtools::install_bitbucket("ibi_group/disgenet2r")
library(disgenet2r)
disgenet_api_key <- get_disgenet_api_key(
email = "wesleycrouse@gmail.com",
password = "uchicago1" )
Sys.setenv(DISGENET_API_KEY= disgenet_api_key)
res_enrich <-disease_enrichment(entities=genes, vocabulary = "HGNC",
database = "CURATED" )
df <- res_enrich@qresult[1:10, c("Description", "FDR", "Ratio", "BgRatio")]
print(df)
#WebGestalt enrichment
library(WebGestaltR)
background <- ctwas_gene_res$genename
#listGeneSet()
databases <- c("pathway_KEGG", "disease_GLAD4U", "disease_OMIM")
enrichResult <- WebGestaltR(enrichMethod="ORA", organism="hsapiens",
interestGene=genes, referenceGene=background,
enrichDatabase=databases, interestGeneType="genesymbol",
referenceGeneType="genesymbol", isOutput=F)
print(enrichResult[,c("description", "size", "overlap", "FDR", "database", "userId")])
}Uploading data to Enrichr... Done.
Querying GO_Biological_Process_2021... Done.
Querying GO_Cellular_Component_2021... Done.
Querying GO_Molecular_Function_2021... Done.
Parsing results... Done.
[1] "GO_Biological_Process_2021"
Term
1 peptidyl-serine phosphorylation (GO:0018105)
2 peptidyl-serine modification (GO:0018209)
3 protein phosphorylation (GO:0006468)
4 positive regulation of cyclin-dependent protein serine/threonine kinase activity (GO:0045737)
5 activin receptor signaling pathway (GO:0032924)
6 positive regulation of cyclin-dependent protein kinase activity (GO:1904031)
Overlap Adjusted.P.value Genes
1 5/156 0.001794023 CSNK1G3;TNKS;PKN3;PRKD2;GAS6
2 5/169 0.001794023 CSNK1G3;TNKS;PKN3;PRKD2;GAS6
3 6/496 0.021266201 CSNK1G3;ACVR1C;TNKS;PKN3;PRKD2;GAS6
4 2/17 0.037009718 CCND2;PSRC1
5 2/19 0.037009718 ACVR1C;INHBB
6 2/20 0.037009718 CCND2;PSRC1
[1] "GO_Cellular_Component_2021"
Term Overlap
1 high-density lipoprotein particle (GO:0034364) 2/19
2 serine/threonine protein kinase complex (GO:1902554) 2/37
Adjusted.P.value Genes
1 0.02262305 HPR;PLTP
2 0.04324534 ACVR1C;CCND2
[1] "GO_Molecular_Function_2021"
Term Overlap Adjusted.P.value
1 cholesterol transfer activity (GO:0120020) 2/18 0.01619101
2 sterol transfer activity (GO:0120015) 2/19 0.01619101
Genes
1 ABCG8;PLTP
2 ABCG8;PLTPRP4-781K5.7 gene(s) from the input list not found in DisGeNET CURATEDPELO gene(s) from the input list not found in DisGeNET CURATEDTNKS gene(s) from the input list not found in DisGeNET CURATEDTTC39B gene(s) from the input list not found in DisGeNET CURATEDPKN3 gene(s) from the input list not found in DisGeNET CURATEDTRIM39 gene(s) from the input list not found in DisGeNET CURATEDADGRL2 gene(s) from the input list not found in DisGeNET CURATEDSPTY2D1 gene(s) from the input list not found in DisGeNET CURATEDCSNK1G3 gene(s) from the input list not found in DisGeNET CURATEDTMEM101 gene(s) from the input list not found in DisGeNET CURATEDDDX56 gene(s) from the input list not found in DisGeNET CURATEDSGMS1 gene(s) from the input list not found in DisGeNET CURATEDHPR gene(s) from the input list not found in DisGeNET CURATEDSYTL1 gene(s) from the input list not found in DisGeNET CURATEDRPS11 gene(s) from the input list not found in DisGeNET CURATEDC10orf88 gene(s) from the input list not found in DisGeNET CURATEDPSRC1 gene(s) from the input list not found in DisGeNET CURATED Description FDR
36 Leukemia, T-Cell, Chronic 0.01502554
47 Opisthorchiasis 0.01502554
72 Caliciviridae Infections 0.01502554
78 Infections, Calicivirus 0.01502554
96 Opisthorchis felineus Infection 0.01502554
97 Opisthorchis viverrini Infection 0.01502554
108 Enteropathy-Associated T-Cell Lymphoma 0.01502554
126 Leukemia, Large Granular Lymphocytic 0.01502554
133 Laron syndrome type 2 0.01502554
139 Leukemia, Natural Killer Cell Large Granular Lymphocytic 0.01502554
Ratio BgRatio
36 1/16 1/9703
47 1/16 1/9703
72 1/16 1/9703
78 1/16 1/9703
96 1/16 1/9703
97 1/16 1/9703
108 1/16 1/9703
126 1/16 1/9703
133 1/16 1/9703
139 1/16 1/9703******************************************
* *
* Welcome to WebGestaltR ! *
* *
******************************************Loading the functional categories...
Loading the ID list...
Loading the reference list...
Performing the enrichment analysis...
description size overlap FDR database
1 Coronary Artery Disease 153 8 0.0001162333 disease_GLAD4U
2 Dyslipidaemia 84 6 0.0005365306 disease_GLAD4U
3 Coronary Disease 171 7 0.0015732547 disease_GLAD4U
4 Arteriosclerosis 173 6 0.0179513410 disease_GLAD4U
5 Myocardial Ischemia 181 6 0.0185276913 disease_GLAD4U
userId
1 PSRC1;ABCG8;INSIG2;TTC39B;SPTY2D1;FADS1;FUT2;PLTP
2 PSRC1;ABCG8;INSIG2;TTC39B;FADS1;PLTP
3 PSRC1;ABCG8;INSIG2;TTC39B;FADS1;FUT2;PLTP
4 PSRC1;ABCG8;TTC39B;FADS1;HPR;PLTP
5 PSRC1;ABCG8;INSIG2;TTC39B;FADS1;PLTP
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