Triglycerides (quantile) - Liver
wesleycrouse
2021-09-09
Last updated: 2021-09-09
Checks: 6 1
Knit directory: ctwas_applied/
This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
The R Markdown file has unstaged changes. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20210726) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 59e5f4d. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Unstaged changes:
Modified: analysis/ukb-d-30500_irnt_Liver.Rmd
Modified: analysis/ukb-d-30600_irnt_Liver.Rmd
Modified: analysis/ukb-d-30610_irnt_Liver.Rmd
Modified: analysis/ukb-d-30620_irnt_Liver.Rmd
Modified: analysis/ukb-d-30630_irnt_Liver.Rmd
Modified: analysis/ukb-d-30640_irnt_Liver.Rmd
Modified: analysis/ukb-d-30650_irnt_Liver.Rmd
Modified: analysis/ukb-d-30660_irnt_Liver.Rmd
Modified: analysis/ukb-d-30670_irnt_Liver.Rmd
Modified: analysis/ukb-d-30680_irnt_Liver.Rmd
Modified: analysis/ukb-d-30690_irnt_Liver.Rmd
Modified: analysis/ukb-d-30700_irnt_Liver.Rmd
Modified: analysis/ukb-d-30710_irnt_Liver.Rmd
Modified: analysis/ukb-d-30720_irnt_Liver.Rmd
Modified: analysis/ukb-d-30730_irnt_Liver.Rmd
Modified: analysis/ukb-d-30740_irnt_Liver.Rmd
Modified: analysis/ukb-d-30750_irnt_Liver.Rmd
Modified: analysis/ukb-d-30760_irnt_Liver.Rmd
Modified: analysis/ukb-d-30770_irnt_Liver.Rmd
Modified: analysis/ukb-d-30780_irnt_Liver.Rmd
Modified: analysis/ukb-d-30790_irnt_Liver.Rmd
Modified: analysis/ukb-d-30800_irnt_Liver.Rmd
Modified: analysis/ukb-d-30810_irnt_Liver.Rmd
Modified: analysis/ukb-d-30820_irnt_Liver.Rmd
Modified: analysis/ukb-d-30830_irnt_Liver.Rmd
Modified: analysis/ukb-d-30840_irnt_Liver.Rmd
Modified: analysis/ukb-d-30850_irnt_Liver.Rmd
Modified: analysis/ukb-d-30860_irnt_Liver.Rmd
Modified: analysis/ukb-d-30870_irnt_Liver.Rmd
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-30870_irnt_Liver.Rmd) and HTML (docs/ukb-d-30870_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 | dfd2b5f | wesleycrouse | 2021-09-07 | regenerating reports |
| html | dfd2b5f | wesleycrouse | 2021-09-07 | regenerating reports |
| Rmd | 61b53b3 | wesleycrouse | 2021-09-06 | updated PVE calculation |
| html | 61b53b3 | wesleycrouse | 2021-09-06 | updated PVE calculation |
| Rmd | 837dd01 | wesleycrouse | 2021-09-01 | adding additional fixedsigma report |
| Rmd | d0a5417 | wesleycrouse | 2021-08-30 | adding new reports to the index |
| Rmd | 0922de7 | wesleycrouse | 2021-08-18 | updating all reports with locus plots |
| html | 0922de7 | wesleycrouse | 2021-08-18 | updating all reports with locus plots |
| html | 1c62980 | wesleycrouse | 2021-08-11 | Updating reports |
| Rmd | 5981e80 | wesleycrouse | 2021-08-11 | Adding more reports |
| html | 5981e80 | wesleycrouse | 2021-08-11 | Adding more reports |
| Rmd | da9f015 | wesleycrouse | 2021-08-07 | adding more reports |
| html | da9f015 | wesleycrouse | 2021-08-07 | adding more reports |
Overview
These are the results of a ctwas analysis of the UK Biobank trait Triglycerides (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-30870_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.0072019892 0.0002214704 #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
81.22873 24.09890 #report sample size
print(sample_size)[1] 343992#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.01853874 0.13494305 #compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))[1] 0.09197643 1.10574326Genes 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
6957 USP1 1_39 1.000 1068.51 3.1e-03 33.23
5389 RPS11 19_34 1.000 23812.81 6.9e-02 -3.44
6391 TTC39B 9_13 0.997 33.57 9.7e-05 -5.47
3212 CCND2 12_4 0.988 109.25 3.1e-04 -10.19
12018 RP11-589P10.5 17_6 0.973 35.30 1.0e-04 -5.31
9052 RMI1 9_41 0.967 76.62 2.2e-04 8.88
6509 NTAN1 16_15 0.962 98.95 2.8e-04 12.14
9390 GAS6 13_62 0.961 111.80 3.1e-04 -11.22
8502 RELA 11_36 0.927 33.57 9.0e-05 -5.28
8803 DLEU1 13_21 0.891 39.39 1.0e-04 6.27
6100 ALLC 2_2 0.874 32.03 8.1e-05 -5.31
3330 SEC16B 1_87 0.873 23.30 5.9e-05 4.47
4395 MICAL2 11_9 0.869 28.04 7.1e-05 -4.84
666 COASY 17_25 0.867 69.31 1.7e-04 7.62
6223 GPR180 13_47 0.807 37.13 8.7e-05 5.86
7794 TMC4 19_37 0.794 23.12 5.3e-05 -4.05
3210 LDAH 2_12 0.781 71.21 1.6e-04 8.00
11210 RPS18 6_28 0.772 27.37 6.1e-05 4.56
8271 INSR 19_7 0.766 30.10 6.7e-05 -4.55
2831 ABTB1 3_79 0.750 37.18 8.1e-05 5.89Genes 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 23812.81 6.9e-02 -3.44
1227 FLT3LG 19_34 0.000 20570.01 1.7e-11 2.72
4556 TMEM60 7_49 0.000 9046.08 0.0e+00 3.14
5393 RCN3 19_34 0.000 7743.94 1.3e-07 4.83
1931 FCGRT 19_34 0.000 7062.72 1.7e-10 3.44
2465 APOA5 11_70 0.000 4023.02 0.0e+00 -55.84
3804 PRRG2 19_34 0.432 3434.97 4.3e-03 5.73
4868 BUD13 11_70 0.000 2617.89 0.0e+00 18.89
6005 SIDT2 11_70 0.000 2433.05 0.0e+00 9.05
3803 PRMT1 19_34 0.000 2371.36 7.3e-09 4.76
3805 SCAF1 19_34 0.000 2324.35 5.8e-11 2.67
839 ZNF37A 10_28 0.000 2299.65 0.0e+00 -0.13
3802 IRF3 19_34 0.000 2260.29 1.2e-10 2.86
1940 SLC17A7 19_34 0.000 1708.06 1.2e-11 -3.62
10903 APTR 7_49 0.000 1641.46 0.0e+00 1.75
6957 USP1 1_39 1.000 1068.51 3.1e-03 33.23
6006 TAGLN 11_70 0.000 1058.31 0.0e+00 -18.77
4317 ANGPTL3 1_39 0.006 1053.25 1.7e-05 32.99
9811 RSBN1L 7_49 0.000 982.83 0.0e+00 1.58
2887 NRBP1 2_16 0.000 816.87 4.0e-07 28.09Genes 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 23812.81 6.9e-02 -3.44
3804 PRRG2 19_34 0.432 3434.97 4.3e-03 5.73
6957 USP1 1_39 1.000 1068.51 3.1e-03 33.23
7656 CATSPER2 15_16 0.662 366.08 7.0e-04 19.03
3212 CCND2 12_4 0.988 109.25 3.1e-04 -10.19
9390 GAS6 13_62 0.961 111.80 3.1e-04 -11.22
6509 NTAN1 16_15 0.962 98.95 2.8e-04 12.14
9052 RMI1 9_41 0.967 76.62 2.2e-04 8.88
1058 GCKR 2_16 0.260 268.62 2.0e-04 -20.17
10987 C2orf16 2_16 0.260 268.62 2.0e-04 -20.17
666 COASY 17_25 0.867 69.31 1.7e-04 7.62
3210 LDAH 2_12 0.781 71.21 1.6e-04 8.00
4925 IFT172 2_16 0.519 99.09 1.5e-04 12.84
10879 TM6SF2 19_15 0.490 95.50 1.4e-04 12.71
11521 GSTA2 6_39 0.612 73.54 1.3e-04 8.48
1231 PABPC4 1_24 0.457 77.20 1.0e-04 -8.72
8803 DLEU1 13_21 0.891 39.39 1.0e-04 6.27
12018 RP11-589P10.5 17_6 0.973 35.30 1.0e-04 -5.31
6391 TTC39B 9_13 0.997 33.57 9.7e-05 -5.47
10606 FKBPL 6_26 0.614 53.31 9.5e-05 6.70Genes 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
2465 APOA5 11_70 0.000 4023.02 0.0e+00 -55.84
6957 USP1 1_39 1.000 1068.51 3.1e-03 33.23
4317 ANGPTL3 1_39 0.006 1053.25 1.7e-05 32.99
11684 RP11-136O12.2 8_83 0.000 767.66 3.1e-11 32.67
2887 NRBP1 2_16 0.000 816.87 4.0e-07 28.09
9720 BACE1 11_70 0.000 771.54 0.0e+00 -23.42
2891 SNX17 2_16 0.013 503.44 1.9e-05 -21.66
8284 RBKS 2_16 0.001 406.40 6.2e-07 21.61
5991 FADS1 11_34 0.006 312.60 5.5e-06 -20.55
1058 GCKR 2_16 0.260 268.62 2.0e-04 -20.17
10987 C2orf16 2_16 0.260 268.62 2.0e-04 -20.17
7656 CATSPER2 15_16 0.662 366.08 7.0e-04 19.03
4868 BUD13 11_70 0.000 2617.89 0.0e+00 18.89
6006 TAGLN 11_70 0.000 1058.31 0.0e+00 -18.77
4507 FADS2 11_34 0.001 283.16 7.4e-07 -18.68
7955 FEN1 11_34 0.001 283.16 7.4e-07 -18.68
8739 LPL 8_21 0.000 379.82 0.0e+00 -18.17
1597 PLTP 20_28 0.009 312.46 8.1e-06 -18.02
4936 SLC5A6 2_16 0.000 310.86 2.0e-07 17.19
4941 ATRAID 2_16 0.000 299.85 1.7e-07 -16.94Comparing 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.02999725#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
2465 APOA5 11_70 0.000 4023.02 0.0e+00 -55.84
6957 USP1 1_39 1.000 1068.51 3.1e-03 33.23
4317 ANGPTL3 1_39 0.006 1053.25 1.7e-05 32.99
11684 RP11-136O12.2 8_83 0.000 767.66 3.1e-11 32.67
2887 NRBP1 2_16 0.000 816.87 4.0e-07 28.09
9720 BACE1 11_70 0.000 771.54 0.0e+00 -23.42
2891 SNX17 2_16 0.013 503.44 1.9e-05 -21.66
8284 RBKS 2_16 0.001 406.40 6.2e-07 21.61
5991 FADS1 11_34 0.006 312.60 5.5e-06 -20.55
1058 GCKR 2_16 0.260 268.62 2.0e-04 -20.17
10987 C2orf16 2_16 0.260 268.62 2.0e-04 -20.17
7656 CATSPER2 15_16 0.662 366.08 7.0e-04 19.03
4868 BUD13 11_70 0.000 2617.89 0.0e+00 18.89
6006 TAGLN 11_70 0.000 1058.31 0.0e+00 -18.77
4507 FADS2 11_34 0.001 283.16 7.4e-07 -18.68
7955 FEN1 11_34 0.001 283.16 7.4e-07 -18.68
8739 LPL 8_21 0.000 379.82 0.0e+00 -18.17
1597 PLTP 20_28 0.009 312.46 8.1e-06 -18.02
4936 SLC5A6 2_16 0.000 310.86 2.0e-07 17.19
4941 ATRAID 2_16 0.000 299.85 1.7e-07 -16.94Locus 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: 11_70"
genename region_tag susie_pip mu2 PVE z
4868 BUD13 11_70 0 2617.89 0 18.89
2465 APOA5 11_70 0 4023.02 0 -55.84
3154 APOA1 11_70 0 245.30 0 3.43
7898 PAFAH1B2 11_70 0 548.14 0 4.17
6005 SIDT2 11_70 0 2433.05 0 9.05
6006 TAGLN 11_70 0 1058.31 0 -18.77
6785 PCSK7 11_70 0 540.25 0 -4.60
7745 RNF214 11_70 0 322.44 0 0.09
2466 CEP164 11_70 0 548.73 0 -8.03
9720 BACE1 11_70 0 771.54 0 -23.42
4881 FXYD2 11_70 0 18.83 0 0.56
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 1_39"
genename region_tag susie_pip mu2 PVE z
6956 TM2D1 1_39 0.017 17.29 8.4e-07 2.77
4316 KANK4 1_39 0.009 10.14 2.5e-07 -0.30
6957 USP1 1_39 1.000 1068.51 3.1e-03 33.23
4317 ANGPTL3 1_39 0.006 1053.25 1.7e-05 32.99
3024 DOCK7 1_39 0.010 97.10 2.8e-06 9.51
3733 ATG4C 1_39 0.010 255.93 7.8e-06 -16.11
| 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 767.66 3.1e-11 32.67
7970 FAM84B 8_83 0 5.01 3.3e-15 -0.09
11702 PCAT1 8_83 0 5.44 3.7e-15 0.26
11735 CASC19 8_83 0 6.81 5.5e-15 0.45
10794 POU5F1B 8_83 0 5.24 3.5e-15 0.22
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 2_16"
genename region_tag susie_pip mu2 PVE z
2881 CENPA 2_16 0.000 11.22 7.2e-09 1.66
11149 OST4 2_16 0.000 6.82 3.0e-09 1.15
4939 EMILIN1 2_16 0.000 41.89 1.7e-08 -8.28
4927 KHK 2_16 0.001 13.42 2.3e-08 0.26
4935 PREB 2_16 0.000 86.22 3.8e-08 8.44
4941 ATRAID 2_16 0.000 299.85 1.7e-07 -16.94
4936 SLC5A6 2_16 0.000 310.86 2.0e-07 17.19
1060 CAD 2_16 0.000 175.13 6.5e-08 10.79
2885 SLC30A3 2_16 0.000 166.42 6.1e-08 13.06
7169 UCN 2_16 0.001 62.16 1.2e-07 10.37
2891 SNX17 2_16 0.013 503.44 1.9e-05 -21.66
7170 ZNF513 2_16 0.000 141.59 5.5e-08 8.93
2887 NRBP1 2_16 0.000 816.87 4.0e-07 28.09
4925 IFT172 2_16 0.519 99.09 1.5e-04 12.84
1058 GCKR 2_16 0.260 268.62 2.0e-04 -20.17
10987 C2orf16 2_16 0.260 268.62 2.0e-04 -20.17
10407 GPN1 2_16 0.000 200.03 8.7e-08 11.18
8847 CCDC121 2_16 0.000 31.06 1.5e-08 -4.53
6575 BRE 2_16 0.003 82.24 6.8e-07 -11.11
8284 RBKS 2_16 0.001 406.40 6.2e-07 21.61
| 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.001 7.20 1.8e-08 -0.95
7662 FAM111A 11_34 0.005 22.34 3.4e-07 -1.70
2444 DTX4 11_34 0.010 28.82 8.6e-07 -2.14
10267 MPEG1 11_34 0.004 21.36 2.5e-07 -2.04
7684 PATL1 11_34 0.002 14.38 6.5e-08 -1.96
7687 STX3 11_34 0.001 5.45 1.2e-08 -0.19
7688 MRPL16 11_34 0.001 7.73 2.0e-08 1.26
5997 MS4A2 11_34 0.087 28.83 7.3e-06 -4.83
2453 MS4A6A 11_34 0.063 27.23 5.0e-06 4.68
10924 MS4A4E 11_34 0.001 5.44 1.1e-08 0.55
7698 MS4A14 11_34 0.001 12.49 4.6e-08 -1.92
7697 MS4A7 11_34 0.001 8.93 2.6e-08 1.24
2455 CCDC86 11_34 0.001 11.47 4.6e-08 -1.19
2456 PRPF19 11_34 0.001 6.00 1.4e-08 0.23
2457 TMEM109 11_34 0.001 7.83 2.3e-08 0.65
2480 SLC15A3 11_34 0.001 5.81 1.2e-08 -0.95
2481 CD5 11_34 0.009 31.12 8.0e-07 -3.58
7874 VPS37C 11_34 0.001 7.37 2.2e-08 -0.14
7875 VWCE 11_34 0.001 8.66 2.9e-08 -0.39
5990 TMEM138 11_34 0.004 18.30 2.3e-07 -0.38
6902 CYB561A3 11_34 0.004 18.30 2.3e-07 -0.38
9789 TMEM216 11_34 0.004 24.31 2.9e-07 2.68
11817 RP11-286N22.8 11_34 0.002 13.60 6.9e-08 1.61
5996 CPSF7 11_34 0.001 9.51 3.2e-08 1.45
6903 PPP1R32 11_34 0.001 9.67 4.0e-08 0.14
11812 RP11-794G24.1 11_34 0.001 7.32 2.0e-08 0.94
4508 TMEM258 11_34 0.001 60.73 1.5e-07 7.94
4507 FADS2 11_34 0.001 283.16 7.4e-07 -18.68
7955 FEN1 11_34 0.001 283.16 7.4e-07 -18.68
5991 FADS1 11_34 0.006 312.60 5.5e-06 -20.55
1196 GANAB 11_34 0.001 155.04 3.3e-07 11.63
11004 FADS3 11_34 0.002 12.10 6.9e-08 -2.59
7876 BEST1 11_34 0.001 26.05 5.7e-08 5.66
3676 DKFZP434K028 11_34 0.001 12.18 3.1e-08 -2.44
5994 INCENP 11_34 0.001 15.10 5.3e-08 2.91
6904 ASRGL1 11_34 0.004 23.77 2.6e-07 2.92
| 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
8327 rs79598313 1_18 1.000 117.69 3.4e-04 10.90
57316 rs878811 1_116 1.000 36.77 1.1e-04 -5.84
57893 rs11122453 1_117 1.000 295.33 8.6e-04 -19.86
70999 rs1042034 2_13 1.000 731.88 2.1e-03 26.49
71000 rs1801699 2_13 1.000 49.23 1.4e-04 -4.43
183004 rs9817452 3_97 1.000 47.73 1.4e-04 -6.86
197598 rs3748034 4_4 1.000 65.13 1.9e-04 10.20
197599 rs3752442 4_4 1.000 50.98 1.5e-04 -8.94
225873 rs7439032 4_58 1.000 37.99 1.1e-04 -6.59
231033 rs35518360 4_67 1.000 52.51 1.5e-04 7.21
286144 rs115912456 5_49 1.000 39.34 1.1e-04 -6.22
297460 rs11064 5_72 1.000 66.25 1.9e-04 -8.31
363030 rs540973884 6_92 1.000 146.06 4.2e-04 -12.42
366968 rs746369662 6_99 1.000 238.98 6.9e-04 -4.33
369796 rs12208357 6_103 1.000 90.52 2.6e-04 9.64
398388 rs13235543 7_47 1.000 722.58 2.1e-03 -37.24
398389 rs12539160 7_47 1.000 90.09 2.6e-04 -5.06
398475 rs199603307 7_48 1.000 87.85 2.6e-04 -9.22
399452 rs761767938 7_49 1.000 9227.51 2.7e-02 -2.70
399460 rs1544459 7_49 1.000 8967.57 2.6e-02 -3.09
433835 rs1495743 8_20 1.000 174.39 5.1e-04 -13.73
434550 rs78963197 8_21 1.000 1378.71 4.0e-03 -50.50
434582 rs75835816 8_21 1.000 483.03 1.4e-03 23.18
434618 rs11986461 8_21 1.000 373.72 1.1e-03 -23.75
445504 rs4738679 8_45 1.000 96.45 2.8e-04 -10.23
465170 rs4604455 8_83 1.000 95.79 2.8e-04 16.84
465208 rs10956254 8_83 1.000 75.73 2.2e-04 12.55
521941 rs71007692 10_28 1.000 6923.17 2.0e-02 2.46
585273 rs1176746 11_67 1.000 3490.14 1.0e-02 -3.60
585275 rs2307599 11_67 1.000 3435.43 1.0e-02 -3.57
586066 rs7927993 11_69 1.000 36.23 1.1e-04 5.98
609638 rs7397189 12_36 1.000 81.54 2.4e-04 -10.46
637860 rs1340819 13_7 1.000 35.26 1.0e-04 -5.52
645501 rs7999449 13_25 1.000 23204.93 6.7e-02 3.35
645503 rs775834524 13_25 1.000 23157.19 6.7e-02 3.46
751402 rs55764662 17_26 1.000 190.88 5.5e-04 13.26
757145 rs1801689 17_38 1.000 91.45 2.7e-04 -9.37
791019 rs111500536 19_8 1.000 91.86 2.7e-04 -8.89
791022 rs116843064 19_8 1.000 705.66 2.1e-03 -27.25
792008 rs1865063 19_10 1.000 141.36 4.1e-04 -4.14
794759 rs3794991 19_15 1.000 393.76 1.1e-03 -21.01
794790 rs113619686 19_15 1.000 67.99 2.0e-04 1.58
799632 rs4806075 19_24 1.000 72.70 2.1e-04 -3.80
802464 rs440446 19_32 1.000 746.08 2.2e-03 24.46
802469 rs113345881 19_32 1.000 221.87 6.4e-04 22.31
820225 rs6066141 20_29 1.000 46.10 1.3e-04 -7.62
883040 rs1260326 2_16 1.000 1590.35 4.6e-03 -43.76
884947 rs11688682 2_70 1.000 45.88 1.3e-04 -6.15
913964 rs112436252 6_25 1.000 59.15 1.7e-04 6.51
921954 rs140570886 6_104 1.000 114.89 3.3e-04 -12.29
922107 rs56393506 6_104 1.000 171.36 5.0e-04 -8.17
927696 rs12113583 7_23 1.000 724.24 2.1e-03 -3.81
927698 rs68112310 7_23 1.000 788.67 2.3e-03 -0.60
927709 rs1534696 7_23 1.000 450.38 1.3e-03 -8.95
973876 rs1784130 11_70 1.000 38345.45 1.1e-01 -15.79
973919 rs535499798 11_70 1.000 38248.41 1.1e-01 -1.43
974351 rs964184 11_70 1.000 7015.20 2.0e-02 -76.81
974598 rs141469619 11_70 1.000 831.79 2.4e-03 27.51
1037837 rs113176985 19_34 1.000 22849.04 6.6e-02 3.75
1037840 rs374141296 19_34 1.000 22946.76 6.7e-02 3.60
14831 rs213501 1_34 0.999 36.99 1.1e-04 -5.89
54781 rs56056346 1_111 0.999 91.75 2.7e-04 -4.78
361290 rs212776 6_88 0.999 41.67 1.2e-04 6.61
434548 rs17091881 8_21 0.999 575.84 1.7e-03 22.21
586062 rs28480969 11_69 0.999 37.57 1.1e-04 6.58
595759 rs12824533 12_11 0.999 32.90 9.6e-05 5.49
598958 rs67981690 12_16 0.999 95.97 2.8e-04 9.34
630011 rs35658692 12_75 0.999 92.80 2.7e-04 11.46
736382 rs12443634 16_45 0.999 106.93 3.1e-04 -11.82
790552 rs10401485 19_7 0.999 40.19 1.2e-04 6.88
802462 rs34878901 19_32 0.999 659.43 1.9e-03 -4.29
1035074 rs117380643 17_25 0.999 92.97 2.7e-04 9.72
308799 rs71591078 5_92 0.998 66.70 1.9e-04 0.64
732356 rs57186116 16_38 0.998 35.59 1.0e-04 3.93
818303 rs56206139 20_24 0.998 35.84 1.0e-04 -5.59
823770 rs6099671 20_33 0.998 42.27 1.2e-04 6.56
465167 rs13252684 8_83 0.997 325.72 9.4e-04 9.33
472428 rs1016565 9_1 0.997 30.65 8.9e-05 5.40
528501 rs2393730 10_42 0.997 34.86 1.0e-04 -6.49
598602 rs66720652 12_15 0.996 31.88 9.2e-05 -5.34
375692 rs852386 7_9 0.995 33.03 9.6e-05 5.52
703520 rs10851699 15_28 0.995 38.16 1.1e-04 6.03
741606 rs11078597 17_2 0.995 48.84 1.4e-04 6.83
278519 rs536916238 5_33 0.994 67.51 2.0e-04 3.25
333396 rs115482652 6_34 0.994 30.53 8.8e-05 5.86
333397 rs9472126 6_34 0.994 30.68 8.9e-05 -5.22
434568 rs11986942 8_21 0.994 718.00 2.1e-03 -42.73
550609 rs75184896 10_84 0.994 29.58 8.5e-05 5.17
803719 rs12978750 19_33 0.994 67.95 2.0e-04 9.00
738594 rs7191098 16_50 0.993 31.43 9.1e-05 -5.37
805833 rs12151142 19_38 0.993 47.33 1.4e-04 8.17
949757 rs74563318 10_70 0.993 59.26 1.7e-04 -8.00
555889 rs34623292 11_10 0.992 38.19 1.1e-04 6.41
55247 rs61830291 1_112 0.991 63.07 1.8e-04 7.31
430182 rs2251473 8_14 0.991 65.40 1.9e-04 -10.13
760145 rs1671012 17_42 0.991 32.18 9.3e-05 -5.91
799825 rs2251125 19_24 0.991 29.53 8.5e-05 5.13
497605 rs2254819 9_53 0.990 35.74 1.0e-04 -5.75
630314 rs11057830 12_76 0.990 29.70 8.5e-05 5.41
327050 rs1233385 6_23 0.989 58.57 1.7e-04 -7.67
328284 rs659445 6_26 0.989 71.93 2.1e-04 12.25
606509 rs73108788 12_28 0.989 27.77 8.0e-05 5.07
799058 rs73035067 19_23 0.989 33.55 9.6e-05 6.82
11935 rs1877454 1_26 0.988 31.93 9.2e-05 5.48
398987 rs6465120 7_48 0.988 37.46 1.1e-04 -6.43
721206 rs34340800 16_12 0.988 39.27 1.1e-04 6.18
55238 rs2642420 1_112 0.986 35.22 1.0e-04 4.32
80110 rs4566412 2_31 0.986 33.16 9.5e-05 5.35
97018 rs1821968 2_66 0.986 31.20 8.9e-05 -5.43
533866 rs2186235 10_51 0.986 28.48 8.2e-05 5.14
226045 rs74678260 4_58 0.985 38.27 1.1e-04 -7.72
363039 rs602261 6_93 0.984 27.07 7.7e-05 -3.82
48315 rs6427759 1_99 0.983 27.69 7.9e-05 -4.96
135424 rs4675812 2_144 0.981 31.86 9.1e-05 -5.58
429989 rs7821812 8_14 0.981 79.91 2.3e-04 10.98
328139 rs3132600 6_24 0.980 148.23 4.2e-04 -9.37
440001 rs12675945 8_34 0.978 26.06 7.4e-05 -4.71
808285 rs6139182 20_3 0.978 27.62 7.9e-05 4.86
398125 rs13227753 7_46 0.976 51.43 1.5e-04 -7.26
568311 rs77377156 11_32 0.976 44.66 1.3e-04 6.26
278501 rs173964 5_33 0.974 279.67 7.9e-04 13.39
414574 rs6961342 7_80 0.974 78.12 2.2e-04 12.34
556161 rs547219635 11_11 0.974 31.47 8.9e-05 -4.99
790959 rs117476590 19_7 0.973 31.78 9.0e-05 -5.67
305004 rs76957426 5_85 0.972 38.84 1.1e-04 5.19
732075 rs71401830 16_37 0.972 25.98 7.3e-05 4.39
527719 rs1171619 10_39 0.970 40.89 1.2e-04 6.18
524576 rs71508062 10_33 0.968 26.53 7.5e-05 5.10
1045041 rs8126001 20_38 0.966 50.17 1.4e-04 -7.76
127797 rs62191851 2_129 0.961 25.82 7.2e-05 4.15
226178 rs74939203 4_59 0.959 49.74 1.4e-04 -4.15
329683 rs181268076 6_27 0.958 69.50 1.9e-04 8.71
40775 rs9425587 1_84 0.957 36.18 1.0e-04 -5.74
7503 rs564646712 1_17 0.954 25.40 7.0e-05 4.52
430407 rs1736062 8_15 0.954 50.01 1.4e-04 -9.57
882197 rs777849586 2_16 0.954 32.14 8.9e-05 -5.26
429305 rs330078 8_12 0.953 31.23 8.6e-05 6.55
450582 rs34582181 8_55 0.952 24.92 6.9e-05 4.66
746905 rs139356332 17_16 0.951 24.15 6.7e-05 4.51
703148 rs72748766 15_27 0.950 25.99 7.2e-05 -4.87
798982 rs55897425 19_23 0.950 25.96 7.2e-05 -4.95
584713 rs117978300 11_66 0.948 26.55 7.3e-05 -5.59
369788 rs8191921 6_103 0.945 50.19 1.4e-04 -7.41
110987 rs6722159 2_96 0.944 30.07 8.3e-05 -5.44
545071 rs2420477 10_73 0.944 25.68 7.1e-05 -4.86
16661 rs141797847 1_38 0.939 24.52 6.7e-05 4.51
762836 rs2279396 17_47 0.938 26.13 7.1e-05 -4.72
140724 rs10602803 3_9 0.937 44.50 1.2e-04 4.67
976005 rs555766713 11_70 0.936 502.88 1.4e-03 -21.24
399456 rs11972122 7_49 0.935 9054.64 2.5e-02 -3.36
433574 rs145696392 8_19 0.935 32.93 9.0e-05 5.22
801517 rs6508974 19_28 0.935 31.23 8.5e-05 -4.80
41276 rs375413887 1_85 0.934 23.88 6.5e-05 -4.34
560832 rs56133711 11_19 0.929 39.75 1.1e-04 6.26
604689 rs117339363 12_25 0.928 23.90 6.4e-05 -4.64
742717 rs9904284 17_4 0.928 26.07 7.0e-05 4.06
231099 rs13140033 4_68 0.925 23.95 6.4e-05 4.44
492489 rs12236183 9_45 0.925 34.85 9.4e-05 7.02
769030 rs9963938 18_11 0.920 42.83 1.1e-04 -6.38
697734 rs511338 15_14 0.918 25.81 6.9e-05 4.63
295003 rs2165929 5_67 0.914 29.12 7.7e-05 5.68
969377 rs71468663 11_36 0.911 119.03 3.2e-04 11.42
111985 rs13389219 2_100 0.910 169.93 4.5e-04 -15.59
321204 rs61439239 6_10 0.903 23.29 6.1e-05 4.37
543863 rs17875416 10_71 0.901 25.39 6.7e-05 -4.64
612083 rs189339 12_40 0.900 24.33 6.4e-05 -4.58
7468 rs2742962 1_16 0.897 35.68 9.3e-05 5.73
288201 rs11741599 5_53 0.897 30.06 7.8e-05 5.36
380453 rs38205 7_16 0.893 40.21 1.0e-04 -6.16
428920 rs7826654 8_11 0.893 73.10 1.9e-04 -10.03
277564 rs1694060 5_31 0.888 26.29 6.8e-05 4.47
351633 rs9496567 6_67 0.885 27.99 7.2e-05 -4.86
769331 rs57440424 18_12 0.885 47.92 1.2e-04 -7.19
389508 rs149901303 7_32 0.877 23.40 6.0e-05 4.06
398387 rs13247874 7_47 0.877 745.09 1.9e-03 -37.21
625592 rs73191121 12_66 0.876 29.66 7.5e-05 -4.78
434724 rs74500831 8_22 0.875 35.48 9.0e-05 -5.70
844470 rs9610329 22_14 0.873 38.32 9.7e-05 6.08
820210 rs1412956 20_29 0.872 31.74 8.0e-05 6.56
54768 rs4846295 1_111 0.871 72.85 1.8e-04 1.38
465166 rs2980858 8_83 0.869 855.54 2.2e-03 -34.02
818163 rs932641 20_24 0.867 31.49 7.9e-05 -5.65
310679 rs56034935 5_96 0.866 23.78 6.0e-05 4.30
566558 rs2167079 11_29 0.863 80.93 2.0e-04 -9.71
359403 rs141281941 6_85 0.861 24.49 6.1e-05 4.41
436357 rs145706224 8_24 0.861 24.51 6.1e-05 4.37
197613 rs36205397 4_4 0.860 48.69 1.2e-04 10.94
792009 rs34692794 19_10 0.859 130.07 3.2e-04 1.32
410620 rs10435378 7_70 0.855 38.50 9.6e-05 5.97
538821 rs7074206 10_60 0.855 24.92 6.2e-05 4.44
366959 rs818442 6_99 0.853 220.15 5.5e-04 1.26
586039 rs10047413 11_69 0.852 31.13 7.7e-05 6.72
226065 rs28529445 4_58 0.847 40.02 9.9e-05 8.05
698306 rs148470008 15_15 0.845 75.38 1.9e-04 8.88
736389 rs11641142 16_45 0.844 43.24 1.1e-04 -8.15
736485 rs8061297 16_45 0.842 26.04 6.4e-05 -4.24
188102 rs234043 3_106 0.841 24.44 6.0e-05 4.46
521368 rs7912430 10_27 0.841 24.06 5.9e-05 4.30
682196 rs1848401 14_36 0.839 26.68 6.5e-05 4.99
414540 rs117948936 7_79 0.838 23.70 5.8e-05 -3.03
673936 rs72681869 14_20 0.836 24.43 5.9e-05 -4.32
281631 rs40087 5_41 0.834 24.51 5.9e-05 -4.33
440258 rs72638505 8_34 0.833 24.39 5.9e-05 4.29
332389 rs3734554 6_32 0.831 24.82 6.0e-05 4.42
566529 rs183830453 11_29 0.831 25.60 6.2e-05 -4.13
845256 rs4821764 22_15 0.831 88.68 2.1e-04 9.53
140701 rs4684848 3_9 0.822 114.60 2.7e-04 -8.10
252036 rs13110927 4_109 0.822 24.72 5.9e-05 4.38
329205 rs9273305 6_26 0.821 53.10 1.3e-04 -5.13
329554 rs9276189 6_27 0.821 47.31 1.1e-04 6.86
732603 rs4318234 16_39 0.815 25.66 6.1e-05 4.42
197169 rs1894425 4_3 0.813 23.66 5.6e-05 -4.02
569123 rs11228939 11_32 0.812 30.13 7.1e-05 4.86
429784 rs13249929 8_13 0.809 45.52 1.1e-04 -9.07
743220 rs4968186 17_7 0.808 28.30 6.6e-05 -5.54
197516 rs112820797 4_4 0.804 52.41 1.2e-04 -7.53
206289 rs56147366 4_22 0.804 52.13 1.2e-04 8.62
675556 rs79509284 14_24 0.802 30.92 7.2e-05 4.95SNPs 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
973876 rs1784130 11_70 1.000 38345.45 1.1e-01 -15.79
973919 rs535499798 11_70 1.000 38248.41 1.1e-01 -1.43
973855 rs1784124 11_70 0.000 37821.67 8.8e-13 -15.88
973906 rs12281217 11_70 0.000 37704.90 0.0e+00 -15.89
973953 rs2846012 11_70 0.000 37701.05 0.0e+00 -15.87
973942 rs1784107 11_70 0.000 37701.03 0.0e+00 -15.84
973954 rs1784104 11_70 0.000 37697.25 0.0e+00 -15.86
973875 rs1784129 11_70 0.000 37695.61 0.0e+00 -15.88
973881 rs1784132 11_70 0.000 37658.75 0.0e+00 -15.91
973825 rs1145189 11_70 0.000 36988.73 0.0e+00 -16.19
973820 rs1145187 11_70 0.000 36897.92 0.0e+00 -16.22
973823 rs1145188 11_70 0.000 36885.71 0.0e+00 -16.23
973848 rs1240777 11_70 0.000 33314.65 0.0e+00 -6.83
973922 rs3016355 11_70 0.000 33259.28 0.0e+00 -6.82
973907 rs12294322 11_70 0.000 33257.89 0.0e+00 -6.82
973895 rs2846016 11_70 0.000 33254.23 0.0e+00 -6.82
973893 rs1784133 11_70 0.000 33254.10 0.0e+00 -6.82
973945 rs1787687 11_70 0.000 33247.94 0.0e+00 -6.74
973931 rs1614444 11_70 0.000 33206.53 0.0e+00 -6.91
974047 rs1787700 11_70 0.000 32310.67 0.0e+00 -16.09
974055 rs1942479 11_70 0.000 32229.48 0.0e+00 -16.07
974059 rs111064597 11_70 0.000 32204.49 0.0e+00 -16.05
974077 rs1145192 11_70 0.000 30914.60 0.0e+00 -16.50
974094 rs1145198 11_70 0.000 29916.45 0.0e+00 -16.56
974001 rs758044300 11_70 0.000 29356.68 0.0e+00 -18.17
973836 rs374292052 11_70 0.000 29262.97 0.0e+00 -13.14
973911 rs528101270 11_70 0.000 28874.34 0.0e+00 -5.08
973944 rs561771129 11_70 0.000 27967.12 0.0e+00 -6.49
973874 rs1787682 11_70 0.000 26271.15 0.0e+00 8.41
973901 rs75930640 11_70 0.000 26264.20 0.0e+00 8.37
973905 rs113264811 11_70 0.000 26261.16 0.0e+00 8.36
973886 rs3016354 11_70 0.000 26260.90 0.0e+00 8.37
973887 rs3017337 11_70 0.000 26260.90 0.0e+00 8.37
973880 rs1784131 11_70 0.000 26260.03 0.0e+00 8.37
973884 rs2846015 11_70 0.000 26259.95 0.0e+00 8.37
973862 rs1784125 11_70 0.000 26167.24 0.0e+00 8.45
973870 rs1632755 11_70 0.000 26015.53 0.0e+00 8.33
974000 rs1145210 11_70 0.000 25008.82 0.0e+00 -23.21
974015 rs1787712 11_70 0.000 24859.86 0.0e+00 0.40
645501 rs7999449 13_25 1.000 23204.93 6.7e-02 3.35
645493 rs7337153 13_25 0.018 23182.87 1.2e-03 3.34
645503 rs775834524 13_25 1.000 23157.19 6.7e-02 3.46
645498 rs9537143 13_25 0.000 23105.79 2.1e-05 -3.41
645494 rs9527399 13_25 0.001 23104.73 6.9e-05 -3.43
645497 rs9597193 13_25 0.000 23104.72 3.3e-05 -3.42
645496 rs9527401 13_25 0.000 23104.57 3.2e-05 -3.42
645492 rs9537125 13_25 0.000 23086.77 4.0e-07 -3.39
645491 rs9527398 13_25 0.000 23086.65 3.8e-07 -3.39
645489 rs9537123 13_25 0.000 23085.02 2.9e-07 -3.39
1037840 rs374141296 19_34 1.000 22946.76 6.7e-02 3.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
973876 rs1784130 11_70 1.000 38345.45 0.11000 -15.79
973919 rs535499798 11_70 1.000 38248.41 0.11000 -1.43
645501 rs7999449 13_25 1.000 23204.93 0.06700 3.35
645503 rs775834524 13_25 1.000 23157.19 0.06700 3.46
1037840 rs374141296 19_34 1.000 22946.76 0.06700 3.60
1037837 rs113176985 19_34 1.000 22849.04 0.06600 3.75
399452 rs761767938 7_49 1.000 9227.51 0.02700 -2.70
399460 rs1544459 7_49 1.000 8967.57 0.02600 -3.09
399456 rs11972122 7_49 0.935 9054.64 0.02500 -3.36
521941 rs71007692 10_28 1.000 6923.17 0.02000 2.46
974351 rs964184 11_70 1.000 7015.20 0.02000 -76.81
521940 rs2474565 10_28 0.660 6956.89 0.01300 2.57
585273 rs1176746 11_67 1.000 3490.14 0.01000 -3.60
585275 rs2307599 11_67 1.000 3435.43 0.01000 -3.57
521947 rs2472183 10_28 0.471 6956.72 0.00950 2.56
521950 rs11011452 10_28 0.408 6956.75 0.00830 2.55
883040 rs1260326 2_16 1.000 1590.35 0.00460 -43.76
434550 rs78963197 8_21 1.000 1378.71 0.00400 -50.50
521938 rs9299760 10_28 0.154 6952.46 0.00310 2.56
974598 rs141469619 11_70 1.000 831.79 0.00240 27.51
927698 rs68112310 7_23 1.000 788.67 0.00230 -0.60
465166 rs2980858 8_83 0.869 855.54 0.00220 -34.02
802464 rs440446 19_32 1.000 746.08 0.00220 24.46
70999 rs1042034 2_13 1.000 731.88 0.00210 26.49
398388 rs13235543 7_47 1.000 722.58 0.00210 -37.24
434568 rs11986942 8_21 0.994 718.00 0.00210 -42.73
791022 rs116843064 19_8 1.000 705.66 0.00210 -27.25
927696 rs12113583 7_23 1.000 724.24 0.00210 -3.81
398387 rs13247874 7_47 0.877 745.09 0.00190 -37.21
802462 rs34878901 19_32 0.999 659.43 0.00190 -4.29
399457 rs11406602 7_49 0.065 9036.44 0.00170 -3.30
434548 rs17091881 8_21 0.999 575.84 0.00170 22.21
434582 rs75835816 8_21 1.000 483.03 0.00140 23.18
976005 rs555766713 11_70 0.936 502.88 0.00140 -21.24
927709 rs1534696 7_23 1.000 450.38 0.00130 -8.95
645493 rs7337153 13_25 0.018 23182.87 0.00120 3.34
434618 rs11986461 8_21 1.000 373.72 0.00110 -23.75
794759 rs3794991 19_15 1.000 393.76 0.00110 -21.01
399446 rs12537244 7_49 0.573 629.16 0.00100 -0.05
465167 rs13252684 8_83 0.997 325.72 0.00094 9.33
57893 rs11122453 1_117 1.000 295.33 0.00086 -19.86
1038101 rs552708909 19_34 0.101 2811.46 0.00082 6.08
278501 rs173964 5_33 0.974 279.67 0.00079 13.39
465151 rs2001844 8_83 0.465 543.28 0.00074 -38.69
465153 rs2980886 8_83 0.457 543.24 0.00072 -38.69
366968 rs746369662 6_99 1.000 238.98 0.00069 -4.33
278517 rs3843467 5_33 0.612 357.17 0.00064 15.05
802469 rs113345881 19_32 1.000 221.87 0.00064 22.31
308808 rs12657266 5_92 0.732 268.82 0.00057 14.97
366959 rs818442 6_99 0.853 220.15 0.00055 1.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
974351 rs964184 11_70 1 7015.20 2.0e-02 -76.81
974372 rs3741298 11_70 0 6067.94 0.0e+00 -59.73
974358 rs11604424 11_70 0 6011.05 0.0e+00 -58.13
974377 rs2266788 11_70 0 3928.00 0.0e+00 -57.38
974363 rs7483863 11_70 0 3910.65 0.0e+00 -57.30
974325 rs6589565 11_70 0 3922.71 0.0e+00 -57.24
974319 rs10790162 11_70 0 3922.35 0.0e+00 -57.23
974371 rs10750096 11_70 0 3879.41 0.0e+00 -57.15
974347 rs2160669 11_70 0 3878.76 0.0e+00 -57.13
974365 rs2075290 11_70 0 3797.07 0.0e+00 -56.56
974293 rs3825041 11_70 0 3769.42 0.0e+00 -56.13
974263 rs6589564 11_70 0 3759.37 0.0e+00 -56.04
974381 rs2072560 11_70 0 3770.85 0.0e+00 -55.88
974384 rs651821 11_70 0 3804.52 0.0e+00 -55.87
974387 rs662799 11_70 0 3800.85 0.0e+00 -55.84
974266 rs7930786 11_70 0 3719.06 0.0e+00 -55.83
974220 rs9326246 11_70 0 3686.26 0.0e+00 -54.85
974210 rs1974718 11_70 0 3664.57 0.0e+00 -54.64
974211 rs1558860 11_70 0 3661.80 0.0e+00 -54.63
974212 rs1558861 11_70 0 3663.99 0.0e+00 -54.63
974146 rs7350481 11_70 0 3436.12 0.0e+00 -51.72
434550 rs78963197 8_21 1 1378.71 4.0e-03 -50.50
434564 rs6999813 8_21 0 1391.94 1.4e-07 -50.45
434543 rs10645926 8_21 0 1341.43 0.0e+00 -50.15
974337 rs7118999 11_70 0 5426.84 0.0e+00 -50.03
434555 rs17489226 8_21 0 1312.98 2.2e-15 -49.74
974390 rs10750097 11_70 0 4675.06 0.0e+00 -49.73
434581 rs7816447 8_21 0 1366.16 8.4e-16 -49.70
434583 rs11984698 8_21 0 1353.62 5.6e-17 -49.64
434576 rs28675909 8_21 0 1346.12 1.0e-17 -49.60
434579 rs11989309 8_21 0 1346.68 1.0e-17 -49.60
434578 rs7004149 8_21 0 1342.81 4.3e-18 -49.58
434577 rs79198716 8_21 0 1342.21 3.5e-18 -49.57
974383 rs3135506 11_70 0 2863.68 0.0e+00 48.30
974370 rs35120633 11_70 0 2827.14 0.0e+00 47.99
434545 rs1569209 8_21 0 1227.52 0.0e+00 -47.91
974354 rs61905116 11_70 0 2804.39 0.0e+00 47.69
974382 rs12287066 11_70 0 2824.56 0.0e+00 47.64
974374 rs12285095 11_70 0 2811.02 0.0e+00 47.41
434547 rs80073370 8_21 0 1157.43 0.0e+00 -47.40
974313 rs12294259 11_70 0 2800.12 0.0e+00 47.30
974308 rs60972755 11_70 0 2796.73 0.0e+00 47.29
974373 rs140050044 11_70 0 2792.50 0.0e+00 47.24
974295 rs10488699 11_70 0 2804.70 0.0e+00 47.19
974270 rs56225305 11_70 0 2793.35 0.0e+00 47.14
974279 rs56224630 11_70 0 2793.62 0.0e+00 47.14
974298 rs57641217 11_70 0 2796.36 0.0e+00 47.02
974301 rs11820589 11_70 0 2772.40 0.0e+00 47.01
434559 rs2410620 8_21 0 788.67 1.5e-10 -46.93
434566 rs1441762 8_21 0 786.19 1.3e-10 -46.91Gene set enrichment for genes with PIP>0.8
#GO enrichment analysis
library(enrichR)Welcome to enrichR
Checking connection ... Enrichr ... Connection is Live!
FlyEnrichr ... Connection is available!
WormEnrichr ... Connection is available!
YeastEnrichr ... Connection is available!
FishEnrichr ... Connection is available!dbs <- c("GO_Biological_Process_2021", "GO_Cellular_Component_2021", "GO_Molecular_Function_2021")
genes <- ctwas_gene_res$genename[ctwas_gene_res$susie_pip>0.8]
#number of genes for gene set enrichment
length(genes)[1] 15if (length(genes)>0){
GO_enrichment <- enrichr(genes, dbs)
for (db in dbs){
print(db)
df <- GO_enrichment[[db]]
df <- df[df$Adjusted.P.value<0.05,c("Term", "Overlap", "Adjusted.P.value", "Genes")]
print(df)
}
#DisGeNET enrichment
# devtools::install_bitbucket("ibi_group/disgenet2r")
library(disgenet2r)
disgenet_api_key <- get_disgenet_api_key(
email = "wesleycrouse@gmail.com",
password = "uchicago1" )
Sys.setenv(DISGENET_API_KEY= disgenet_api_key)
res_enrich <-disease_enrichment(entities=genes, vocabulary = "HGNC",
database = "CURATED" )
df <- res_enrich@qresult[1:10, c("Description", "FDR", "Ratio", "BgRatio")]
print(df)
#WebGestalt enrichment
library(WebGestaltR)
background <- ctwas_gene_res$genename
#listGeneSet()
databases <- c("pathway_KEGG", "disease_GLAD4U", "disease_OMIM")
enrichResult <- WebGestaltR(enrichMethod="ORA", organism="hsapiens",
interestGene=genes, referenceGene=background,
enrichDatabase=databases, interestGeneType="genesymbol",
referenceGeneType="genesymbol", isOutput=F)
print(enrichResult[,c("description", "size", "overlap", "FDR", "database", "userId")])
}Uploading data to Enrichr... Done.
Querying GO_Biological_Process_2021... Done.
Querying GO_Cellular_Component_2021... Done.
Querying GO_Molecular_Function_2021... Done.
Parsing results... Done.
[1] "GO_Biological_Process_2021"
[1] Term Overlap Adjusted.P.value Genes
<0 rows> (or 0-length row.names)
[1] "GO_Cellular_Component_2021"
[1] Term Overlap Adjusted.P.value Genes
<0 rows> (or 0-length row.names)
[1] "GO_Molecular_Function_2021"
[1] Term Overlap Adjusted.P.value Genes
<0 rows> (or 0-length row.names)DLEU1 gene(s) from the input list not found in DisGeNET CURATEDTTC39B gene(s) from the input list not found in DisGeNET CURATEDALLC gene(s) from the input list not found in DisGeNET CURATEDRP11-589P10.5 gene(s) from the input list not found in DisGeNET CURATEDUSP1 gene(s) from the input list not found in DisGeNET CURATEDMICAL2 gene(s) from the input list not found in DisGeNET CURATEDRMI1 gene(s) from the input list not found in DisGeNET CURATEDRPS11 gene(s) from the input list not found in DisGeNET CURATEDNTAN1 gene(s) from the input list not found in DisGeNET CURATED Description
28 Diabetic Nephropathy
31 Nodular glomerulosclerosis
105 NEURODEGENERATION WITH BRAIN IRON ACCUMULATION 6
108 MEGALENCEPHALY-POLYMICROGYRIA-POLYDACTYLY-HYDROCEPHALUS SYNDROME 3
111 Coenzyme A synthase protein associated neurodegeneration
114 PONTOCEREBELLAR HYPOPLASIA, TYPE 12
32 Cardiomegaly
91 Cardiac Hypertrophy
101 Alcohol Toxicity
110 RELA fusion-positive ependymoma
FDR Ratio BgRatio
28 0.01175015 2/6 44/9703
31 0.01175015 2/6 41/9703
105 0.01175015 1/6 1/9703
108 0.01175015 1/6 1/9703
111 0.01175015 1/6 1/9703
114 0.01175015 1/6 1/9703
32 0.01409655 2/6 82/9703
91 0.01409655 2/6 82/9703
101 0.01409655 1/6 2/9703
110 0.01409655 1/6 2/9703******************************************
* *
* Welcome to WebGestaltR ! *
* *
******************************************Loading the functional categories...
Loading the ID list...
Loading the reference list...
Performing the enrichment analysis...Warning in oraEnrichment(interestGeneList, referenceGeneList, geneSet,
minNum = minNum, : No significant gene set is identified based on FDR 0.05!NULL
sessionInfo()R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)
Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] WebGestaltR_0.4.4 disgenet2r_0.99.2 enrichR_3.0 cowplot_1.0.0
[5] ggplot2_3.3.3
loaded via a namespace (and not attached):
[1] bitops_1.0-6 matrixStats_0.57.0
[3] fs_1.3.1 bit64_4.0.5
[5] doParallel_1.0.16 progress_1.2.2
[7] httr_1.4.1 rprojroot_2.0.2
[9] GenomeInfoDb_1.20.0 doRNG_1.8.2
[11] tools_3.6.1 utf8_1.2.1
[13] R6_2.5.0 DBI_1.1.1
[15] BiocGenerics_0.30.0 colorspace_1.4-1
[17] withr_2.4.1 tidyselect_1.1.0
[19] prettyunits_1.0.2 bit_4.0.4
[21] curl_3.3 compiler_3.6.1
[23] git2r_0.26.1 Biobase_2.44.0
[25] DelayedArray_0.10.0 rtracklayer_1.44.0
[27] labeling_0.3 scales_1.1.0
[29] readr_1.4.0 apcluster_1.4.8
[31] stringr_1.4.0 digest_0.6.20
[33] Rsamtools_2.0.0 svglite_1.2.2
[35] rmarkdown_1.13 XVector_0.24.0
[37] pkgconfig_2.0.3 htmltools_0.3.6
[39] fastmap_1.1.0 BSgenome_1.52.0
[41] rlang_0.4.11 RSQLite_2.2.7
[43] generics_0.0.2 farver_2.1.0
[45] jsonlite_1.6 BiocParallel_1.18.0
[47] dplyr_1.0.7 VariantAnnotation_1.30.1
[49] RCurl_1.98-1.1 magrittr_2.0.1
[51] GenomeInfoDbData_1.2.1 Matrix_1.2-18
[53] Rcpp_1.0.6 munsell_0.5.0
[55] S4Vectors_0.22.1 fansi_0.5.0
[57] gdtools_0.1.9 lifecycle_1.0.0
[59] stringi_1.4.3 whisker_0.3-2
[61] yaml_2.2.0 SummarizedExperiment_1.14.1
[63] zlibbioc_1.30.0 plyr_1.8.4
[65] grid_3.6.1 blob_1.2.1
[67] parallel_3.6.1 promises_1.0.1
[69] crayon_1.4.1 lattice_0.20-38
[71] Biostrings_2.52.0 GenomicFeatures_1.36.3
[73] hms_1.1.0 knitr_1.23
[75] pillar_1.6.1 igraph_1.2.4.1
[77] GenomicRanges_1.36.0 rjson_0.2.20
[79] rngtools_1.5 codetools_0.2-16
[81] reshape2_1.4.3 biomaRt_2.40.1
[83] stats4_3.6.1 XML_3.98-1.20
[85] glue_1.4.2 evaluate_0.14
[87] data.table_1.14.0 foreach_1.5.1
[89] vctrs_0.3.8 httpuv_1.5.1
[91] gtable_0.3.0 purrr_0.3.4
[93] assertthat_0.2.1 cachem_1.0.5
[95] xfun_0.8 later_0.8.0
[97] tibble_3.1.2 iterators_1.0.13
[99] GenomicAlignments_1.20.1 AnnotationDbi_1.46.0
[101] memoise_2.0.0 IRanges_2.18.1
[103] workflowr_1.6.2 ellipsis_0.3.2