Alanine aminotransferase (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
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-30620_irnt_Liver.Rmd) and HTML (docs/ukb-d-30620_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 Alanine aminotransferase (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-30620_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.0123608448 0.0002075549 #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
29.46131 11.32261 #report sample size
print(sample_size)[1] 344136#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.01153550 0.05939304 #compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))[1] 0.0273382 0.5979644Genes 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
6291 JAZF1 7_23 1.000 69.93 2.0e-04 -8.32
2463 PANX1 11_53 1.000 1457.10 4.2e-03 16.23
12135 S1PR2 19_9 1.000 69.50 2.0e-04 -7.49
2924 EFHD1 2_136 0.999 160.01 4.6e-04 12.80
2204 AKNA 9_59 0.998 223.50 6.5e-04 19.72
3212 CCND2 12_4 0.990 41.74 1.2e-04 -6.27
11790 CYP2A6 19_28 0.990 42.75 1.2e-04 6.55
9482 ACTG1 17_46 0.979 108.09 3.1e-04 10.51
5563 ABCG8 2_27 0.977 55.84 1.6e-04 8.07
10912 STARD10 11_41 0.975 73.00 2.1e-04 9.11
9404 PTTG1IP 21_23 0.974 48.92 1.4e-04 7.00
2985 PLEKHA3 2_108 0.971 37.68 1.1e-04 -5.69
5415 SYTL1 1_19 0.970 32.49 9.2e-05 -5.44
8803 DLEU1 13_21 0.970 44.65 1.3e-04 6.43
9985 LITAF 16_12 0.962 23.98 6.7e-05 -4.19
12467 RP11-219B17.3 15_27 0.948 54.97 1.5e-04 -7.27
6171 ARL14EP 11_21 0.946 278.61 7.7e-04 -5.33
7706 EVA1C 21_13 0.934 34.20 9.3e-05 5.55
1848 CD276 15_35 0.932 69.59 1.9e-04 8.07
4239 TRIM5 11_4 0.930 43.36 1.2e-04 -5.20
676 IFT80 3_99 0.911 81.60 2.2e-04 8.77
2475 NECTIN1 11_72 0.906 21.51 5.7e-05 -4.15
2004 TGFB1 19_28 0.897 26.70 7.0e-05 4.67
9783 RHD 1_18 0.886 29.25 7.5e-05 5.07
2801 TFDP2 3_87 0.865 33.24 8.4e-05 5.66
474 TAB2 6_97 0.864 20.85 5.2e-05 3.92
2341 DDX5 17_37 0.863 20.58 5.2e-05 -3.68
5769 MLIP 6_40 0.831 35.15 8.5e-05 -5.85
11198 RP6-109B7.2 22_20 0.826 21.70 5.2e-05 -4.08
3291 SLF2 10_64 0.818 29.34 7.0e-05 4.78
5436 PSMA5 1_67 0.810 32.94 7.8e-05 -5.07
11072 PTPRD-AS1 9_9 0.807 21.33 5.0e-05 -3.88Genes 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
4556 TMEM60 7_49 0.000 3027.40 1.6e-14 3.85
2463 PANX1 11_53 1.000 1457.10 4.2e-03 16.23
9370 C11orf54 11_53 0.000 1015.23 0.0e+00 -1.41
10840 PPP1CB 2_17 0.000 631.98 1.0e-10 -2.87
10903 APTR 7_49 0.000 600.12 0.0e+00 2.36
8270 TRMT61B 2_17 0.000 488.64 2.5e-15 2.54
1320 CWF19L1 10_64 0.001 457.96 1.4e-06 -24.16
7172 SPDYA 2_17 0.000 425.56 1.5e-18 2.33
258 MRE11 11_53 0.000 396.86 0.0e+00 2.05
9811 RSBN1L 7_49 0.000 362.21 0.0e+00 2.97
6171 ARL14EP 11_21 0.946 278.61 7.7e-04 -5.33
2204 AKNA 9_59 0.998 223.50 6.5e-04 19.72
3308 CPN1 10_64 0.000 214.32 8.1e-08 17.98
92 PHTF2 7_49 0.000 196.06 0.0e+00 -0.08
11684 RP11-136O12.2 8_83 0.002 171.41 9.6e-07 15.32
2924 EFHD1 2_136 0.999 160.01 4.6e-04 12.80
5400 EPHA2 1_11 0.011 158.11 5.0e-06 -13.66
402 USP28 11_67 0.002 139.86 9.3e-07 -0.70
11699 RP11-10A14.4 8_11 0.000 136.93 1.4e-10 3.60
5872 NAPRT 8_94 0.001 126.21 2.2e-07 -5.21Genes 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
2463 PANX1 11_53 1.000 1457.10 0.00420 16.23
6171 ARL14EP 11_21 0.946 278.61 0.00077 -5.33
2204 AKNA 9_59 0.998 223.50 0.00065 19.72
2924 EFHD1 2_136 0.999 160.01 0.00046 12.80
9482 ACTG1 17_46 0.979 108.09 0.00031 10.51
676 IFT80 3_99 0.911 81.60 0.00022 8.77
10912 STARD10 11_41 0.975 73.00 0.00021 9.11
6291 JAZF1 7_23 1.000 69.93 0.00020 -8.32
12135 S1PR2 19_9 1.000 69.50 0.00020 -7.49
11738 RP11-115J16.2 8_12 0.725 89.96 0.00019 -10.54
1848 CD276 15_35 0.932 69.59 0.00019 8.07
5563 ABCG8 2_27 0.977 55.84 0.00016 8.07
12467 RP11-219B17.3 15_27 0.948 54.97 0.00015 -7.27
9404 PTTG1IP 21_23 0.974 48.92 0.00014 7.00
6831 RPL8 8_94 0.514 90.13 0.00013 -8.73
8803 DLEU1 13_21 0.970 44.65 0.00013 6.43
4239 TRIM5 11_4 0.930 43.36 0.00012 -5.20
3212 CCND2 12_4 0.990 41.74 0.00012 -6.27
11790 CYP2A6 19_28 0.990 42.75 0.00012 6.55
8148 SPDYE5 7_48 0.787 48.20 0.00011 6.83Genes 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
1320 CWF19L1 10_64 0.001 457.96 1.4e-06 -24.16
2204 AKNA 9_59 0.998 223.50 6.5e-04 19.72
3308 CPN1 10_64 0.000 214.32 8.1e-08 17.98
2463 PANX1 11_53 1.000 1457.10 4.2e-03 16.23
11684 RP11-136O12.2 8_83 0.002 171.41 9.6e-07 15.32
5400 EPHA2 1_11 0.011 158.11 5.0e-06 -13.66
2924 EFHD1 2_136 0.999 160.01 4.6e-04 12.80
281 ABCC2 10_64 0.111 78.90 2.5e-05 11.84
6478 HKDC1 10_46 0.012 119.36 4.2e-06 -11.01
11738 RP11-115J16.2 8_12 0.725 89.96 1.9e-04 -10.54
9482 ACTG1 17_46 0.979 108.09 3.1e-04 10.51
6086 DLG5 10_50 0.014 77.33 3.1e-06 10.46
9761 FSCN2 17_46 0.027 101.11 8.0e-06 10.07
4319 RSG1 1_11 0.006 86.43 1.5e-06 -9.38
10912 STARD10 11_41 0.975 73.00 2.1e-04 9.11
8267 SOX7 8_14 0.551 57.93 9.3e-05 -9.08
10066 ZNF34 8_94 0.020 84.38 5.0e-06 9.04
2536 SH2B3 12_67 0.575 47.72 8.0e-05 9.03
2237 ERLIN1 10_64 0.002 96.09 5.7e-07 -8.84
676 IFT80 3_99 0.911 81.60 2.2e-04 8.77Comparing 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.02054857#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
1320 CWF19L1 10_64 0.001 457.96 1.4e-06 -24.16
2204 AKNA 9_59 0.998 223.50 6.5e-04 19.72
3308 CPN1 10_64 0.000 214.32 8.1e-08 17.98
2463 PANX1 11_53 1.000 1457.10 4.2e-03 16.23
11684 RP11-136O12.2 8_83 0.002 171.41 9.6e-07 15.32
5400 EPHA2 1_11 0.011 158.11 5.0e-06 -13.66
2924 EFHD1 2_136 0.999 160.01 4.6e-04 12.80
281 ABCC2 10_64 0.111 78.90 2.5e-05 11.84
6478 HKDC1 10_46 0.012 119.36 4.2e-06 -11.01
11738 RP11-115J16.2 8_12 0.725 89.96 1.9e-04 -10.54
9482 ACTG1 17_46 0.979 108.09 3.1e-04 10.51
6086 DLG5 10_50 0.014 77.33 3.1e-06 10.46
9761 FSCN2 17_46 0.027 101.11 8.0e-06 10.07
4319 RSG1 1_11 0.006 86.43 1.5e-06 -9.38
10912 STARD10 11_41 0.975 73.00 2.1e-04 9.11
8267 SOX7 8_14 0.551 57.93 9.3e-05 -9.08
10066 ZNF34 8_94 0.020 84.38 5.0e-06 9.04
2536 SH2B3 12_67 0.575 47.72 8.0e-05 9.03
2237 ERLIN1 10_64 0.002 96.09 5.7e-07 -8.84
676 IFT80 3_99 0.911 81.60 2.2e-04 8.77Locus 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: 10_64"
genename region_tag susie_pip mu2 PVE z
3299 CNNM1 10_64 0.000 11.36 4.5e-09 1.94
3307 GOT1 10_64 0.000 19.63 9.6e-09 4.11
11056 RP11-441O15.3 10_64 0.000 19.63 9.6e-09 -4.11
11947 RP11-85A1.3 10_64 0.000 11.08 2.6e-09 2.31
10330 ENTPD7 10_64 0.000 28.10 1.9e-08 -2.81
3296 CUTC 10_64 0.000 27.01 3.2e-08 -2.25
228 COX15 10_64 0.000 18.56 2.5e-08 1.34
281 ABCC2 10_64 0.111 78.90 2.5e-05 11.84
2234 DNMBP 10_64 0.000 41.70 2.4e-08 -7.04
3308 CPN1 10_64 0.000 214.32 8.1e-08 17.98
2237 ERLIN1 10_64 0.002 96.09 5.7e-07 -8.84
10819 CHUK 10_64 0.001 59.68 1.1e-07 -7.11
1320 CWF19L1 10_64 0.001 457.96 1.4e-06 -24.16
10014 BLOC1S2 10_64 0.001 87.02 2.1e-07 -8.73
11326 OLMALINC 10_64 0.000 18.04 8.9e-09 -2.20
12405 RP11-285F16.1 10_64 0.000 5.40 1.3e-09 -0.35
7557 NDUFB8 10_64 0.000 5.45 1.3e-09 -0.10
3291 SLF2 10_64 0.818 29.34 7.0e-05 4.78
1321 SEMA4G 10_64 0.000 8.05 2.4e-09 -0.70
2256 LZTS2 10_64 0.043 23.05 2.8e-06 3.82
9772 PDZD7 10_64 0.000 5.28 1.2e-09 0.50
2254 TLX1 10_64 0.000 7.78 2.5e-09 0.23
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 9_59"
genename region_tag susie_pip mu2 PVE z
1319 WHRN 9_59 0.007 25.13 4.8e-07 -5.37
2204 AKNA 9_59 0.998 223.50 6.5e-04 19.72
4774 ATP6V1G1 9_59 0.008 7.55 1.7e-07 0.01
6550 TMEM268 9_59 0.006 16.08 2.8e-07 3.91
375 TNC 9_59 0.007 8.28 1.7e-07 1.15
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 11_53"
genename region_tag susie_pip mu2 PVE z
7542 CEP295 11_53 0 53.80 0.0000 0.30
9370 C11orf54 11_53 0 1015.23 0.0000 -1.41
381 MED17 11_53 0 54.71 0.0000 1.61
2463 PANX1 11_53 1 1457.10 0.0042 16.23
258 MRE11 11_53 0 396.86 0.0000 2.05
| 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.002 171.41 9.6e-07 15.32
7970 FAM84B 8_83 0.021 31.81 1.9e-06 -2.47
11702 PCAT1 8_83 0.001 5.15 1.8e-08 0.21
11735 CASC19 8_83 0.001 7.20 3.1e-08 0.67
10794 POU5F1B 8_83 0.001 4.91 1.6e-08 0.12
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 1_11"
genename region_tag susie_pip mu2 PVE z
8350 TMEM51 1_11 0.007 6.82 1.5e-07 -0.62
5402 EFHD2 1_11 0.006 5.11 9.1e-08 -0.13
5398 CELA2A 1_11 0.009 9.38 2.5e-07 1.16
4320 CASP9 1_11 0.007 7.33 1.5e-07 1.00
3043 AGMAT 1_11 0.006 6.09 1.1e-07 0.89
3047 PLEKHM2 1_11 0.006 5.17 9.2e-08 -0.24
11270 UQCRHL 1_11 0.007 6.02 1.3e-07 -0.36
599 SPEN 1_11 0.007 6.15 1.3e-07 -0.28
3050 ZBTB17 1_11 0.008 12.97 3.0e-07 -2.56
9739 CLCNKA 1_11 0.006 4.86 8.3e-08 -0.17
8571 HSPB7 1_11 0.010 17.81 5.2e-07 -3.40
9630 FAM131C 1_11 0.020 22.43 1.3e-06 3.21
5400 EPHA2 1_11 0.011 158.11 5.0e-06 -13.66
5401 ARHGEF19 1_11 0.050 27.93 4.1e-06 -2.85
4319 RSG1 1_11 0.006 86.43 1.5e-06 -9.38
352 FBXO42 1_11 0.012 27.56 9.6e-07 4.59
9800 SPATA21 1_11 0.006 6.40 1.1e-07 -1.45
6519 NECAP2 1_11 0.028 18.02 1.5e-06 1.02
11088 LINC01772 1_11 0.007 13.57 2.6e-07 -2.67
10977 NBPF1 1_11 0.066 25.95 5.0e-06 -2.01
11259 LINC01783 1_11 0.013 10.17 3.8e-07 -0.02
| 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
5089 rs4336844 1_11 1.000 181.35 5.3e-04 14.42
36086 rs189026820 1_79 1.000 124.90 3.6e-04 8.41
36140 rs146423333 1_79 1.000 63.50 1.8e-04 -5.39
63263 rs12239046 1_131 1.000 48.22 1.4e-04 -7.26
73290 rs569546056 2_17 1.000 792.78 2.3e-03 3.68
113678 rs1862069 2_102 1.000 50.45 1.5e-04 8.48
148625 rs2649750 3_28 1.000 34.37 1.0e-04 -6.03
173223 rs9870956 3_77 1.000 49.64 1.4e-04 7.14
182009 rs9817452 3_97 1.000 44.64 1.3e-04 6.81
186502 rs149368105 3_105 1.000 104.47 3.0e-04 -11.00
242750 rs11727676 4_94 1.000 42.68 1.2e-04 6.79
291483 rs163895 5_63 1.000 34.89 1.0e-04 -5.66
398421 rs10277379 7_49 1.000 2444.74 7.1e-03 -4.68
398424 rs761767938 7_49 1.000 3143.00 9.1e-03 -3.42
398432 rs1544459 7_49 1.000 3083.39 9.0e-03 -3.87
428030 rs2428 8_11 1.000 556.65 1.6e-03 8.64
428035 rs758184196 8_11 1.000 601.85 1.7e-03 -2.30
434664 rs2293400 8_23 1.000 46.19 1.3e-04 5.84
443080 rs140753685 8_42 1.000 41.72 1.2e-04 6.62
470820 rs4251692 8_94 1.000 57.18 1.7e-04 -9.06
530518 rs9645500 10_46 1.000 148.56 4.3e-04 12.94
568346 rs17157266 11_34 1.000 47.26 1.4e-04 -7.06
577033 rs144988974 11_52 1.000 36.83 1.1e-04 6.20
583173 rs1176746 11_67 1.000 1060.26 3.1e-03 -2.94
583175 rs2307599 11_67 1.000 1060.48 3.1e-03 -2.75
593719 rs12824533 12_11 1.000 63.58 1.8e-04 8.23
596562 rs66720652 12_15 1.000 48.70 1.4e-04 -6.83
673324 rs72681869 14_20 1.000 54.60 1.6e-04 -9.53
687249 rs1243165 14_48 1.000 49.98 1.5e-04 3.73
697122 rs511338 15_14 1.000 44.51 1.3e-04 7.09
702505 rs2070895 15_26 1.000 52.05 1.5e-04 -7.39
734679 rs11645522 16_45 1.000 63.22 1.8e-04 7.57
755357 rs1801689 17_38 1.000 125.02 3.6e-04 11.78
757136 rs1477066 17_41 1.000 46.40 1.3e-04 6.87
777175 rs62094231 18_31 1.000 67.92 2.0e-04 -8.33
798107 rs814573 19_32 1.000 97.36 2.8e-04 -11.60
842322 rs139050 22_19 1.000 157.19 4.6e-04 -12.57
842323 rs6006585 22_19 1.000 62.60 1.8e-04 -5.48
852593 rs35130213 1_19 1.000 428.28 1.2e-03 2.49
867030 rs140584594 1_67 1.000 73.54 2.1e-04 8.43
962588 rs11601507 11_4 1.000 95.75 2.8e-04 9.71
973496 rs145982925 11_21 1.000 4582.54 1.3e-02 3.27
973497 rs35827570 11_21 1.000 4581.14 1.3e-02 3.23
978227 rs2511241 11_41 1.000 38.70 1.1e-04 -6.84
983276 rs111443113 11_53 1.000 27631.80 8.0e-02 -0.61
54555 rs2642420 1_112 0.999 42.47 1.2e-04 -6.69
78667 rs72800939 2_28 0.999 32.21 9.3e-05 5.73
276940 rs536916238 5_33 0.999 41.77 1.2e-04 -0.43
470843 rs74735293 8_94 0.999 61.36 1.8e-04 8.02
795814 rs2251125 19_24 0.999 37.34 1.1e-04 -6.18
842324 rs11090617 22_19 0.999 591.92 1.7e-03 30.16
983266 rs148050219 11_53 0.999 27841.47 8.1e-02 -15.40
1046029 rs117643180 17_6 0.999 78.80 2.3e-04 9.29
54533 rs337171 1_112 0.998 48.53 1.4e-04 7.55
92467 rs894194 2_55 0.998 39.13 1.1e-04 -6.40
470850 rs111470088 8_94 0.998 44.20 1.3e-04 -6.21
324709 rs115740542 6_20 0.997 30.83 8.9e-05 5.67
397361 rs12539160 7_47 0.997 33.27 9.6e-05 5.46
302353 rs769204262 5_84 0.996 30.00 8.7e-05 5.47
409592 rs10435378 7_70 0.996 45.94 1.3e-04 6.83
556853 rs7481951 11_15 0.996 47.40 1.4e-04 7.20
567798 rs1593480 11_34 0.996 31.84 9.2e-05 -5.53
748382 rs3760511 17_22 0.996 28.29 8.2e-05 -5.16
839989 rs132642 22_14 0.996 163.27 4.7e-04 13.69
186523 rs234043 3_106 0.995 28.81 8.3e-05 -5.30
224679 rs77094191 4_59 0.995 59.26 1.7e-04 -4.92
428010 rs11774455 8_11 0.995 107.99 3.1e-04 7.51
795028 rs889140 19_23 0.995 41.53 1.2e-04 -8.34
830968 rs113617088 21_19 0.995 36.74 1.1e-04 -6.67
49488 rs4951163 1_104 0.993 31.90 9.2e-05 5.08
777468 rs12373325 18_31 0.992 165.97 4.8e-04 -14.96
39689 rs9425587 1_84 0.991 29.73 8.6e-05 -5.39
464180 rs10956254 8_83 0.991 31.40 9.0e-05 7.08
795020 rs12610925 19_23 0.991 49.57 1.4e-04 -8.65
275983 rs28499105 5_31 0.990 31.34 9.0e-05 5.56
242958 rs4552481 4_95 0.989 240.07 6.9e-04 16.59
428774 rs11777976 8_13 0.988 57.58 1.7e-04 -8.60
814763 rs6129631 20_24 0.988 45.35 1.3e-04 -5.34
16654 rs12140153 1_39 0.987 27.22 7.8e-05 -5.46
300004 rs72791146 5_79 0.987 29.38 8.4e-05 5.35
816668 rs6066141 20_29 0.987 32.13 9.2e-05 -5.84
325846 rs2235698 6_23 0.986 30.22 8.7e-05 5.96
326664 rs28780090 6_24 0.986 37.52 1.1e-04 6.04
312296 rs9313604 5_103 0.985 28.90 8.3e-05 5.38
661729 rs7318424 13_59 0.985 34.38 9.8e-05 -5.90
471400 rs1016565 9_1 0.984 25.66 7.3e-05 4.91
1037667 rs4782568 16_48 0.983 93.19 2.7e-04 -9.44
679209 rs12894822 14_32 0.981 42.66 1.2e-04 6.62
728197 rs72784008 16_31 0.980 38.45 1.1e-04 -6.24
464139 rs13252684 8_83 0.979 52.65 1.5e-04 3.57
673407 rs2883893 14_20 0.979 29.22 8.3e-05 -5.00
847948 rs114165349 1_18 0.972 120.81 3.4e-04 11.61
790733 rs11668601 19_14 0.967 84.32 2.4e-04 9.97
493528 rs113609637 9_47 0.966 30.79 8.6e-05 -5.79
504741 rs199755552 9_67 0.966 26.02 7.3e-05 -5.18
77708 rs11124740 2_26 0.963 29.20 8.2e-05 -5.31
734678 rs13334801 16_45 0.963 33.26 9.3e-05 4.78
596918 rs67981690 12_16 0.962 34.43 9.6e-05 5.54
550678 rs7939634 11_2 0.960 27.04 7.5e-05 -5.03
484157 rs11557154 9_26 0.958 33.04 9.2e-05 6.11
327877 rs3020583 6_25 0.957 57.47 1.6e-04 9.26
823108 rs2823025 21_2 0.957 28.85 8.0e-05 5.29
809978 rs11557577 20_13 0.956 26.92 7.5e-05 -4.99
188624 rs13089089 3_110 0.950 24.18 6.7e-05 -4.47
815919 rs4810422 20_28 0.948 29.67 8.2e-05 5.40
276922 rs173964 5_33 0.942 105.62 2.9e-04 8.09
536934 rs2243548 10_56 0.939 32.49 8.9e-05 -5.74
586414 rs11220136 11_77 0.933 36.20 9.8e-05 6.02
304888 rs35341726 5_88 0.930 24.83 6.7e-05 4.76
1070201 rs58542926 19_15 0.926 267.31 7.2e-04 19.31
687245 rs941594 14_48 0.922 55.63 1.5e-04 4.99
155602 rs62247577 3_43 0.920 23.51 6.3e-05 4.37
697114 rs17723097 15_13 0.918 44.51 1.2e-04 7.20
786890 rs627379 19_4 0.915 27.78 7.4e-05 -5.10
353667 rs7752846 6_75 0.914 24.33 6.5e-05 -4.77
779149 rs73963711 18_35 0.914 24.78 6.6e-05 -5.02
139758 rs115532219 3_9 0.913 25.44 6.8e-05 3.99
173243 rs141809192 3_78 0.903 24.03 6.3e-05 4.15
571504 rs11236797 11_42 0.901 27.95 7.3e-05 -4.85
269960 rs13172112 5_21 0.899 46.08 1.2e-04 9.04
74618 rs72787520 2_20 0.890 23.81 6.2e-05 -4.53
329826 rs1126511 6_27 0.889 40.81 1.1e-04 7.43
985349 rs72963839 11_54 0.889 40.02 1.0e-04 -6.45
607598 rs7397189 12_36 0.884 22.49 5.8e-05 4.05
139576 rs13085211 3_9 0.882 149.80 3.8e-04 -10.80
30482 rs325937 1_69 0.881 25.89 6.6e-05 -4.74
281126 rs3010273 5_43 0.875 35.72 9.1e-05 -6.07
73083 rs937813 2_16 0.873 27.96 7.1e-05 5.15
428051 rs13265731 8_11 0.867 545.78 1.4e-03 8.08
673372 rs142004400 14_20 0.867 48.06 1.2e-04 -8.89
830896 rs2836882 21_18 0.866 24.11 6.1e-05 4.48
151163 rs7650241 3_35 0.861 26.21 6.6e-05 -4.98
300160 rs6894249 5_79 0.859 28.56 7.1e-05 -5.04
47499 rs12144388 1_99 0.848 27.78 6.8e-05 4.99
464128 rs79658059 8_83 0.848 174.91 4.3e-04 -15.98
673839 rs17252546 14_21 0.844 24.43 6.0e-05 -4.43
473010 rs6915 9_5 0.842 26.26 6.4e-05 -4.88
174605 rs4073154 3_80 0.840 24.39 6.0e-05 -4.44
10529 rs260970 1_24 0.838 24.41 5.9e-05 4.55
690522 rs17617994 14_54 0.835 31.12 7.6e-05 -5.37
73289 rs7562170 2_17 0.832 762.73 1.8e-03 3.52
409 rs10910028 1_2 0.830 42.16 1.0e-04 7.16
635568 rs1778790 13_7 0.829 57.12 1.4e-04 -7.85
833939 rs181391 22_2 0.824 28.32 6.8e-05 5.13
196014 rs13116176 4_4 0.823 30.49 7.3e-05 -6.04
685612 rs243215 14_45 0.820 26.12 6.2e-05 4.77
709725 rs35408448 15_41 0.819 71.26 1.7e-04 -8.67
180512 rs7610095 3_94 0.816 36.17 8.6e-05 -6.13
837883 rs8138401 22_10 0.815 32.10 7.6e-05 5.91
99654 rs192728998 2_70 0.812 24.53 5.8e-05 -4.50
774137 rs9953845 18_26 0.808 30.75 7.2e-05 5.33
73293 rs4580350 2_17 0.807 763.86 1.8e-03 -3.46
737349 rs75079463 16_51 0.806 25.95 6.1e-05 4.61SNPs 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
983266 rs148050219 11_53 0.999 27841.47 8.1e-02 -15.40
983275 rs60550219 11_53 0.077 27806.06 6.2e-03 -15.42
983262 rs7105405 11_53 0.086 27799.88 7.0e-03 -15.44
983300 rs67167563 11_53 0.016 27795.41 1.3e-03 -15.41
983307 rs113426210 11_53 0.114 27786.48 9.2e-03 -15.46
983258 rs9888156 11_53 0.000 27757.14 4.0e-07 -15.40
983313 rs950878 11_53 0.000 27750.73 1.4e-05 -15.44
983256 rs67232024 11_53 0.000 27707.74 6.0e-12 -15.34
983235 rs7927828 11_53 0.000 27703.37 5.9e-14 -15.29
983253 rs9888266 11_53 0.000 27666.45 1.2e-15 -15.31
983241 rs67812366 11_53 0.000 27665.84 7.7e-16 -15.31
983244 rs7109132 11_53 0.000 27665.65 5.9e-16 -15.31
983257 rs16919533 11_53 0.000 27657.35 5.5e-15 -15.37
983236 rs57856352 11_53 0.000 27655.63 8.9e-18 -15.27
983276 rs111443113 11_53 1.000 27631.80 8.0e-02 -0.61
983255 rs67549397 11_53 0.000 27616.41 0.0e+00 -15.29
983254 rs9888143 11_53 0.000 27571.40 0.0e+00 -15.27
983245 rs60351354 11_53 0.000 27568.33 0.0e+00 -15.27
983246 rs60546087 11_53 0.000 27568.29 0.0e+00 -15.27
983250 rs1573567 11_53 0.000 27568.20 0.0e+00 -15.27
983247 rs7109819 11_53 0.000 27568.14 0.0e+00 -15.27
983213 rs7932290 11_53 0.000 27437.31 0.0e+00 -15.28
983181 rs7934467 11_53 0.000 27206.09 0.0e+00 -15.03
983576 rs72966603 11_53 0.000 22977.96 0.0e+00 -16.38
983706 rs12419615 11_53 0.000 21762.21 0.0e+00 -16.55
983757 rs58964858 11_53 0.000 18641.51 0.0e+00 -16.59
983759 rs72968738 11_53 0.000 18604.92 0.0e+00 -16.54
983783 rs138626734 11_53 0.000 18383.34 0.0e+00 -16.62
983801 rs4408267 11_53 0.000 18373.02 0.0e+00 -16.60
983769 rs72968745 11_53 0.000 18362.03 0.0e+00 -16.47
983768 rs4491178 11_53 0.000 18360.68 0.0e+00 -16.46
983829 rs11604580 11_53 0.000 18353.14 0.0e+00 -16.62
983834 rs4342991 11_53 0.000 18351.02 0.0e+00 -16.62
983773 rs4753124 11_53 0.000 18308.87 0.0e+00 -16.63
983806 rs16919942 11_53 0.000 18295.99 0.0e+00 -16.64
983400 rs72962880 11_53 0.000 18212.52 0.0e+00 -12.36
983687 rs7945841 11_53 0.000 18211.46 0.0e+00 -14.94
983389 rs55659547 11_53 0.000 18159.48 0.0e+00 -12.33
983388 rs7950356 11_53 0.000 18156.26 0.0e+00 -12.32
983399 rs56359140 11_53 0.000 18135.77 0.0e+00 -12.29
983392 rs72962872 11_53 0.000 18134.54 0.0e+00 -12.30
983394 rs140989262 11_53 0.000 17898.12 0.0e+00 -12.26
983725 rs7119800 11_53 0.000 17858.28 0.0e+00 -15.03
983414 rs72962891 11_53 0.000 17752.12 0.0e+00 -12.23
983433 rs72964604 11_53 0.000 17715.69 0.0e+00 -12.26
983729 rs2176565 11_53 0.000 17632.70 0.0e+00 -15.17
983730 rs7949551 11_53 0.000 17116.80 0.0e+00 -15.51
983196 rs1506657 11_53 0.000 16789.18 0.0e+00 12.23
983733 rs72968710 11_53 0.000 16731.69 0.0e+00 -15.74
983736 rs16919917 11_53 0.000 16585.55 0.0e+00 -15.60SNPs 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
983266 rs148050219 11_53 0.999 27841.47 0.08100 -15.40
983276 rs111443113 11_53 1.000 27631.80 0.08000 -0.61
973496 rs145982925 11_21 1.000 4582.54 0.01300 3.27
973497 rs35827570 11_21 1.000 4581.14 0.01300 3.23
983307 rs113426210 11_53 0.114 27786.48 0.00920 -15.46
398424 rs761767938 7_49 1.000 3143.00 0.00910 -3.42
398432 rs1544459 7_49 1.000 3083.39 0.00900 -3.87
398421 rs10277379 7_49 1.000 2444.74 0.00710 -4.68
983262 rs7105405 11_53 0.086 27799.88 0.00700 -15.44
398428 rs11972122 7_49 0.773 2879.65 0.00650 -4.21
983275 rs60550219 11_53 0.077 27806.06 0.00620 -15.42
583173 rs1176746 11_67 1.000 1060.26 0.00310 -2.94
583175 rs2307599 11_67 1.000 1060.48 0.00310 -2.75
73290 rs569546056 2_17 1.000 792.78 0.00230 3.68
398429 rs11406602 7_49 0.227 2875.81 0.00190 -4.15
73289 rs7562170 2_17 0.832 762.73 0.00180 3.52
73293 rs4580350 2_17 0.807 763.86 0.00180 -3.46
428035 rs758184196 8_11 1.000 601.85 0.00170 -2.30
842324 rs11090617 22_19 0.999 591.92 0.00170 30.16
428030 rs2428 8_11 1.000 556.65 0.00160 8.64
428051 rs13265731 8_11 0.867 545.78 0.00140 8.08
983300 rs67167563 11_53 0.016 27795.41 0.00130 -15.41
852593 rs35130213 1_19 1.000 428.28 0.00120 2.49
1070201 rs58542926 19_15 0.926 267.31 0.00072 19.31
242958 rs4552481 4_95 0.989 240.07 0.00069 16.59
73292 rs2169748 2_17 0.295 757.07 0.00065 3.44
852599 rs34563832 1_19 0.529 416.25 0.00064 2.81
948714 rs10883451 10_64 0.376 531.89 0.00058 -26.78
948691 rs2862954 10_64 0.360 531.80 0.00056 -26.79
5089 rs4336844 1_11 1.000 181.35 0.00053 14.42
777468 rs12373325 18_31 0.992 165.97 0.00048 -14.96
839989 rs132642 22_14 0.996 163.27 0.00047 13.69
942494 rs7041363 9_59 0.615 261.05 0.00047 -20.99
842322 rs139050 22_19 1.000 157.19 0.00046 -12.57
224684 rs10440365 4_59 0.636 233.98 0.00043 -15.52
464128 rs79658059 8_83 0.848 174.91 0.00043 -15.98
530518 rs9645500 10_46 1.000 148.56 0.00043 12.94
470889 rs7812366 8_94 0.757 191.39 0.00042 -14.09
948692 rs1408579 10_64 0.264 530.93 0.00041 -26.77
139576 rs13085211 3_9 0.882 149.80 0.00038 -10.80
842340 rs9626064 22_19 0.644 203.71 0.00038 14.99
36086 rs189026820 1_79 1.000 124.90 0.00036 8.41
755357 rs1801689 17_38 1.000 125.02 0.00036 11.78
847948 rs114165349 1_18 0.972 120.81 0.00034 11.61
852604 rs2234918 1_19 0.269 414.46 0.00032 2.77
428010 rs11774455 8_11 0.995 107.99 0.00031 7.51
186502 rs149368105 3_105 1.000 104.47 0.00030 -11.00
276922 rs173964 5_33 0.942 105.62 0.00029 8.09
798107 rs814573 19_32 1.000 97.36 0.00028 -11.60
962588 rs11601507 11_4 1.000 95.75 0.00028 9.71SNPs 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
842324 rs11090617 22_19 0.999 591.92 1.7e-03 30.16
842327 rs1977081 22_19 0.001 581.47 1.1e-06 29.82
842334 rs13056555 22_19 0.000 570.82 4.5e-09 29.50
842331 rs2401512 22_19 0.000 569.42 3.3e-09 29.48
842332 rs4823176 22_19 0.000 569.50 3.4e-09 29.48
842333 rs4823178 22_19 0.000 569.59 3.5e-09 29.48
842330 rs2072905 22_19 0.000 569.32 3.3e-09 29.47
842329 rs1883348 22_19 0.000 560.03 8.4e-10 29.26
948691 rs2862954 10_64 0.360 531.80 5.6e-04 -26.79
948714 rs10883451 10_64 0.376 531.89 5.8e-04 -26.78
948692 rs1408579 10_64 0.264 530.93 4.1e-04 -26.77
948773 rs11597086 10_64 0.000 423.32 2.9e-08 -24.56
948808 rs1327578 10_64 0.000 422.70 2.9e-08 -24.55
948832 rs11591741 10_64 0.000 422.72 2.9e-08 -24.55
948844 rs17882431 10_64 0.000 423.06 2.9e-08 -24.55
948848 rs17882802 10_64 0.000 422.70 2.9e-08 -24.55
948850 rs34027394 10_64 0.000 422.47 2.9e-08 -24.54
948857 rs17885645 10_64 0.000 421.47 2.7e-08 -24.52
948887 rs3904036 10_64 0.000 420.98 2.7e-08 -24.51
948893 rs17729876 10_64 0.000 421.09 2.7e-08 -24.51
948881 rs11596950 10_64 0.000 421.04 2.8e-08 -24.50
948895 rs17668255 10_64 0.000 420.26 2.6e-08 -24.50
948909 rs17668357 10_64 0.000 419.79 2.6e-08 -24.49
948968 rs34762508 10_64 0.000 419.44 2.6e-08 -24.48
948975 rs2039015 10_64 0.000 419.22 2.6e-08 -24.47
948734 rs35614792 10_64 0.000 417.93 2.4e-08 -24.45
949019 rs34955138 10_64 0.000 418.48 2.5e-08 -24.45
948978 rs12784396 10_64 0.000 415.90 2.4e-08 -24.40
949015 rs61871342 10_64 0.000 378.85 1.3e-08 -23.52
948999 rs34539738 10_64 0.000 371.06 1.2e-08 -23.32
949006 rs61871341 10_64 0.000 369.97 1.2e-08 -23.31
948669 rs148631880 10_64 0.000 360.69 1.6e-08 -23.08
948996 rs138052038 10_64 0.000 361.83 8.8e-09 -23.05
949000 rs35696979 10_64 0.000 348.51 8.4e-09 -22.25
948997 rs112468457 10_64 0.000 330.14 5.8e-09 -22.20
948955 rs568600628 10_64 0.000 331.61 4.9e-09 -21.98
948670 rs10883448 10_64 0.000 318.14 7.7e-09 -21.96
948679 rs4919426 10_64 0.000 336.10 2.5e-09 -21.94
948666 rs10883447 10_64 0.000 315.25 6.6e-09 -21.91
948700 rs3927496 10_64 0.000 336.18 2.3e-09 -21.89
942494 rs7041363 9_59 0.615 261.05 4.7e-04 -20.99
942491 rs4979373 9_59 0.143 256.98 1.1e-04 -20.91
942485 rs4979372 9_59 0.154 256.89 1.2e-04 -20.85
942501 rs10733608 9_59 0.013 252.34 9.6e-06 -20.83
949035 rs200718402 10_64 0.000 290.92 3.6e-09 -20.74
942483 rs4979371 9_59 0.072 256.25 5.4e-05 -20.72
948624 rs12781812 10_64 0.000 275.40 3.5e-09 -20.60
948617 rs11598495 10_64 0.000 274.30 3.3e-09 -20.58
948973 rs765549407 10_64 0.000 286.06 3.3e-09 -20.56
948622 rs1324693 10_64 0.000 272.13 3.0e-09 -20.53Gene 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] 32if (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 epithelial to mesenchymal transition (GO:0001837)
2 mesenchymal cell differentiation (GO:0048762)
3 positive regulation of production of miRNAs involved in gene silencing by miRNA (GO:1903800)
Overlap Adjusted.P.value Genes
1 3/47 0.02036974 DDX5;AKNA;TGFB1
2 3/51 0.02036974 DDX5;AKNA;TGFB1
3 2/11 0.02496309 DDX5;TGFB1
[1] "GO_Cellular_Component_2021"
Term Overlap Adjusted.P.value Genes
1 microvillus (GO:0005902) 3/57 0.006056438 TGFB1;STARD10;SYTL1
[1] "GO_Molecular_Function_2021"
[1] Term Overlap Adjusted.P.value Genes
<0 rows> (or 0-length row.names)RP11-219B17.3 gene(s) from the input list not found in DisGeNET CURATEDSYTL1 gene(s) from the input list not found in DisGeNET CURATEDPLEKHA3 gene(s) from the input list not found in DisGeNET CURATEDPTPRD-AS1 gene(s) from the input list not found in DisGeNET CURATEDSLF2 gene(s) from the input list not found in DisGeNET CURATEDMLIP gene(s) from the input list not found in DisGeNET CURATEDDLEU1 gene(s) from the input list not found in DisGeNET CURATEDEVA1C gene(s) from the input list not found in DisGeNET CURATEDRP6-109B7.2 gene(s) from the input list not found in DisGeNET CURATEDAKNA gene(s) from the input list not found in DisGeNET CURATEDPTTG1IP gene(s) from the input list not found in DisGeNET CURATED Description FDR Ratio BgRatio
5 Osteoporosis, Age-Related 0.01245970 3/21 61/9703
114 Osteoporosis 0.01245970 3/21 63/9703
116 Osteoporosis, Senile 0.01245970 3/21 61/9703
124 Prostatic Neoplasms 0.01245970 7/21 616/9703
174 Centriacinar Emphysema 0.01245970 2/21 12/9703
191 Panacinar Emphysema 0.01245970 2/21 12/9703
217 Malignant neoplasm of prostate 0.01245970 7/21 616/9703
242 Post-Traumatic Osteoporosis 0.01245970 3/21 61/9703
299 Focal Emphysema 0.01245970 2/21 12/9703
9 Asphyxia 0.01287695 1/21 1/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