Red blood cell (erythrocyte) count - Whole_Blood
wesleycrouse
2021-09-7
Last updated: 2021-09-13
Checks: 7 0
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.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
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 72c8ef7. 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:
working directory clean
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-30010_irnt_Whole_Blood_known.Rmd) and HTML (docs/ukb-d-30010_irnt_Whole_Blood_known.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 | 72c8ef7 | wesleycrouse | 2021-09-13 | changing mart for biomart |
| Rmd | a49c40e | wesleycrouse | 2021-09-13 | updating with bystander analysis |
| html | 7e22565 | wesleycrouse | 2021-09-13 | updating reports |
| Rmd | 3a7fbc1 | wesleycrouse | 2021-09-08 | generating reports for known annotations |
| html | 3a7fbc1 | wesleycrouse | 2021-09-08 | generating reports for known annotations |
| Rmd | cbf7408 | wesleycrouse | 2021-09-08 | adding enrichment to reports |
Overview
These are the results of a ctwas analysis of the UK Biobank trait Red blood cell (erythrocyte) count using Whole_Blood gene weights.
The GWAS was conducted by the Neale Lab, and the biomarker traits we analyzed are discussed here. Summary statistics were obtained from IEU OpenGWAS using GWAS ID: ukb-d-30010_irnt. Results were obtained from from IEU rather than Neale Lab because they are in a standardard format (GWAS VCF). Note that 3 of the 34 biomarker traits were not available from IEU and were excluded from analysis.
The weights are mashr GTEx v8 models on Whole_Blood eQTL obtained from PredictDB. We performed a full harmonization of the variants, including recovering strand ambiguous variants. This procedure is discussed in a separate document. (TO-DO: add report that describes harmonization)
LD matrices were computed from a 10% subset of Neale lab subjects. Subjects were matched using the plate and well information from genotyping. We included only biallelic variants with MAF>0.01 in the original Neale Lab GWAS. (TO-DO: add more details [number of subjects, variants, etc])
Weight QC
TO-DO: add enhanced QC reporting (total number of weights, why each variant was missing for all genes)
qclist_all <- list()
qc_files <- paste0(results_dir, "/", list.files(results_dir, pattern="exprqc.Rd"))
for (i in 1:length(qc_files)){
load(qc_files[i])
chr <- unlist(strsplit(rev(unlist(strsplit(qc_files[i], "_")))[1], "[.]"))[1]
qclist_all[[chr]] <- cbind(do.call(rbind, lapply(qclist,unlist)), as.numeric(substring(chr,4)))
}
qclist_all <- data.frame(do.call(rbind, qclist_all))
colnames(qclist_all)[ncol(qclist_all)] <- "chr"
rm(qclist, wgtlist, z_gene_chr)
#number of imputed weights
nrow(qclist_all)[1] 11095#number of imputed weights by chromosome
table(qclist_all$chr)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
1129 747 624 400 479 621 560 383 404 430 682 652 192 362 331
16 17 18 19 20 21 22
551 725 159 911 313 130 310 #proportion of imputed weights without missing variants
mean(qclist_all$nmiss==0)[1] 0.762776Load ctwas results
#load ctwas results
ctwas_res <- data.table::fread(paste0(results_dir, "/", analysis_id, "_ctwas.susieIrss.txt"))
#make unique identifier for regions
ctwas_res$region_tag <- paste(ctwas_res$region_tag1, ctwas_res$region_tag2, sep="_")
#load z scores for SNPs and collect sample size
load(paste0(results_dir, "/", analysis_id, "_expr_z_snp.Rd"))
sample_size <- z_snp$ss
sample_size <- as.numeric(names(which.max(table(sample_size))))
#compute PVE for each gene/SNP
ctwas_res$PVE = ctwas_res$susie_pip*ctwas_res$mu2/sample_size
#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 scores to results
load(paste0(results_dir, "/", analysis_id, "_expr_z_gene.Rd"))
ctwas_gene_res$z <- z_gene[ctwas_gene_res$id,]$z
z_snp <- z_snp[z_snp$id %in% ctwas_snp_res$id,]
ctwas_snp_res$z <- z_snp$z[match(ctwas_snp_res$id, z_snp$id)]
#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_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 |
|---|---|---|
| 3a7fbc1 | wesleycrouse | 2021-09-08 |
#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.0167529576 0.0002048743 #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.41100 22.58389 #report sample size
print(sample_size)[1] 350475#report group size
group_size <- c(nrow(ctwas_gene_res), n_snps)
print(group_size)[1] 11095 8696600#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.01559809 0.11480972 #compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))[1] 0.270587 2.989793Genes with highest PIPs
#distribution of PIPs
hist(ctwas_gene_res$susie_pip, xlim=c(0,1), main="Distribution of Gene PIPs")
| Version | Author | Date |
|---|---|---|
| 3a7fbc1 | wesleycrouse | 2021-09-08 |
#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
3193 PIK3R3 1_28 1.000 196.02 5.6e-04 -13.79
3208 DARS2 1_85 1.000 38989.57 1.1e-01 5.12
797 TFRC 3_120 1.000 153.92 4.4e-04 -18.27
12332 U91328.22 6_20 1.000 95.06 2.7e-04 -6.80
11668 CEBPA 19_23 1.000 69.02 2.0e-04 -6.91
1980 FCGRT 19_34 1.000 37858.86 1.1e-01 -5.07
301 TYMP 22_24 1.000 336.08 9.6e-04 -20.84
8670 MTX1 1_77 0.999 82.58 2.4e-04 11.39
5665 CNIH4 1_114 0.999 37.26 1.1e-04 -6.64
2106 KLF1 19_11 0.991 279.50 7.9e-04 -20.21
12494 DUSP14 17_22 0.979 42.22 1.2e-04 6.33
7345 RNF168 3_121 0.976 57.70 1.6e-04 3.42
8999 LPCAT4 15_10 0.973 22.61 6.3e-05 -4.38
12181 EGLN2 19_30 0.967 99.52 2.7e-04 -10.90
6290 ZFP36L2 2_27 0.964 41.27 1.1e-04 6.07
9284 SERTAD2 2_42 0.962 33.97 9.3e-05 -5.67
8523 ZNF217 20_31 0.959 36.96 1.0e-04 5.89
1523 CHKB 22_24 0.957 94.97 2.6e-04 -13.37
4915 TEX10 9_50 0.956 25.62 7.0e-05 4.71
1803 DHODH 16_38 0.955 25.86 7.0e-05 3.99
7460 ANKRD55 5_33 0.949 27.69 7.5e-05 -4.83
11449 SMIM1 1_3 0.947 81.55 2.2e-04 9.06
12180 CTC-490E21.10 19_30 0.945 47.54 1.3e-04 5.39
4283 TRAF3 14_54 0.942 28.58 7.7e-05 -4.91
9567 AFMID 17_43 0.936 37.46 1.0e-04 -5.63
11168 FAM228B 2_14 0.931 54.18 1.4e-04 10.19
84 OSBPL7 17_28 0.923 37.30 9.8e-05 6.72
3569 NT5C3A 7_25 0.919 40.58 1.1e-04 6.43
12615 EXOC3L2 19_31 0.919 37.77 9.9e-05 6.70
8170 NPIPB12 16_24 0.918 67.69 1.8e-04 -11.53
433 VAMP3 1_6 0.912 25.85 6.7e-05 4.80
2686 SERINC1 6_82 0.909 43.61 1.1e-04 -6.51
8816 UGT8 4_73 0.906 47.63 1.2e-04 -7.02
6813 PSKH1 16_36 0.898 107.96 2.8e-04 -10.97
6089 FADS1 11_34 0.897 143.45 3.7e-04 13.04
9692 ACTRT3 3_104 0.883 51.53 1.3e-04 6.74
6255 ARL11 13_21 0.882 22.10 5.6e-05 -3.89
4631 DAGLA 11_34 0.880 33.28 8.4e-05 -7.36
5130 SPPL2A 15_20 0.877 33.62 8.4e-05 -5.82
6590 NTAN1 16_15 0.877 26.07 6.5e-05 5.03
3867 MMP24 20_21 0.873 25.85 6.4e-05 3.99
9830 PBX1 1_80 0.870 28.44 7.1e-05 5.01
8466 LRRC8C 1_55 0.861 21.05 5.2e-05 -4.02
2682 MCM9 6_79 0.858 31.85 7.8e-05 5.38
615 NTHL1 16_2 0.845 37.19 9.0e-05 -5.50
11041 LAT 16_23 0.840 26.62 6.4e-05 -6.03
6253 RNF219 13_39 0.839 21.44 5.1e-05 -4.25
1415 SETD1A 16_24 0.834 146.99 3.5e-04 -10.46
9813 MUC1 1_77 0.825 102.55 2.4e-04 10.52
2612 CDKN1B 12_12 0.820 22.21 5.2e-05 4.33
7506 CDK5 7_94 0.801 68.78 1.6e-04 8.46Genes 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 |
|---|---|---|
| 3a7fbc1 | wesleycrouse | 2021-09-08 |
#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
168 SPRTN 1_118 0.000 60519.03 0.0e+00 -8.13
3138 EXOC8 1_118 0.000 43517.88 0.0e+00 -7.29
3208 DARS2 1_85 1.000 38989.57 1.1e-01 5.12
1980 FCGRT 19_34 1.000 37858.86 1.1e-01 -5.07
9785 ZBTB37 1_85 0.000 26931.89 0.0e+00 6.10
3207 PRDX6 1_85 0.000 16172.75 0.0e+00 -3.63
5520 RCN3 19_34 0.000 12016.60 1.2e-14 -5.07
5911 PURB 7_33 0.312 11350.90 1.0e-02 -12.87
3361 ZNF410 14_34 0.000 10269.34 0.0e+00 3.20
3211 SERPINC1 1_85 0.000 9500.63 0.0e+00 -2.54
6245 RABGAP1L 1_85 0.000 9432.83 0.0e+00 -1.93
905 KLHL20 1_85 0.000 8700.39 0.0e+00 0.59
881 ZNF37A 10_28 0.000 8610.60 0.0e+00 3.25
3140 TSNAX 1_118 0.000 8371.43 0.0e+00 0.62
8678 MFSD4B 6_74 0.000 7618.08 6.5e-09 -5.77
11290 RP11-160H22.5 1_85 0.000 7552.41 0.0e+00 -4.12
2741 SLC16A10 6_74 0.000 7453.44 2.1e-07 4.57
4733 AHI1 6_89 0.000 7026.10 0.0e+00 3.54
7145 DISC1 1_118 0.000 5088.00 0.0e+00 -2.69
5284 PTGR2 14_34 0.000 4594.80 0.0e+00 0.32Genes 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
3208 DARS2 1_85 1.000 38989.57 0.11000 5.12
1980 FCGRT 19_34 1.000 37858.86 0.11000 -5.07
5911 PURB 7_33 0.312 11350.90 0.01000 -12.87
3268 ALDH8A1 6_89 0.724 1380.15 0.00290 13.02
11039 PPP1CB 2_19 0.432 1637.28 0.00200 -5.57
301 TYMP 22_24 1.000 336.08 0.00096 -20.84
8411 TRMT61B 2_19 0.198 1652.35 0.00093 5.28
2106 KLF1 19_11 0.991 279.50 0.00079 -20.21
3193 PIK3R3 1_28 1.000 196.02 0.00056 -13.79
797 TFRC 3_120 1.000 153.92 0.00044 -18.27
6089 FADS1 11_34 0.897 143.45 0.00037 13.04
9223 ZBTB7A 19_4 0.725 174.32 0.00036 13.56
1415 SETD1A 16_24 0.834 146.99 0.00035 -10.46
6813 PSKH1 16_36 0.898 107.96 0.00028 -10.97
9605 FAM46C 1_72 0.665 140.55 0.00027 -6.92
12332 U91328.22 6_20 1.000 95.06 0.00027 -6.80
12181 EGLN2 19_30 0.967 99.52 0.00027 -10.90
1523 CHKB 22_24 0.957 94.97 0.00026 -13.37
8670 MTX1 1_77 0.999 82.58 0.00024 11.39
9813 MUC1 1_77 0.825 102.55 0.00024 10.52Genes 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
2736 HBS1L 6_89 0.000 2366.86 0.0e+00 38.09
12556 PRICKLE4 6_32 0.000 699.62 4.4e-19 -25.00
3734 RP11-298J23.10 6_32 0.000 693.34 4.4e-19 24.93
2751 BYSL 6_32 0.000 679.00 0.0e+00 23.99
8866 ABO 9_70 0.000 275.33 7.0e-08 21.31
301 TYMP 22_24 1.000 336.08 9.6e-04 -20.84
2106 KLF1 19_11 0.991 279.50 7.9e-04 -20.21
9142 ODF3B 22_24 0.004 306.50 3.2e-06 -20.20
2194 MOSPD3 7_62 0.000 159.76 2.2e-14 -19.61
9969 MAPT 17_27 0.000 345.47 1.1e-08 -18.36
797 TFRC 3_120 1.000 153.92 4.4e-04 -18.27
6997 SYCE2 19_11 0.001 214.67 3.6e-07 17.71
9242 FARSA 19_11 0.001 214.65 3.6e-07 -17.64
2195 PCOLCE 7_62 0.000 149.61 1.4e-12 -17.15
4731 CD164 6_73 0.003 195.08 1.8e-06 -16.82
2363 WNT3 17_27 0.000 268.25 4.9e-09 16.67
3742 MED20 6_32 0.000 353.61 1.1e-19 -15.60
11112 TOMM6 6_32 0.000 361.32 1.1e-19 -15.36
9017 LRRC37A 17_27 0.000 161.66 1.8e-07 15.31
6117 TBX6 16_24 0.058 196.36 3.2e-05 -14.78Comparing 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 |
|---|---|---|
| 3a7fbc1 | wesleycrouse | 2021-09-08 |
#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 |
|---|---|---|
| 3a7fbc1 | wesleycrouse | 2021-09-08 |
#proportion of significant z scores
mean(abs(ctwas_gene_res$z) > sig_thresh)[1] 0.04849031#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
2736 HBS1L 6_89 0.000 2366.86 0.0e+00 38.09
12556 PRICKLE4 6_32 0.000 699.62 4.4e-19 -25.00
3734 RP11-298J23.10 6_32 0.000 693.34 4.4e-19 24.93
2751 BYSL 6_32 0.000 679.00 0.0e+00 23.99
8866 ABO 9_70 0.000 275.33 7.0e-08 21.31
301 TYMP 22_24 1.000 336.08 9.6e-04 -20.84
2106 KLF1 19_11 0.991 279.50 7.9e-04 -20.21
9142 ODF3B 22_24 0.004 306.50 3.2e-06 -20.20
2194 MOSPD3 7_62 0.000 159.76 2.2e-14 -19.61
9969 MAPT 17_27 0.000 345.47 1.1e-08 -18.36
797 TFRC 3_120 1.000 153.92 4.4e-04 -18.27
6997 SYCE2 19_11 0.001 214.67 3.6e-07 17.71
9242 FARSA 19_11 0.001 214.65 3.6e-07 -17.64
2195 PCOLCE 7_62 0.000 149.61 1.4e-12 -17.15
4731 CD164 6_73 0.003 195.08 1.8e-06 -16.82
2363 WNT3 17_27 0.000 268.25 4.9e-09 16.67
3742 MED20 6_32 0.000 353.61 1.1e-19 -15.60
11112 TOMM6 6_32 0.000 361.32 1.1e-19 -15.36
9017 LRRC37A 17_27 0.000 161.66 1.8e-07 15.31
6117 TBX6 16_24 0.058 196.36 3.2e-05 -14.78Locus 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: 6_89"
genename region_tag susie_pip mu2 PVE z
3268 ALDH8A1 6_89 0.724 1380.15 0.0029 13.02
2736 HBS1L 6_89 0.000 2366.86 0.0000 38.09
4733 AHI1 6_89 0.000 7026.10 0.0000 3.54
| Version | Author | Date |
|---|---|---|
| 3a7fbc1 | wesleycrouse | 2021-09-08 |
[1] "Region: 6_32"
genename region_tag susie_pip mu2 PVE z
6555 LRFN2 6_32 0.000 24.02 2.3e-20 1.92
3736 UNC5CL 6_32 0.000 7.44 0.0e+00 0.65
2715 TSPO2 6_32 0.000 13.61 4.3e-21 0.99
3746 APOBEC2 6_32 0.000 9.71 3.1e-21 -0.93
7 NFYA 6_32 0.000 22.35 1.4e-20 1.83
1377 TREM2 6_32 0.000 18.10 5.7e-21 -1.28
2712 TREML2 6_32 0.000 8.14 2.6e-21 -0.69
10066 TREML4 6_32 0.000 7.98 0.0e+00 -0.79
3749 TREM1 6_32 0.000 5.46 0.0e+00 0.74
3734 RP11-298J23.10 6_32 0.000 693.34 4.4e-19 24.93
12556 PRICKLE4 6_32 0.000 699.62 4.4e-19 -25.00
4956 FRS3 6_32 0.000 79.56 5.0e-20 -4.48
11112 TOMM6 6_32 0.000 361.32 1.1e-19 -15.36
7478 USP49 6_32 0.699 100.06 2.0e-04 -10.33
2751 BYSL 6_32 0.000 679.00 0.0e+00 23.99
3742 MED20 6_32 0.000 353.61 1.1e-19 -15.60
| Version | Author | Date |
|---|---|---|
| 3a7fbc1 | wesleycrouse | 2021-09-08 |
[1] "Region: 9_70"
genename region_tag susie_pip mu2 PVE z
3800 DDX31 9_70 0.000 6.73 1.5e-09 0.75
7622 SPACA9 9_70 0.000 18.70 1.3e-08 -2.39
7623 TSC1 9_70 0.000 8.77 2.7e-09 -0.77
7624 GFI1B 9_70 0.357 44.44 4.5e-05 -3.91
6002 GTF3C5 9_70 0.000 12.17 5.2e-09 -1.32
5995 GBGT1 9_70 0.000 32.47 1.5e-08 5.35
8866 ABO 9_70 0.000 275.33 7.0e-08 21.31
5998 SURF6 9_70 0.000 35.51 1.3e-08 -6.39
6001 RPL7A 9_70 0.000 44.55 1.3e-08 -7.85
5996 SURF1 9_70 0.000 18.09 4.6e-09 1.59
5997 SURF2 9_70 0.000 6.58 1.3e-09 -0.52
5994 SURF4 9_70 0.000 57.28 7.1e-08 8.36
10687 STKLD1 9_70 0.000 6.77 1.3e-09 0.65
6876 ADAMTS13 9_70 0.000 33.95 7.3e-09 -4.03
6000 REXO4 9_70 0.000 55.64 3.2e-08 7.96
3633 DBH 9_70 0.000 12.03 5.2e-09 -1.11
3632 SARDH 9_70 0.000 14.86 8.3e-09 1.50
6868 VAV2 9_70 0.000 5.25 1.1e-09 0.02
11444 LINC00094 9_70 0.251 26.02 1.9e-05 -4.65
8258 BRD3 9_70 0.002 30.10 1.7e-07 -3.17
10249 WDR5 9_70 0.000 4.89 9.8e-10 0.13
| Version | Author | Date |
|---|---|---|
| 3a7fbc1 | wesleycrouse | 2021-09-08 |
[1] "Region: 22_24"
genename region_tag susie_pip mu2 PVE z
11011 RP1-29C18.9 22_24 0.009 18.44 4.9e-07 -2.72
10106 C22orf34 22_24 0.002 5.09 2.3e-08 -0.16
1570 BRD1 22_24 0.006 13.52 2.3e-07 -0.87
1571 ZBED4 22_24 0.002 4.99 2.2e-08 -0.43
9670 CRELD2 22_24 0.007 15.42 3.1e-07 0.89
10584 PIM3 22_24 0.002 7.52 4.2e-08 -1.06
9539 ALG12 22_24 0.003 8.83 6.9e-08 0.40
12575 CITF22-49E9.3 22_24 0.003 9.35 8.4e-08 0.32
1572 MLC1 22_24 0.003 12.61 9.5e-08 -2.27
8363 TRABD 22_24 0.002 6.72 3.6e-08 -0.65
824 PANX2 22_24 0.005 14.13 2.0e-07 -1.48
825 SELENOO 22_24 0.020 17.27 9.7e-07 -1.55
4006 TUBGCP6 22_24 0.111 31.80 1.0e-05 -3.27
1573 HDAC10 22_24 0.005 19.60 2.9e-07 4.49
10069 MAPK12 22_24 0.008 23.09 5.1e-07 4.92
10932 DENND6B 22_24 0.007 16.42 3.1e-07 -1.20
10299 PLXNB2 22_24 0.002 5.80 2.9e-08 -0.67
1508 PPP6R2 22_24 0.002 5.87 2.6e-08 -0.73
1509 SBF1 22_24 0.009 18.70 4.9e-07 0.53
4007 ADM2 22_24 0.002 13.56 8.0e-08 -3.51
1513 MIOX 22_24 0.002 9.98 4.4e-08 2.45
1514 LMF2 22_24 0.111 104.82 3.3e-05 -9.24
302 NCAPH2 22_24 0.002 15.44 1.1e-07 2.66
4194 SCO2 22_24 0.011 31.90 9.6e-07 2.33
301 TYMP 22_24 1.000 336.08 9.6e-04 -20.84
9142 ODF3B 22_24 0.004 306.50 3.2e-06 -20.20
1523 CHKB 22_24 0.957 94.97 2.6e-04 -13.37
4193 KLHDC7B 22_24 0.005 38.06 5.0e-07 -7.40
12317 CTA-384D8.35 22_24 0.211 33.11 2.0e-05 -2.96
11155 SYCE3 22_24 0.021 25.64 1.6e-06 -3.11
10928 CPT1B 22_24 0.002 43.76 1.9e-07 -8.83
153 MAPK8IP2 22_24 0.003 69.60 6.2e-07 -10.28
1530 ARSA 22_24 0.003 21.00 1.6e-07 -2.44
| Version | Author | Date |
|---|---|---|
| 3a7fbc1 | wesleycrouse | 2021-09-08 |
[1] "Region: 19_11"
genename region_tag susie_pip mu2 PVE z
10404 KANK2 19_11 0.001 4.85 8.0e-09 -0.08
4154 ANGPTL8 19_11 0.001 4.91 8.1e-09 0.24
4149 DOCK6 19_11 0.001 10.27 2.7e-08 1.33
4153 TSPAN16 19_11 0.001 5.49 9.5e-09 -0.43
2093 RAB3D 19_11 0.001 8.43 1.9e-08 0.97
2094 TMEM205 19_11 0.001 7.79 1.7e-08 0.87
2096 PLPPR2 19_11 0.001 7.36 1.5e-08 0.84
8757 SWSAP1 19_11 0.001 8.14 1.8e-08 0.93
9595 CCDC159 19_11 0.001 9.28 2.4e-08 1.09
10003 EPOR 19_11 0.001 8.43 1.9e-08 0.97
10925 RGL3 19_11 0.003 18.05 1.6e-07 1.84
10530 CCDC151 19_11 0.001 5.75 1.1e-08 -0.18
4155 PRKCSH 19_11 0.011 28.32 9.2e-07 -2.59
7002 ZNF653 19_11 0.002 11.97 5.3e-08 -1.42
4150 ECSIT 19_11 0.001 6.15 1.1e-08 -0.74
4152 ELOF1 19_11 0.003 19.07 1.8e-07 1.48
1770 ACP5 19_11 0.001 6.02 1.0e-08 -0.93
10514 ZNF823 19_11 0.001 8.87 2.2e-08 1.09
9100 ZNF491 19_11 0.001 6.45 1.3e-08 -0.34
8440 ZNF440 19_11 0.001 5.33 9.2e-09 -0.25
8439 ZNF439 19_11 0.001 8.43 1.8e-08 -1.32
10596 ZNF69 19_11 0.002 15.62 1.1e-07 -0.71
10371 ZNF763 19_11 0.002 14.27 8.3e-08 -0.61
10470 ZNF433 19_11 0.002 13.79 6.1e-08 1.42
11203 ZNF844 19_11 0.001 13.06 5.3e-08 -1.36
4343 ZNF20 19_11 0.001 5.70 1.0e-08 -0.10
11889 ZNF625 19_11 0.001 4.87 8.0e-09 -0.19
10497 ZNF44 19_11 0.001 6.55 1.1e-08 1.00
10143 ZNF563 19_11 0.001 5.29 8.8e-09 0.13
10582 ZNF442 19_11 0.001 5.28 9.1e-09 -0.06
10278 ZNF799 19_11 0.002 30.05 1.9e-07 4.65
9366 ZNF443 19_11 0.001 6.74 1.1e-08 2.04
11602 ZNF709 19_11 0.001 11.78 3.9e-08 0.10
11718 ZNF564 19_11 0.001 16.40 6.1e-08 0.43
10062 ZNF490 19_11 0.001 6.61 1.2e-08 -1.67
1965 MAN2B1 19_11 0.004 53.10 5.5e-07 6.97
2103 WDR83OS 19_11 0.001 20.94 8.3e-08 3.66
3608 WDR83 19_11 0.001 15.93 3.1e-08 -4.19
1354 DHPS 19_11 0.001 7.61 1.3e-08 2.05
4342 FBXW9 19_11 0.001 11.47 2.2e-08 -3.48
2102 TNPO2 19_11 0.001 19.38 3.6e-08 5.58
3606 C19orf43 19_11 0.001 8.89 1.5e-08 -1.49
7981 PRDX2 19_11 0.001 28.41 4.8e-08 1.51
1355 HOOK2 19_11 0.001 38.73 1.2e-07 5.05
8434 JUNB 19_11 0.001 5.23 8.7e-09 0.21
1987 RNASEH2A 19_11 0.001 31.81 6.9e-08 5.73
2106 KLF1 19_11 0.991 279.50 7.9e-04 -20.21
2104 GCDH 19_11 0.001 113.01 2.0e-07 -6.18
6997 SYCE2 19_11 0.001 214.67 3.6e-07 17.71
9252 CALR 19_11 0.001 105.24 2.9e-07 9.23
9242 FARSA 19_11 0.001 214.65 3.6e-07 -17.64
9256 RAD23A 19_11 0.001 185.35 3.4e-07 13.38
148 NFIX 19_11 0.003 21.73 1.6e-07 2.37
6930 NACC1 19_11 0.001 8.88 1.5e-08 -0.56
1993 STX10 19_11 0.001 7.27 1.4e-08 -1.02
6933 IER2 19_11 0.001 11.30 1.9e-08 3.26
5469 CACNA1A 19_11 0.001 12.60 4.7e-08 -1.25
2000 CCDC130 19_11 0.002 15.18 7.4e-08 -1.51
369 MRI1 19_11 0.001 4.99 8.4e-09 0.05
2009 C19orf53 19_11 0.002 15.98 8.7e-08 1.49
4341 ZSWIM4 19_11 0.002 16.26 8.7e-08 1.62
4345 DCAF15 19_11 0.003 22.53 1.8e-07 -2.45
8420 RLN3 19_11 0.001 9.39 2.3e-08 -1.27
2012 IL27RA 19_11 0.001 9.00 1.9e-08 -1.63
5470 MISP3 19_11 0.054 31.97 4.9e-06 4.65
789 PRKACA 19_11 0.001 13.56 2.8e-08 2.94
2013 ASF1B 19_11 0.001 6.08 1.1e-08 0.89
790 ADGRL1 19_11 0.001 11.33 3.3e-08 1.46
| Version | Author | Date |
|---|---|---|
| 3a7fbc1 | wesleycrouse | 2021-09-08 |
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
17346 rs4655791 1_43 1.000 41.66 1.2e-04 -6.28
44303 rs322933 1_99 1.000 118.28 3.4e-04 -7.46
48555 rs3754140 1_108 1.000 122.51 3.5e-04 -12.18
53168 rs766167074 1_118 1.000 63985.54 1.8e-01 8.67
59322 rs12027924 1_131 1.000 69.30 2.0e-04 -2.03
59328 rs66858280 1_131 1.000 241.14 6.9e-04 -14.00
68664 rs569546056 2_19 1.000 2827.95 8.1e-03 5.27
74730 rs71422190 2_28 1.000 240.31 6.9e-04 18.35
74756 rs10495928 2_28 1.000 359.66 1.0e-03 24.83
79902 rs168565 2_40 1.000 49.41 1.4e-04 -6.52
110723 rs863678 2_106 1.000 74.19 2.1e-04 -9.82
137263 rs3933701 3_12 1.000 114.94 3.3e-04 -10.89
178606 rs9817452 3_97 1.000 34.27 9.8e-05 -5.44
208609 rs145740858 4_39 1.000 90.48 2.6e-04 6.18
208663 rs115306053 4_40 1.000 53.71 1.5e-04 -5.20
208668 rs28651467 4_40 1.000 63.37 1.8e-04 -4.12
224620 rs35518360 4_67 1.000 73.08 2.1e-04 -8.57
224686 rs13140033 4_68 1.000 54.00 1.5e-04 -7.53
225354 rs7699743 4_69 1.000 40.83 1.2e-04 5.79
275831 rs2928169 5_45 1.000 53.35 1.5e-04 -6.79
299759 rs78112077 5_91 1.000 41.94 1.2e-04 -2.84
309372 rs70995354 6_2 1.000 3215.42 9.2e-03 -3.25
309375 rs3929752 6_2 1.000 3211.78 9.2e-03 3.30
318681 rs1419642 6_23 1.000 70.48 2.0e-04 -1.74
320627 rs2394982 6_26 1.000 208.83 6.0e-04 15.40
324628 rs12664662 6_34 1.000 74.79 2.1e-04 6.77
351322 rs113033196 6_89 1.000 416.08 1.2e-03 -32.54
351481 rs199804242 6_89 1.000 31488.87 9.0e-02 5.33
352570 rs590325 6_92 1.000 88.09 2.5e-04 11.89
352580 rs151288714 6_92 1.000 230.49 6.6e-04 -16.83
352587 rs602261 6_93 1.000 61.35 1.8e-04 -5.66
363660 rs13241427 7_3 1.000 56.50 1.6e-04 7.13
381040 rs74850530 7_36 1.000 55.02 1.6e-04 -4.11
393285 rs2237572 7_56 1.000 50.47 1.4e-04 6.80
395244 rs141862851 7_61 1.000 166.80 4.8e-04 13.20
395248 rs3807471 7_61 1.000 59.95 1.7e-04 3.87
395293 rs6465770 7_61 1.000 67.01 1.9e-04 8.20
395321 rs3197597 7_61 1.000 271.94 7.8e-04 -18.46
396391 rs763798411 7_65 1.000 16070.00 4.6e-02 2.68
396394 rs10274607 7_65 1.000 16081.91 4.6e-02 -3.08
422991 rs17296501 8_23 1.000 109.00 3.1e-04 10.58
461372 rs6415788 9_4 1.000 74.62 2.1e-04 8.96
465567 rs17278406 9_12 1.000 103.86 3.0e-04 10.34
472420 rs10813816 9_25 1.000 55.94 1.6e-04 -7.38
508527 rs71007692 10_28 1.000 20197.65 5.8e-02 4.54
509988 rs11239271 10_31 1.000 237.45 6.8e-04 15.84
517623 rs6480402 10_46 1.000 543.34 1.6e-03 -9.16
517643 rs73267649 10_46 1.000 371.95 1.1e-03 5.63
548552 rs369062552 11_21 1.000 80.35 2.3e-04 7.93
548562 rs34830202 11_21 1.000 109.78 3.1e-04 -10.42
571155 rs10891252 11_66 1.000 41.89 1.2e-04 4.92
571841 rs1176746 11_67 1.000 4453.77 1.3e-02 4.15
571843 rs2307599 11_67 1.000 4454.93 1.3e-02 -3.98
579853 rs35167730 12_2 1.000 143.11 4.1e-04 -12.29
593747 rs1859440 12_30 1.000 73.64 2.1e-04 10.58
605245 rs4842610 12_53 1.000 100.25 2.9e-04 10.30
665422 rs17766755 14_29 1.000 63.47 1.8e-04 7.63
668490 rs7156583 14_34 1.000 44671.86 1.3e-01 2.48
668493 rs369107859 14_34 1.000 44644.92 1.3e-01 -6.40
688397 rs2070895 15_27 1.000 55.55 1.6e-04 7.43
690595 rs11857609 15_30 1.000 112.14 3.2e-04 11.14
691187 rs137913180 15_31 1.000 82.18 2.3e-04 -0.81
691200 rs66461959 15_31 1.000 218.54 6.2e-04 2.76
691214 rs67453880 15_31 1.000 216.76 6.2e-04 -4.87
693179 rs2955742 15_36 1.000 49.37 1.4e-04 7.04
693753 rs4886565 15_37 1.000 105.04 3.0e-04 10.50
701953 rs2238368 16_1 1.000 207.42 5.9e-04 -14.88
701981 rs59428751 16_1 1.000 106.19 3.0e-04 11.04
712590 rs62039688 16_27 1.000 80.62 2.3e-04 8.38
712672 rs17616063 16_27 1.000 40.40 1.2e-04 5.42
727852 rs57828851 17_7 1.000 62.68 1.8e-04 8.51
739577 rs11650989 17_36 1.000 82.82 2.4e-04 9.43
763882 rs144596877 18_35 1.000 49.35 1.4e-04 6.99
770609 rs34181610 19_4 1.000 37.43 1.1e-04 -5.84
774343 rs146213062 19_12 1.000 57.08 1.6e-04 7.80
795783 rs73124945 20_24 1.000 40.56 1.2e-04 6.09
800073 rs6099616 20_33 1.000 62.34 1.8e-04 7.80
801868 rs2427539 20_38 1.000 58.77 1.7e-04 5.62
802624 rs2823025 21_2 1.000 46.39 1.3e-04 6.36
810089 rs9981948 21_17 1.000 40.22 1.1e-04 -6.26
823424 rs75792643 22_20 1.000 125.47 3.6e-04 -12.33
879771 rs199677830 1_85 1.000 38386.38 1.1e-01 -5.02
909598 rs62160676 2_66 1.000 196.82 5.6e-04 18.60
948129 rs760400154 5_2 1.000 59574.69 1.7e-01 -5.88
948131 rs563200821 5_2 1.000 59576.38 1.7e-01 6.01
964181 rs1480380 6_27 1.000 105.39 3.0e-04 11.83
972348 rs1319012 6_32 1.000 160.13 4.6e-04 10.10
972499 rs112233623 6_32 1.000 635.22 1.8e-03 19.90
972500 rs9349205 6_32 1.000 928.37 2.6e-03 -29.58
986317 rs140312767 6_74 1.000 21873.31 6.2e-02 4.70
986345 rs189426171 6_74 1.000 21850.02 6.2e-02 -3.36
999480 rs2983226 6_112 1.000 11990.15 3.4e-02 2.75
999490 rs139588569 6_112 1.000 12054.11 3.4e-02 -3.30
1010160 rs200742690 7_33 1.000 21690.04 6.2e-02 10.40
1013915 rs61739556 7_62 1.000 130.76 3.7e-04 17.14
1014029 rs78054447 7_62 1.000 85.75 2.4e-04 11.25
1015361 rs147330150 7_62 1.000 88.91 2.5e-04 14.16
1023126 rs145700013 7_94 1.000 223.43 6.4e-04 4.79
1031741 rs12005199 9_5 1.000 104.38 3.0e-04 -13.13
1146684 rs9302635 16_38 1.000 95.61 2.7e-04 7.44
1153961 rs61745086 16_54 1.000 51.89 1.5e-04 -5.75
1215095 rs8107162 19_23 1.000 106.85 3.0e-04 4.95
1215120 rs7255661 19_23 1.000 100.95 2.9e-04 -12.41
1215150 rs78744187 19_23 1.000 470.91 1.3e-03 -21.23
1218281 rs61750953 19_30 1.000 58.60 1.7e-04 6.63
1234140 rs61371437 19_34 1.000 37160.08 1.1e-01 -4.94
1234152 rs374141296 19_34 1.000 37343.97 1.1e-01 4.86
1253372 rs151305716 20_31 1.000 74.08 2.1e-04 8.75
68667 rs4580350 2_19 0.999 2796.64 8.0e-03 5.25
173646 rs79308158 3_86 0.999 102.60 2.9e-04 12.22
204698 rs10029009 4_32 0.999 45.54 1.3e-04 -6.45
396402 rs4997569 7_65 0.999 16088.06 4.6e-02 -2.99
494158 rs1886296 9_73 0.999 33.05 9.4e-05 -5.48
724288 rs565911571 16_51 0.999 46.60 1.3e-04 6.74
750496 rs35954636 18_9 0.999 45.60 1.3e-04 -6.06
759130 rs2878889 18_27 0.999 37.07 1.1e-04 6.07
828732 rs1569419 1_2 0.999 63.62 1.8e-04 -8.34
887422 rs17534202 1_102 0.999 63.43 1.8e-04 -7.99
32118 rs76016895 1_73 0.998 51.39 1.5e-04 -7.21
44371 rs74448403 1_100 0.998 31.11 8.9e-05 -4.77
279486 rs244760 5_52 0.998 71.74 2.0e-04 9.82
286194 rs61215818 5_66 0.998 35.88 1.0e-04 6.08
361490 rs1472267 6_108 0.998 31.56 9.0e-05 -5.49
543692 rs72864006 11_12 0.998 32.26 9.2e-05 4.63
612538 rs7970581 12_68 0.998 31.51 9.0e-05 5.43
661108 rs2883893 14_20 0.998 47.99 1.4e-04 6.28
683716 rs28714278 15_13 0.998 69.89 2.0e-04 7.92
348687 rs10782229 6_84 0.997 41.03 1.2e-04 -3.96
348691 rs9388461 6_84 0.997 53.72 1.5e-04 5.16
809218 rs147050753 21_15 0.997 34.26 9.8e-05 5.97
859580 rs10923397 1_72 0.997 114.56 3.3e-04 4.76
934152 rs11712192 3_120 0.997 42.98 1.2e-04 -6.19
1023135 rs2536072 7_94 0.997 122.13 3.5e-04 -4.94
395007 rs149211972 7_60 0.996 29.84 8.5e-05 -4.94
509831 rs3780890 10_31 0.996 36.72 1.0e-04 -6.19
587187 rs2344157 12_18 0.996 31.46 8.9e-05 5.34
593859 rs567199163 12_30 0.996 42.69 1.2e-04 8.52
1135900 rs3809627 16_24 0.996 198.54 5.6e-04 -16.34
36857 rs9425587 1_84 0.995 35.18 1.0e-04 6.60
48758 rs6673799 1_109 0.995 35.83 1.0e-04 5.63
59656 rs4335411 1_131 0.995 34.96 9.9e-05 -4.96
487903 rs117731582 9_57 0.995 41.50 1.2e-04 -6.48
318178 rs114563774 6_22 0.994 33.83 9.6e-05 -5.41
1180199 rs2748427 17_43 0.994 130.56 3.7e-04 -13.23
273740 rs451157 5_41 0.993 41.29 1.2e-04 -6.28
317920 rs113014817 6_21 0.993 38.38 1.1e-04 -5.82
487962 rs7034523 9_57 0.993 32.83 9.3e-05 5.73
809131 rs2834259 21_14 0.993 57.46 1.6e-04 -7.56
961274 rs3892710 6_27 0.993 50.24 1.4e-04 -7.54
1180215 rs383603 17_43 0.993 52.25 1.5e-04 -8.95
48502 rs114775675 1_108 0.992 34.10 9.7e-05 -6.94
94156 rs141849010 2_69 0.992 30.80 8.7e-05 -5.44
311281 rs585312 6_7 0.992 40.45 1.1e-04 -5.80
701978 rs1203981 16_1 0.992 69.05 2.0e-04 2.42
712565 rs11641493 16_27 0.992 30.25 8.6e-05 4.89
736520 rs8073834 17_27 0.992 41.53 1.2e-04 4.59
208664 rs4864911 4_40 0.991 49.97 1.4e-04 -5.46
318900 rs3873238 6_23 0.991 37.06 1.0e-04 -6.32
351325 rs7775698 6_89 0.991 4699.90 1.3e-02 65.30
571888 rs17116384 11_67 0.990 77.11 2.2e-04 -7.88
153945 rs13071298 3_48 0.989 58.49 1.6e-04 -7.34
488452 rs12375564 9_58 0.989 27.45 7.7e-05 4.23
612008 rs12423752 12_66 0.989 40.67 1.1e-04 6.13
5037 rs603412 1_14 0.988 30.37 8.6e-05 5.25
311301 rs55792466 6_7 0.988 82.14 2.3e-04 -8.82
619614 rs4883568 12_82 0.988 57.62 1.6e-04 -8.51
736610 rs11658398 17_29 0.988 32.04 9.0e-05 -5.60
761904 rs62094231 18_31 0.988 31.23 8.8e-05 5.62
114659 rs62181701 2_112 0.987 30.26 8.5e-05 -5.26
139820 rs1667747 3_18 0.986 41.73 1.2e-04 -4.16
1154067 rs2306050 16_54 0.986 49.44 1.4e-04 -7.07
316340 rs78381359 6_16 0.985 55.29 1.6e-04 -7.73
62812 rs4669306 2_6 0.984 58.00 1.6e-04 7.59
105774 rs115200150 2_96 0.984 29.59 8.3e-05 -6.06
500061 rs737376 10_12 0.984 37.51 1.1e-04 -5.83
676897 rs35604463 14_52 0.984 43.49 1.2e-04 -6.41
33203 rs138055271 1_79 0.983 32.05 9.0e-05 -4.94
778392 rs9304832 19_22 0.983 31.67 8.9e-05 5.40
1014258 rs7778978 7_62 0.983 267.33 7.5e-04 -25.01
139640 rs9819064 3_17 0.982 26.32 7.4e-05 4.82
538888 rs2239681 11_2 0.981 31.06 8.7e-05 -5.48
279490 rs304151 5_52 0.980 60.45 1.7e-04 10.70
324783 rs758853393 6_34 0.980 35.39 9.9e-05 -5.74
11330 rs11589584 1_27 0.979 52.14 1.5e-04 -7.73
517624 rs10998724 10_46 0.979 223.67 6.2e-04 2.57
565918 rs3016490 11_55 0.979 26.60 7.4e-05 -4.83
800133 rs1328757 20_33 0.978 39.09 1.1e-04 -5.90
1010131 rs6965751 7_33 0.978 21566.80 6.0e-02 -10.60
54370 rs2884425 1_122 0.977 66.40 1.9e-04 6.51
717698 rs4396539 16_37 0.976 27.69 7.7e-05 5.09
977481 rs3025012 6_33 0.976 54.27 1.5e-04 -9.60
74757 rs17034605 2_28 0.975 107.41 3.0e-04 14.01
736234 rs2532379 17_27 0.975 58.29 1.6e-04 -1.65
692773 rs876383 15_35 0.973 29.02 8.1e-05 5.17
542824 rs11022762 11_10 0.972 46.97 1.3e-04 -6.82
794295 rs3746615 20_19 0.972 77.61 2.2e-04 -7.95
1195273 rs12609178 19_11 0.971 82.62 2.3e-04 -11.67
74704 rs35617557 2_28 0.970 63.03 1.7e-04 -5.58
615203 rs56179458 12_74 0.970 153.34 4.2e-04 -14.39
79370 rs1641155 2_39 0.968 27.84 7.7e-05 -5.41
275689 rs10060848 5_45 0.968 40.34 1.1e-04 6.16
305567 rs322351 5_103 0.968 28.78 7.9e-05 4.40
614351 rs12830847 12_73 0.968 29.63 8.2e-05 -5.50
312068 rs79248374 6_9 0.967 25.40 7.0e-05 -5.00
972451 rs3218092 6_32 0.967 287.13 7.9e-04 22.02
445653 rs72671791 8_69 0.963 28.12 7.7e-05 -5.00
612350 rs653178 12_67 0.963 232.04 6.4e-04 17.84
117704 rs67409803 2_120 0.961 27.71 7.6e-05 -5.10
394980 rs117645417 7_60 0.961 25.61 7.0e-05 -4.34
510177 rs78553679 10_31 0.961 31.05 8.5e-05 4.75
810726 rs1984021 21_18 0.959 52.14 1.4e-04 6.75
724133 rs1915035 16_50 0.958 24.89 6.8e-05 4.60
593689 rs11608702 12_30 0.957 26.28 7.2e-05 -4.94
323241 rs4713943 6_29 0.956 34.27 9.3e-05 5.89
428187 rs13274178 8_35 0.956 28.15 7.7e-05 -5.02
299725 rs74417235 5_91 0.955 36.12 9.8e-05 7.30
425083 rs34795012 8_27 0.954 33.17 9.0e-05 5.41
683715 rs11853517 15_13 0.954 32.15 8.8e-05 -4.37
147139 rs4974078 3_34 0.953 31.94 8.7e-05 -4.78
805806 rs118097076 21_8 0.953 24.67 6.7e-05 4.45
809186 rs2834321 21_15 0.953 31.43 8.5e-05 5.93
208695 rs74995222 4_40 0.951 47.85 1.3e-04 4.77
810019 rs8128857 21_16 0.951 30.71 8.3e-05 5.53
313025 rs7453037 6_11 0.949 27.16 7.4e-05 4.97
15237 rs2261980 1_38 0.948 25.17 6.8e-05 4.63
305972 rs189650 5_104 0.948 48.61 1.3e-04 -6.00
810117 rs9305604 21_17 0.948 28.54 7.7e-05 -5.13
334406 rs144917451 6_51 0.947 26.80 7.2e-05 4.95
360818 rs4709820 6_107 0.947 163.72 4.4e-04 -12.14
44332 rs567188561 1_99 0.943 36.09 9.7e-05 -7.72
701950 rs2238366 16_1 0.940 35.83 9.6e-05 -4.67
353387 rs13208190 6_94 0.937 27.52 7.4e-05 -4.94
1048366 rs782134971 9_70 0.936 591.36 1.6e-03 25.45
574401 rs497879 11_74 0.935 25.25 6.7e-05 -4.66
525602 rs11462611 10_61 0.934 50.42 1.3e-04 7.38
771368 rs72990643 19_5 0.934 36.50 9.7e-05 -6.69
527770 rs111604275 10_65 0.933 26.36 7.0e-05 -4.84
616270 rs4765552 12_75 0.932 28.79 7.7e-05 4.73
177099 rs370612 3_94 0.931 28.24 7.5e-05 5.09
736542 rs1808192 17_27 0.931 44.13 1.2e-04 4.62
446164 rs2570952 8_69 0.928 31.14 8.2e-05 5.26
318733 rs1233385 6_23 0.927 88.13 2.3e-04 7.11
114984 rs192394817 2_113 0.925 24.83 6.6e-05 -4.89
30963 rs60897900 1_71 0.924 40.01 1.1e-04 5.69
363659 rs78148157 7_3 0.922 26.68 7.0e-05 -3.90
819219 rs9609565 22_12 0.922 62.49 1.6e-04 -10.53
1042196 rs72752884 9_66 0.921 75.49 2.0e-04 -8.74
20362 rs12134068 1_48 0.920 26.43 6.9e-05 -4.74
532272 rs1925282 10_73 0.917 25.97 6.8e-05 4.73
721156 rs34252270 16_46 0.917 25.96 6.8e-05 -4.68
216523 rs11725113 4_52 0.915 27.35 7.1e-05 3.41
75004 rs2166475 2_29 0.914 40.66 1.1e-04 -5.90
571198 rs1123066 11_66 0.914 41.93 1.1e-04 -5.39
101699 rs10201197 2_86 0.913 44.82 1.2e-04 -6.93
1031738 rs409950 9_5 0.909 34.73 9.0e-05 -0.27
314276 rs7767676 6_13 0.908 57.12 1.5e-04 -7.53
1052443 rs11245982 11_1 0.908 43.14 1.1e-04 -5.83
122689 rs115997322 2_129 0.907 27.34 7.1e-05 5.63
237319 rs13106087 4_94 0.907 25.79 6.7e-05 4.61
139836 rs12632081 3_18 0.905 63.91 1.7e-04 5.66
574540 rs1945397 11_75 0.905 27.41 7.1e-05 4.90
798852 rs2072241 20_32 0.905 28.99 7.5e-05 4.90
339865 rs11757155 6_61 0.903 36.57 9.4e-05 6.05
619379 rs57565032 12_81 0.902 27.33 7.0e-05 4.90
1062308 rs10840114 11_7 0.901 69.54 1.8e-04 8.52
693077 rs11072556 15_35 0.900 30.33 7.8e-05 4.33
694404 rs4605119 15_37 0.899 31.21 8.0e-05 5.46
104992 rs78236918 2_94 0.898 25.10 6.4e-05 4.57
352588 rs113883411 6_93 0.898 82.59 2.1e-04 7.54
1107724 rs12885878 14_54 0.898 41.98 1.1e-04 -6.85
676471 rs7159884 14_51 0.896 29.41 7.5e-05 5.18
580372 rs594088 12_4 0.895 25.23 6.4e-05 4.58
403524 rs9649546 7_80 0.894 26.82 6.8e-05 4.86
461618 rs6476934 9_6 0.890 25.53 6.5e-05 4.86
722140 rs11862358 16_46 0.889 25.60 6.5e-05 -4.40
494651 rs113790047 10_3 0.888 32.48 8.2e-05 5.51
63258 rs139638572 2_6 0.887 30.90 7.8e-05 5.29
317559 rs145151542 6_19 0.887 24.89 6.3e-05 3.87
94499 rs35883727 2_70 0.884 64.51 1.6e-04 -8.23
284855 rs6884480 5_63 0.883 23.53 5.9e-05 4.37
48250 rs12408694 1_108 0.882 41.09 1.0e-04 -6.28
149079 rs137980154 3_39 0.881 23.82 6.0e-05 4.09
381038 rs12718598 7_36 0.881 297.95 7.5e-04 16.82
170005 rs2811366 3_79 0.880 24.32 6.1e-05 -4.41
286207 rs62371569 5_66 0.880 25.39 6.4e-05 -4.80
7231 rs7525656 1_19 0.879 27.63 6.9e-05 4.74
697649 rs2518967 15_43 0.879 26.75 6.7e-05 4.79
731914 rs8078135 17_16 0.879 25.43 6.4e-05 -4.79
92044 rs144442782 2_65 0.878 25.58 6.4e-05 4.85
33707 rs2494211 1_79 0.877 29.81 7.5e-05 5.05
44675 rs2737649 1_100 0.875 25.80 6.4e-05 -4.76
972997 rs2492939 6_33 0.875 49.62 1.2e-04 -6.50
351497 rs6923513 6_89 0.874 31493.25 7.9e-02 -6.32
489355 rs1861882 9_60 0.871 25.41 6.3e-05 -4.30
199066 rs963209 4_20 0.870 36.93 9.2e-05 6.64
655672 rs797342 14_8 0.870 47.24 1.2e-04 -6.12
770114 rs2238586 19_2 0.870 30.87 7.7e-05 -5.11
809183 rs4817606 21_15 0.870 33.32 8.3e-05 6.04
781008 rs309182 19_33 0.869 63.30 1.6e-04 8.14
972940 rs73735020 6_33 0.869 81.11 2.0e-04 -8.91
741156 rs147300783 17_39 0.866 25.91 6.4e-05 4.65
279262 rs7700278 5_51 0.864 26.22 6.5e-05 -4.81
138085 rs2948097 3_15 0.863 31.75 7.8e-05 -4.78
15 rs28562941 1_1 0.862 24.07 5.9e-05 4.26
613275 rs7972038 12_70 0.862 28.83 7.1e-05 4.94
12202 rs57401684 1_32 0.855 39.93 9.7e-05 5.66
322396 rs115401532 6_28 0.855 26.40 6.4e-05 4.49
700882 rs11635896 15_48 0.855 24.36 5.9e-05 -4.37
738098 rs754267414 17_32 0.855 23.95 5.8e-05 4.30
74845 rs11539645 2_29 0.853 25.20 6.1e-05 -4.73
44275 rs10489689 1_99 0.852 77.40 1.9e-04 5.80
24346 rs12130551 1_56 0.851 25.08 6.1e-05 4.70
665679 rs55792096 14_30 0.850 60.28 1.5e-04 9.83
958406 rs487624 6_20 0.850 50.90 1.2e-04 -3.46
144377 rs111998223 3_27 0.849 55.32 1.3e-04 -7.13
667061 rs1275824 14_33 0.849 28.48 6.9e-05 5.14
580334 rs61907084 12_4 0.846 25.75 6.2e-05 -4.57
544490 rs11025069 11_13 0.843 32.44 7.8e-05 -6.25
717529 rs11075747 16_37 0.843 36.73 8.8e-05 -6.34
847645 rs61784835 1_29 0.841 68.63 1.6e-04 -2.66
154007 rs7640269 3_48 0.838 27.76 6.6e-05 4.14
619598 rs4883580 12_82 0.838 24.92 6.0e-05 4.76
798793 rs2145697 20_30 0.837 59.02 1.4e-04 10.42
691506 rs16952867 15_32 0.832 31.04 7.4e-05 5.38
544240 rs66514853 11_13 0.827 27.05 6.4e-05 -4.85
727939 rs7224566 17_7 0.827 37.70 8.9e-05 6.45
809487 rs2834711 21_15 0.827 32.45 7.7e-05 -5.75
320907 rs77424484 6_26 0.826 78.45 1.8e-04 -9.32
1852 rs205474 1_9 0.822 25.07 5.9e-05 4.49
31683 rs839605 1_73 0.822 51.56 1.2e-04 -6.55
909586 rs372515880 2_66 0.820 192.26 4.5e-04 17.77
383794 rs189049917 7_40 0.819 25.08 5.9e-05 4.64
775034 rs113809670 19_14 0.817 27.54 6.4e-05 4.80
153917 rs13079951 3_48 0.816 25.33 5.9e-05 -4.54
291760 rs4463178 5_77 0.816 27.26 6.3e-05 -4.88
345085 rs13196629 6_73 0.816 80.67 1.9e-04 12.98
802662 rs926167 21_2 0.814 29.65 6.9e-05 5.75
44567 rs4915386 1_100 0.812 25.27 5.9e-05 4.29
453609 rs10956395 8_84 0.808 26.43 6.1e-05 4.85
241450 rs1830456 4_102 0.806 25.61 5.9e-05 -4.67
148999 rs6445819 3_39 0.805 144.48 3.3e-04 12.20
540861 rs2659865 11_5 0.804 25.06 5.7e-05 4.57
738256 rs3794748 17_32 0.803 28.29 6.5e-05 4.82
665671 rs8006419 14_30 0.801 63.44 1.4e-04 10.26SNPs 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 |
|---|---|---|
| 3a7fbc1 | wesleycrouse | 2021-09-08 |
#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
53168 rs766167074 1_118 1.000 63985.54 1.8e-01 8.67
53167 rs971534 1_118 0.709 63595.21 1.3e-01 -8.84
53166 rs2486737 1_118 0.570 63595.04 1.0e-01 -8.83
53165 rs10489611 1_118 0.149 63593.50 2.7e-02 -8.82
53162 rs2790891 1_118 0.131 63588.78 2.4e-02 -8.83
53163 rs2491405 1_118 0.131 63588.78 2.4e-02 -8.83
53159 rs2256908 1_118 0.074 63588.49 1.3e-02 -8.82
53175 rs2211176 1_118 0.104 63570.19 1.9e-02 -8.82
53176 rs2790882 1_118 0.104 63570.19 1.9e-02 -8.82
53174 rs2248646 1_118 0.590 63568.76 1.1e-01 -8.84
53155 rs1076804 1_118 0.000 63492.56 2.7e-05 -8.83
53177 rs1416913 1_118 0.002 63491.03 3.7e-04 -8.84
53180 rs2790874 1_118 0.020 63485.88 3.7e-03 -8.87
53156 rs910824 1_118 0.000 63329.55 2.5e-12 -8.78
53179 rs2808603 1_118 0.000 63324.80 9.5e-08 -8.89
53154 rs2474635 1_118 0.000 63137.96 0.0e+00 -8.75
53173 rs143905935 1_118 0.000 63047.79 0.0e+00 -8.65
53151 rs722302 1_118 0.000 63013.14 0.0e+00 -8.61
53171 rs2739509 1_118 0.000 62378.40 0.0e+00 -8.66
53181 rs2474631 1_118 0.000 61280.60 0.0e+00 -8.30
53178 rs3071894 1_118 0.000 61230.60 0.0e+00 -8.73
53183 rs2790871 1_118 0.000 60606.01 0.0e+00 -8.90
53184 rs2790870 1_118 0.000 60594.37 0.0e+00 -8.91
53169 rs2739512 1_118 0.000 60500.91 0.0e+00 -8.77
948131 rs563200821 5_2 1.000 59576.38 1.7e-01 6.01
948129 rs760400154 5_2 1.000 59574.69 1.7e-01 -5.88
948099 rs116876144 5_2 0.000 59489.92 6.8e-06 5.99
948095 rs79843325 5_2 0.000 59483.74 2.0e-06 5.99
948101 rs113261319 5_2 0.000 59479.17 1.7e-08 5.95
948141 rs118079687 5_2 0.000 59477.55 5.2e-07 6.00
948089 rs79635912 5_2 0.000 59440.08 5.4e-09 5.99
948080 rs117110053 5_2 0.000 59298.67 2.4e-15 6.07
53149 rs7414807 1_118 0.000 59291.85 0.0e+00 -8.11
948179 rs112060190 5_2 0.000 59265.58 0.0e+00 6.01
53186 rs371976741 1_118 0.000 59262.53 0.0e+00 8.41
948059 rs117324785 5_2 0.000 59236.40 0.0e+00 6.08
948183 rs117279826 5_2 0.000 59216.04 0.0e+00 6.02
948057 rs117169356 5_2 0.000 59214.67 0.0e+00 6.07
948056 rs117971540 5_2 0.000 59210.09 0.0e+00 6.07
948194 rs111378768 5_2 0.000 59183.75 0.0e+00 6.01
948159 rs576289416 5_2 0.000 59123.32 0.0e+00 5.94
948156 rs111206023 5_2 0.000 58951.33 0.0e+00 6.00
948124 rs554884394 5_2 0.000 58936.35 0.0e+00 6.01
948070 rs117025467 5_2 0.000 58918.52 0.0e+00 6.01
948197 rs117576772 5_2 0.000 58830.16 0.0e+00 5.98
948062 rs116875834 5_2 0.000 58823.28 0.0e+00 6.09
948139 rs372122514 5_2 0.000 58792.37 0.0e+00 -4.98
948181 rs117589613 5_2 0.000 58768.88 0.0e+00 6.06
948190 rs111460998 5_2 0.000 58735.90 0.0e+00 6.06
948012 rs116961707 5_2 0.000 58733.04 0.0e+00 6.06SNPs 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
53168 rs766167074 1_118 1.000 63985.54 0.1800 8.67
948129 rs760400154 5_2 1.000 59574.69 0.1700 -5.88
948131 rs563200821 5_2 1.000 59576.38 0.1700 6.01
53167 rs971534 1_118 0.709 63595.21 0.1300 -8.84
668490 rs7156583 14_34 1.000 44671.86 0.1300 2.48
668493 rs369107859 14_34 1.000 44644.92 0.1300 -6.40
53174 rs2248646 1_118 0.590 63568.76 0.1100 -8.84
879771 rs199677830 1_85 1.000 38386.38 0.1100 -5.02
1234140 rs61371437 19_34 1.000 37160.08 0.1100 -4.94
1234152 rs374141296 19_34 1.000 37343.97 0.1100 4.86
53166 rs2486737 1_118 0.570 63595.04 0.1000 -8.83
351481 rs199804242 6_89 1.000 31488.87 0.0900 5.33
351497 rs6923513 6_89 0.874 31493.25 0.0790 -6.32
986317 rs140312767 6_74 1.000 21873.31 0.0620 4.70
986345 rs189426171 6_74 1.000 21850.02 0.0620 -3.36
1010160 rs200742690 7_33 1.000 21690.04 0.0620 10.40
1010131 rs6965751 7_33 0.978 21566.80 0.0600 -10.60
508527 rs71007692 10_28 1.000 20197.65 0.0580 4.54
396391 rs763798411 7_65 1.000 16070.00 0.0460 2.68
396394 rs10274607 7_65 1.000 16081.91 0.0460 -3.08
396402 rs4997569 7_65 0.999 16088.06 0.0460 -2.99
508536 rs11011452 10_28 0.715 20297.86 0.0410 -4.44
999480 rs2983226 6_112 1.000 11990.15 0.0340 2.75
999490 rs139588569 6_112 1.000 12054.11 0.0340 -3.30
508526 rs2474565 10_28 0.539 20296.20 0.0310 -4.45
53165 rs10489611 1_118 0.149 63593.50 0.0270 -8.82
999501 rs59421548 6_112 0.756 11934.62 0.0260 2.83
53162 rs2790891 1_118 0.131 63588.78 0.0240 -8.83
53163 rs2491405 1_118 0.131 63588.78 0.0240 -8.83
508533 rs2472183 10_28 0.376 20296.24 0.0220 -4.45
53175 rs2211176 1_118 0.104 63570.19 0.0190 -8.82
53176 rs2790882 1_118 0.104 63570.19 0.0190 -8.82
986264 rs9400465 6_74 0.608 7880.65 0.0140 5.83
53159 rs2256908 1_118 0.074 63588.49 0.0130 -8.82
351325 rs7775698 6_89 0.991 4699.90 0.0130 65.30
571841 rs1176746 11_67 1.000 4453.77 0.0130 4.15
571843 rs2307599 11_67 1.000 4454.93 0.0130 -3.98
1010161 rs6958902 7_33 0.218 21559.32 0.0130 -10.59
351480 rs2327654 6_89 0.126 31489.40 0.0110 -6.31
1010229 rs7811411 7_33 0.412 9377.50 0.0110 11.99
309372 rs70995354 6_2 1.000 3215.42 0.0092 -3.25
309375 rs3929752 6_2 1.000 3211.78 0.0092 3.30
999496 rs62424359 6_112 0.244 11980.07 0.0083 2.89
68664 rs569546056 2_19 1.000 2827.95 0.0081 5.27
68667 rs4580350 2_19 0.999 2796.64 0.0080 5.25
986360 rs783084 6_74 0.316 8613.86 0.0078 -5.58
68663 rs7562170 2_19 0.776 2785.67 0.0062 -5.15
53180 rs2790874 1_118 0.020 63485.88 0.0037 -8.87
1010068 rs7016 7_33 0.049 21540.77 0.0030 -10.58
1010234 rs757694 7_33 0.111 9328.73 0.0030 11.95SNPs 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
351325 rs7775698 6_89 0.991 4699.90 1.3e-02 65.30
351327 rs9402685 6_89 0.009 4694.16 1.2e-04 65.22
351326 rs56293029 6_89 0.000 4589.14 0.0e+00 64.57
351332 rs6925841 6_89 0.000 910.97 4.6e-18 -39.29
351333 rs9321485 6_89 0.000 888.19 5.6e-19 -38.59
351337 rs6909822 6_89 0.000 767.67 0.0e+00 -37.04
351340 rs12204380 6_89 0.000 750.95 0.0e+00 -36.83
208656 rs218251 4_39 0.526 1062.69 1.6e-03 34.44
208653 rs218237 4_39 0.333 1061.77 1.0e-03 34.42
208655 rs218240 4_39 0.145 1059.75 4.4e-04 34.40
208654 rs218239 4_39 0.002 652.71 3.6e-06 34.35
351322 rs113033196 6_89 1.000 416.08 1.2e-03 -32.54
351278 rs6899500 6_89 0.000 1519.25 0.0e+00 -32.02
1013987 rs2075672 7_62 0.642 537.30 9.8e-04 30.80
1013979 rs9801017 7_62 0.209 534.26 3.2e-04 30.75
1013977 rs7385804 7_62 0.149 533.30 2.3e-04 30.74
1014051 rs221790 7_62 0.000 514.22 4.1e-10 30.65
1014034 rs221788 7_62 0.000 509.12 1.1e-10 30.63
1014078 rs221799 7_62 0.000 511.29 3.5e-10 30.59
1014106 rs2432930 7_62 0.000 511.80 3.9e-10 30.59
1014114 rs221772 7_62 0.000 511.63 3.8e-10 30.59
1014115 rs221771 7_62 0.000 511.82 3.9e-10 30.59
1014082 rs221801 7_62 0.000 510.13 3.0e-10 30.57
1014165 rs1617640 7_62 0.000 509.39 4.1e-10 30.54
1014128 rs10277087 7_62 0.000 507.48 2.6e-10 30.53
1014131 rs504141 7_62 0.000 507.29 2.6e-10 30.53
1014134 rs1734909 7_62 0.000 507.23 2.6e-10 30.53
1014164 rs576236 7_62 0.000 507.93 3.3e-10 30.53
1014153 rs554055 7_62 0.000 507.21 2.9e-10 30.52
1014155 rs56189543 7_62 0.000 507.30 2.9e-10 30.52
1014174 rs551238 7_62 0.000 508.13 4.1e-10 30.52
1014142 rs568733 7_62 0.000 505.83 2.2e-10 30.50
1014150 rs1734908 7_62 0.000 505.50 2.1e-10 30.50
1014162 rs35495521 7_62 0.000 506.83 2.7e-10 30.49
1014167 rs507392 7_62 0.000 504.74 2.3e-10 30.48
1014159 rs475339 7_62 0.000 504.09 1.9e-10 30.47
1014133 rs1734910 7_62 0.000 500.32 1.0e-10 30.42
351324 rs55731938 6_89 0.000 325.54 0.0e+00 -30.01
972500 rs9349205 6_32 1.000 928.37 2.6e-03 -29.58
1014027 rs371423364 7_62 0.000 475.59 3.4e-12 29.42
1013963 rs4548095 7_62 0.000 479.06 1.3e-13 29.39
1013952 rs4729597 7_62 0.000 486.04 9.7e-13 29.38
1014101 rs66481397 7_62 0.000 459.83 3.2e-12 29.38
1013956 rs35142847 7_62 0.000 477.95 4.6e-14 29.20
972452 rs1410492 6_32 0.000 887.86 1.2e-12 -29.12
972449 rs3218097 6_32 0.000 885.31 3.7e-13 -29.06
972456 rs66717417 6_32 0.000 887.16 8.7e-13 -29.05
972507 rs9394841 6_32 0.000 880.18 1.0e-13 -28.92
972596 rs9471708 6_32 0.000 702.97 0.0e+00 -28.89
972599 rs9471709 6_32 0.000 702.64 0.0e+00 -28.88Gene 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] 51if (length(genes)>0){
GO_enrichment <- enrichr(genes, dbs)
for (db in dbs){
print(db)
df <- GO_enrichment[[db]]
df <- df[df$Adjusted.P.value<0.05,c("Term", "Overlap", "Adjusted.P.value", "Genes")]
print(df)
}
#DisGeNET enrichment
# devtools::install_bitbucket("ibi_group/disgenet2r")
library(disgenet2r)
disgenet_api_key <- get_disgenet_api_key(
email = "wesleycrouse@gmail.com",
password = "uchicago1" )
Sys.setenv(DISGENET_API_KEY= disgenet_api_key)
res_enrich <-disease_enrichment(entities=genes, vocabulary = "HGNC",
database = "CURATED" )
df <- res_enrich@qresult[1:10, c("Description", "FDR", "Ratio", "BgRatio")]
print(df)
#WebGestalt enrichment
library(WebGestaltR)
background <- ctwas_gene_res$genename
#listGeneSet()
databases <- c("pathway_KEGG", "disease_GLAD4U", "disease_OMIM")
enrichResult <- WebGestaltR(enrichMethod="ORA", organism="hsapiens",
interestGene=genes, referenceGene=background,
enrichDatabase=databases, interestGeneType="genesymbol",
referenceGeneType="genesymbol", isOutput=F)
print(enrichResult[,c("description", "size", "overlap", "FDR", "database", "userId")])
}Uploading data to Enrichr... Done.
Querying GO_Biological_Process_2021... Done.
Querying GO_Cellular_Component_2021... Done.
Querying GO_Molecular_Function_2021... Done.
Parsing results... Done.
[1] "GO_Biological_Process_2021"
Term Overlap
1 pyrimidine nucleoside metabolic process (GO:0006213) 3/15
2 nucleoside metabolic process (GO:0009116) 2/10
3 pyrimidine nucleoside catabolic process (GO:0046135) 2/11
4 pyrimidine-containing compound metabolic process (GO:0072527) 2/12
5 nucleoside catabolic process (GO:0009164) 2/12
Adjusted.P.value Genes
1 0.003922708 NT5C3A;DHODH;TYMP
2 0.046692545 NT5C3A;TYMP
3 0.046692545 NT5C3A;TYMP
4 0.046692545 NT5C3A;TYMP
5 0.046692545 NT5C3A;TYMP
[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)CTC-490E21.10 gene(s) from the input list not found in DisGeNET CURATEDACTRT3 gene(s) from the input list not found in DisGeNET CURATEDCNIH4 gene(s) from the input list not found in DisGeNET CURATEDAFMID gene(s) from the input list not found in DisGeNET CURATEDSPPL2A gene(s) from the input list not found in DisGeNET CURATEDMTX1 gene(s) from the input list not found in DisGeNET CURATEDPIK3R3 gene(s) from the input list not found in DisGeNET CURATEDRNF219 gene(s) from the input list not found in DisGeNET CURATEDTEX10 gene(s) from the input list not found in DisGeNET CURATEDFCGRT gene(s) from the input list not found in DisGeNET CURATEDVAMP3 gene(s) from the input list not found in DisGeNET CURATEDFAM228B gene(s) from the input list not found in DisGeNET CURATEDSMIM1 gene(s) from the input list not found in DisGeNET CURATEDSERINC1 gene(s) from the input list not found in DisGeNET CURATEDDAGLA gene(s) from the input list not found in DisGeNET CURATEDSERTAD2 gene(s) from the input list not found in DisGeNET CURATEDLRRC8C gene(s) from the input list not found in DisGeNET CURATEDU91328.22 gene(s) from the input list not found in DisGeNET CURATEDNPIPB12 gene(s) from the input list not found in DisGeNET CURATEDDUSP14 gene(s) from the input list not found in DisGeNET CURATEDUGT8 gene(s) from the input list not found in DisGeNET CURATEDLPCAT4 gene(s) from the input list not found in DisGeNET CURATEDOSBPL7 gene(s) from the input list not found in DisGeNET CURATEDNTAN1 gene(s) from the input list not found in DisGeNET CURATEDPSKH1 gene(s) from the input list not found in DisGeNET CURATED Description FDR Ratio BgRatio
38 Hemoglobin F Disease 0.001820019 3/26 15/9703
125 Chloracne 0.010568196 3/26 33/9703
93 Helicobacter Infections 0.011461378 2/26 7/9703
4 Cooley's anemia 0.015310974 2/26 12/9703
10 beta Thalassemia 0.015310974 2/26 12/9703
100 Thalassemia Minor 0.015310974 2/26 12/9703
130 Thalassemia Intermedia 0.015310974 2/26 12/9703
8 Arenaviridae Infections 0.017791292 1/26 1/9703
108 Infections, Arenavirus 0.017791292 1/26 1/9703
116 Adult Acute Myeloblastic Leukemia 0.017791292 1/26 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!NULLSensitivity, specificity and precision for silver standard genes
library("readxl")
known_annotations <- read_xlsx("data/summary_known_genes_annotations.xlsx", sheet="RBC count")
known_annotations <- unique(known_annotations$`Gene Symbol`)
unrelated_genes <- ctwas_gene_res$genename[!(ctwas_gene_res$genename %in% known_annotations)]
#number of genes in known annotations
print(length(known_annotations))[1] 71#number of genes in known annotations with imputed expression
print(sum(known_annotations %in% ctwas_gene_res$genename))[1] 39#assign ctwas, TWAS, and bystander genes
ctwas_genes <- ctwas_gene_res$genename[ctwas_gene_res$susie_pip>0.8]
twas_genes <- ctwas_gene_res$genename[abs(ctwas_gene_res$z)>sig_thresh]
novel_genes <- ctwas_genes[!(ctwas_genes %in% twas_genes)]
#significance threshold for TWAS
print(sig_thresh)[1] 4.586539#number of ctwas genes
length(ctwas_genes)[1] 51#number of TWAS genes
length(twas_genes)[1] 538#show novel genes (ctwas genes with not in TWAS genes)
ctwas_gene_res[ctwas_gene_res$genename %in% novel_genes,report_cols] genename region_tag susie_pip mu2 PVE z
8466 LRRC8C 1_55 0.861 21.05 5.2e-05 -4.02
7345 RNF168 3_121 0.976 57.70 1.6e-04 3.42
2612 CDKN1B 12_12 0.820 22.21 5.2e-05 4.33
6255 ARL11 13_21 0.882 22.10 5.6e-05 -3.89
6253 RNF219 13_39 0.839 21.44 5.1e-05 -4.25
8999 LPCAT4 15_10 0.973 22.61 6.3e-05 -4.38
1803 DHODH 16_38 0.955 25.86 7.0e-05 3.99
3867 MMP24 20_21 0.873 25.85 6.4e-05 3.99#sensitivity / recall
sensitivity <- rep(NA,2)
names(sensitivity) <- c("ctwas", "TWAS")
sensitivity["ctwas"] <- sum(ctwas_genes %in% known_annotations)/length(known_annotations)
sensitivity["TWAS"] <- sum(twas_genes %in% known_annotations)/length(known_annotations)
sensitivity ctwas TWAS
0.04225352 0.05633803 #specificity
specificity <- rep(NA,2)
names(specificity) <- c("ctwas", "TWAS")
specificity["ctwas"] <- sum(!(unrelated_genes %in% ctwas_genes))/length(unrelated_genes)
specificity["TWAS"] <- sum(!(unrelated_genes %in% twas_genes))/length(unrelated_genes)
specificity ctwas TWAS
0.9956585 0.9517004 #precision / PPV
precision <- rep(NA,2)
names(precision) <- c("ctwas", "TWAS")
precision["ctwas"] <- sum(ctwas_genes %in% known_annotations)/length(ctwas_genes)
precision["TWAS"] <- sum(twas_genes %in% known_annotations)/length(twas_genes)
precision ctwas TWAS
0.058823529 0.007434944 #ROC curves
pip_range <- (0:1000)/1000
sensitivity <- rep(NA, length(pip_range))
specificity <- rep(NA, length(pip_range))
for (index in 1:length(pip_range)){
pip <- pip_range[index]
ctwas_genes <- ctwas_gene_res$genename[ctwas_gene_res$susie_pip>=pip]
sensitivity[index] <- sum(ctwas_genes %in% known_annotations)/length(known_annotations)
specificity[index] <- sum(!(unrelated_genes %in% ctwas_genes))/length(unrelated_genes)
}
plot(1-specificity, sensitivity, type="l", xlim=c(0,1), ylim=c(0,1))
sig_thresh_range <- seq(from=0, to=max(abs(ctwas_gene_res$z)), length.out=length(pip_range))
for (index in 1:length(sig_thresh_range)){
sig_thresh_plot <- sig_thresh_range[index]
twas_genes <- ctwas_gene_res$genename[abs(ctwas_gene_res$z)>=sig_thresh_plot]
sensitivity[index] <- sum(twas_genes %in% known_annotations)/length(known_annotations)
specificity[index] <- sum(!(unrelated_genes %in% twas_genes))/length(unrelated_genes)
}
lines(1-specificity, sensitivity, xlim=c(0,1), ylim=c(0,1), col="red", lty=2)
| Version | Author | Date |
|---|---|---|
| 3a7fbc1 | wesleycrouse | 2021-09-08 |
Sensitivity, specificity and precision for silver standard genes - bystanders only
This section first uses all silver standard genes to identify bystander genes within 1Mb. The silver standard and bystander gene lists are then subset to only genes with imputed expression in this analysis. Then, the ctwas and TWAS gene lists from this analysis are subset to only genes that are in the (subset) silver standard and bystander genes. These gene lists are then used to compute sensitivity, specificity and precision for ctwas and TWAS.
library(biomaRt)
library(GenomicRanges)Loading required package: stats4Loading required package: BiocGenericsLoading required package: parallel
Attaching package: 'BiocGenerics'The following objects are masked from 'package:parallel':
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, parApply, parCapply, parLapply,
parLapplyLB, parRapply, parSapply, parSapplyLBThe following objects are masked from 'package:stats':
IQR, mad, sd, var, xtabsThe following objects are masked from 'package:base':
anyDuplicated, append, as.data.frame, basename, cbind,
colnames, dirname, do.call, duplicated, eval, evalq, Filter,
Find, get, grep, grepl, intersect, is.unsorted, lapply, Map,
mapply, match, mget, order, paste, pmax, pmax.int, pmin,
pmin.int, Position, rank, rbind, Reduce, rownames, sapply,
setdiff, sort, table, tapply, union, unique, unsplit, which,
which.max, which.minLoading required package: S4Vectors
Attaching package: 'S4Vectors'The following object is masked from 'package:base':
expand.gridLoading required package: IRangesLoading required package: GenomeInfoDbensembl <- useEnsembl(biomart="ENSEMBL_MART_ENSEMBL", dataset="hsapiens_gene_ensembl")
G_list <- getBM(filters= "chromosome_name", attributes= c("hgnc_symbol","chromosome_name","start_position","end_position","gene_biotype"), values=1:22, mart=ensembl)
G_list <- G_list[G_list$hgnc_symbol!="",]
G_list <- G_list[G_list$gene_biotype %in% c("protein_coding","lncRNA"),]
G_list$start <- G_list$start_position
G_list$end <- G_list$end_position
G_list_granges <- makeGRangesFromDataFrame(G_list, keep.extra.columns=T)
known_annotations_positions <- G_list[G_list$hgnc_symbol %in% known_annotations,]
half_window <- 1000000
known_annotations_positions$start <- known_annotations_positions$start_position - half_window
known_annotations_positions$end <- known_annotations_positions$end_position + half_window
known_annotations_positions$start[known_annotations_positions$start<1] <- 1
known_annotations_granges <- makeGRangesFromDataFrame(known_annotations_positions, keep.extra.columns=T)
bystanders <- findOverlaps(known_annotations_granges,G_list_granges)
bystanders <- unique(subjectHits(bystanders))
bystanders <- G_list$hgnc_symbol[bystanders]
bystanders <- bystanders[!(bystanders %in% known_annotations)]
unrelated_genes <- bystanders
#number of genes in known annotations
print(length(known_annotations))[1] 71#number of genes in known annotations with imputed expression
print(sum(known_annotations %in% ctwas_gene_res$genename))[1] 39#number of bystander genes
print(length(unrelated_genes))[1] 1909#number of bystander genes with imputed expression
print(sum(unrelated_genes %in% ctwas_gene_res$genename))[1] 968#remove genes without imputed expression from gene lists
known_annotations <- known_annotations[known_annotations %in% ctwas_gene_res$genename]
unrelated_genes <- unrelated_genes[unrelated_genes %in% ctwas_gene_res$genename]
#assign ctwas and TWAS genes
ctwas_genes <- ctwas_gene_res$genename[ctwas_gene_res$susie_pip>0.8]
twas_genes <- ctwas_gene_res$genename[abs(ctwas_gene_res$z)>sig_thresh]
#significance threshold for TWAS
print(sig_thresh)[1] 4.586539#number of ctwas genes
length(ctwas_genes)[1] 51#number of ctwas genes in known annotations or bystanders
sum(ctwas_genes %in% c(known_annotations, unrelated_genes))[1] 7#number of ctwas genes
length(twas_genes)[1] 538#number of TWAS genes
sum(twas_genes %in% c(known_annotations, unrelated_genes))[1] 77#remove genes not in known or bystander lists from results
ctwas_genes <- ctwas_genes[ctwas_genes %in% c(known_annotations, unrelated_genes)]
twas_genes <- twas_genes[twas_genes %in% c(known_annotations, unrelated_genes)]
#sensitivity / recall
sensitivity <- rep(NA,2)
names(sensitivity) <- c("ctwas", "TWAS")
sensitivity["ctwas"] <- sum(ctwas_genes %in% known_annotations)/length(known_annotations)
sensitivity["TWAS"] <- sum(twas_genes %in% known_annotations)/length(known_annotations)
sensitivity ctwas TWAS
0.07692308 0.10256410 #specificity
specificity <- rep(NA,2)
names(specificity) <- c("ctwas", "TWAS")
specificity["ctwas"] <- sum(!(unrelated_genes %in% ctwas_genes))/length(unrelated_genes)
specificity["TWAS"] <- sum(!(unrelated_genes %in% twas_genes))/length(unrelated_genes)
specificity ctwas TWAS
0.9958678 0.9245868 #precision / PPV
precision <- rep(NA,2)
names(precision) <- c("ctwas", "TWAS")
precision["ctwas"] <- sum(ctwas_genes %in% known_annotations)/length(ctwas_genes)
precision["TWAS"] <- sum(twas_genes %in% known_annotations)/length(twas_genes)
precision ctwas TWAS
0.42857143 0.05194805
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] parallel stats4 stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] GenomicRanges_1.36.0 GenomeInfoDb_1.20.0 IRanges_2.18.1
[4] S4Vectors_0.22.1 BiocGenerics_0.30.0 biomaRt_2.40.1
[7] readxl_1.3.1 WebGestaltR_0.4.4 disgenet2r_0.99.2
[10] enrichR_3.0 cowplot_1.0.0 ggplot2_3.3.3
loaded via a namespace (and not attached):
[1] bitops_1.0-6 fs_1.3.1 bit64_4.0.5
[4] doParallel_1.0.16 progress_1.2.2 httr_1.4.1
[7] rprojroot_2.0.2 tools_3.6.1 doRNG_1.8.2
[10] utf8_1.2.1 R6_2.5.0 DBI_1.1.1
[13] colorspace_1.4-1 withr_2.4.1 tidyselect_1.1.0
[16] prettyunits_1.0.2 bit_4.0.4 curl_3.3
[19] compiler_3.6.1 git2r_0.26.1 Biobase_2.44.0
[22] labeling_0.3 scales_1.1.0 readr_1.4.0
[25] stringr_1.4.0 apcluster_1.4.8 digest_0.6.20
[28] rmarkdown_1.13 svglite_1.2.2 XVector_0.24.0
[31] pkgconfig_2.0.3 htmltools_0.3.6 fastmap_1.1.0
[34] rlang_0.4.11 RSQLite_2.2.7 farver_2.1.0
[37] generics_0.0.2 jsonlite_1.6 dplyr_1.0.7
[40] RCurl_1.98-1.1 magrittr_2.0.1 GenomeInfoDbData_1.2.1
[43] Matrix_1.2-18 Rcpp_1.0.6 munsell_0.5.0
[46] fansi_0.5.0 gdtools_0.1.9 lifecycle_1.0.0
[49] stringi_1.4.3 whisker_0.3-2 yaml_2.2.0
[52] zlibbioc_1.30.0 plyr_1.8.4 grid_3.6.1
[55] blob_1.2.1 promises_1.0.1 crayon_1.4.1
[58] lattice_0.20-38 hms_1.1.0 knitr_1.23
[61] pillar_1.6.1 igraph_1.2.4.1 rjson_0.2.20
[64] rngtools_1.5 reshape2_1.4.3 codetools_0.2-16
[67] XML_3.98-1.20 glue_1.4.2 evaluate_0.14
[70] data.table_1.14.0 vctrs_0.3.8 httpuv_1.5.1
[73] foreach_1.5.1 cellranger_1.1.0 gtable_0.3.0
[76] purrr_0.3.4 assertthat_0.2.1 cachem_1.0.5
[79] xfun_0.8 later_0.8.0 tibble_3.1.2
[82] iterators_1.0.13 AnnotationDbi_1.46.0 memoise_2.0.0
[85] workflowr_1.6.2 ellipsis_0.3.2