Last updated: 2021-09-09

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 59e5f4d. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Unstaged changes:
    Modified:   analysis/ukb-d-30500_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30500_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30600_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30600_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30610_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30610_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30620_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30620_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30630_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30630_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30640_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30640_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30650_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30650_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30660_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30660_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30670_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30670_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30680_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30690_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30690_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30700_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30700_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30710_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30710_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30720_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30720_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30730_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30740_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30740_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30750_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30750_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30760_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30760_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30770_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30770_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30780_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30780_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30790_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30800_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30800_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30810_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30820_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30820_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30830_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30830_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30840_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30840_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30850_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30850_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30860_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30860_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30870_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30870_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30880_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30880_irnt_Whole_Blood.Rmd
    Modified:   analysis/ukb-d-30890_irnt_Liver.Rmd
    Modified:   analysis/ukb-d-30890_irnt_Whole_Blood.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-30660_irnt_Liver_trim.Rmd) and HTML (docs/ukb-d-30660_irnt_Liver_trim.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
html cbf7408 wesleycrouse 2021-09-08 adding enrichment to reports
html 4970e3e wesleycrouse 2021-09-08 updating reports
html dfd2b5f wesleycrouse 2021-09-07 regenerating reports
html 47f58ac wesleycrouse 2021-09-06 fixing thin argument for fixed pi results
Rmd 209346f wesleycrouse 2021-09-06 updating additional analyses
html 209346f wesleycrouse 2021-09-06 updating additional analyses
html b14741c wesleycrouse 2021-09-06 switching from render to wflow_build
html 61b53b3 wesleycrouse 2021-09-06 updated PVE calculation
html 837dd01 wesleycrouse 2021-09-01 adding additional fixedsigma report
html 2739f4f wesleycrouse 2021-08-30 fixing typo
html b1e6b7e wesleycrouse 2021-08-30 fixing alignment on index
html d7dfe76 wesleycrouse 2021-08-30 Adding detailed reports for 30660
Rmd ea2e654 wesleycrouse 2021-08-30 Exploring fixed priors and trimming large z scores

Overview

These are the results of a ctwas analysis of the UK Biobank trait Direct bilirubin (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-30660_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] 10896
#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  763  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.8367291

Load 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$mu/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
b14741c wesleycrouse 2021-09-06
#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 
3.276241e-02 4.791019e-05 
#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 
 5.907624 26.209767 
#report sample size
print(sample_size)
[1] 292933
#report group size
group_size <- c(nrow(ctwas_gene_res), n_snps)
print(group_size)
[1]   10896 8696790
#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.007199253 0.037280498 
#compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))
[1] 0.0489526 0.5016989

Genes with highest PIPs

#distribution of PIPs
hist(ctwas_gene_res$susie_pip, xlim=c(0,1), main="Distribution of Gene PIPs")

Version Author Date
b14741c wesleycrouse 2021-09-06
#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
3212          CCND2       12_4     0.996 26.57 9.0e-05  5.34
7040          INHBB       2_70     0.992 22.43 7.6e-05  4.81
12467 RP11-219B17.3      15_27     0.985 44.27 1.5e-04  7.18
3562         ACVR1C       2_94     0.981 21.16 7.1e-05  4.62
1320        CWF19L1      10_64     0.981 27.59 9.2e-05 -7.09
11790        CYP2A6      19_28     0.975 21.17 7.0e-05 -4.73
5563          ABCG8       2_27     0.966 31.05 1.0e-04  5.88
2359          ABCC3      17_29     0.960 18.71 6.1e-05  4.38
4269          ITGB4      17_42     0.955 20.21 6.6e-05 -4.91
12687   RP4-781K5.7      1_121     0.938 18.26 5.8e-05 -4.17
5978        ZC3H12C      11_65     0.912 18.40 5.7e-05 -4.19
10495         PRMT6       1_66     0.911 24.80 7.7e-05  5.14
2924          EFHD1      2_136     0.910 26.06 8.1e-05  6.05
1146         DNMT3B      20_19     0.910 17.17 5.3e-05 -3.98
1848          CD276      15_35     0.908 32.47 1.0e-04  6.13
1120           CETP      16_31     0.904 18.25 5.6e-05 -4.03
1231         PABPC4       1_24     0.901 20.87 6.4e-05  4.52
1153           TGDS      13_47     0.885 17.41 5.3e-05 -4.00
6682         CYB5R1      1_102     0.868 17.41 5.2e-05 -3.95
5089         SCAF11      12_29     0.856 18.24 5.3e-05  4.13
6093        CSNK1G3       5_75     0.855 21.30 6.2e-05  4.74
10212          IL27      16_23     0.849 21.58 6.3e-05 -4.76
666           COASY      17_25     0.829 17.68 5.0e-05 -3.97
11669 RP11-452H21.4      11_43     0.816 29.71 8.3e-05  5.78
5544          CNIH4      1_114     0.810 17.41 4.8e-05 -3.78
1494           TXN2      22_14     0.803 16.69 4.6e-05  3.67
7547           LIPC      15_26     0.802 18.54 5.1e-05  3.99

Genes 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
b14741c wesleycrouse 2021-09-06
#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
12683    HCP5B       6_24     0.000 10425.85 2.7e-07  -3.52
10663   TRIM31       6_24     0.000  5437.97 9.8e-11   1.07
4833     FLOT1       6_24     0.000  5168.50 8.4e-11  -1.07
1088     USP40      2_137     0.000  3395.97 0.0e+00 -46.64
10651    ABCF1       6_24     0.000  2380.05 7.7e-09  -3.76
5766   PPP1R18       6_24     0.000  2074.46 1.3e-08  -3.94
4836       NRM       6_24     0.000  1159.31 6.1e-12  -0.40
10667    HLA-G       6_24     0.762  1150.05 3.0e-03  -6.69
624      ZNRD1       6_24     0.000   816.36 5.0e-12   0.19
879       DGKD      2_137     0.000   796.64 0.0e+00  -0.07
10747  SLCO1B7      12_16     0.000   624.38 0.0e+00  12.26
3556     HJURP      2_137     0.000   373.30 0.0e+00  10.96
10661   TRIM10       6_24     0.000   346.69 2.0e-12  -0.49
2584   SLCO1B3      12_16     0.000   287.51 0.0e+00   9.93
10648 C6orf136       6_24     0.000   155.03 8.4e-13   0.11
11136    HCG20       6_24     0.000   150.97 4.5e-12  -1.98
4482       SPX      12_16     0.000   150.77 0.0e+00   4.60
10664    RNF39       6_24     0.000   134.08 1.1e-12  -1.02
36       RECQL      12_16     0.000   129.06 0.0e+00   3.81
823      TACR2      10_46     0.014   107.78 5.3e-06   3.90

Genes 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
10667         HLA-G       6_24     0.762 1150.05 3.0e-03 -6.69
12467 RP11-219B17.3      15_27     0.985   44.27 1.5e-04  7.18
5563          ABCG8       2_27     0.966   31.05 1.0e-04  5.88
1848          CD276      15_35     0.908   32.47 1.0e-04  6.13
537           DGAT2      11_42     0.608   46.46 9.6e-05 -7.51
1320        CWF19L1      10_64     0.981   27.59 9.2e-05 -7.09
3212          CCND2       12_4     0.996   26.57 9.0e-05  5.34
11669 RP11-452H21.4      11_43     0.816   29.71 8.3e-05  5.78
2924          EFHD1      2_136     0.910   26.06 8.1e-05  6.05
2004          TGFB1      19_28     0.731   31.30 7.8e-05  5.64
10495         PRMT6       1_66     0.911   24.80 7.7e-05  5.14
7040          INHBB       2_70     0.992   22.43 7.6e-05  4.81
3562         ACVR1C       2_94     0.981   21.16 7.1e-05  4.62
11790        CYP2A6      19_28     0.975   21.17 7.0e-05 -4.73
4269          ITGB4      17_42     0.955   20.21 6.6e-05 -4.91
8142         CNTROB       17_7     0.685   27.81 6.5e-05 -5.71
1231         PABPC4       1_24     0.901   20.87 6.4e-05  4.52
10212          IL27      16_23     0.849   21.58 6.3e-05 -4.76
6093        CSNK1G3       5_75     0.855   21.30 6.2e-05  4.74
2359          ABCC3      17_29     0.960   18.71 6.1e-05  4.38

Genes 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
1088          USP40      2_137     0.000 3395.97 0.0e+00 -46.64
10747       SLCO1B7      12_16     0.000  624.38 0.0e+00  12.26
3556          HJURP      2_137     0.000  373.30 0.0e+00  10.96
8651           MSL2       3_84     0.091   81.97 2.5e-05  10.28
2584        SLCO1B3      12_16     0.000  287.51 0.0e+00   9.93
2586         GOLT1B      12_16     0.000   62.28 0.0e+00   7.53
537           DGAT2      11_42     0.608   46.46 9.6e-05  -7.51
11290  MAPKAPK5-AS1      12_67     0.127   37.01 1.6e-05  -7.21
12467 RP11-219B17.3      15_27     0.985   44.27 1.5e-04   7.18
2541          ALDH2      12_67     0.107   35.10 1.3e-05   7.10
1320        CWF19L1      10_64     0.981   27.59 9.2e-05  -7.09
2536          SH2B3      12_67     0.055   30.20 5.7e-06   6.80
10667         HLA-G       6_24     0.762 1150.05 3.0e-03  -6.69
2170            AHR       7_17     0.052   30.99 5.5e-06  -6.58
4962          EXOC6      10_59     0.099   42.61 1.4e-05  -6.37
8505         HECTD4      12_67     0.209   32.59 2.3e-05   6.33
1191          ERP29      12_67     0.185   29.97 1.9e-05   6.25
10370       TMEM116      12_67     0.185   29.97 1.9e-05  -6.25
9829        ZKSCAN4       6_22     0.095   21.43 6.9e-06  -6.24
10425        AKR1C4       10_6     0.198   31.54 2.1e-05   6.16

Comparing 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
b14741c wesleycrouse 2021-09-06
#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
b14741c wesleycrouse 2021-09-06
#proportion of significant z scores
mean(abs(ctwas_gene_res$z) > sig_thresh)
[1] 0.006240822
#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
1088          USP40      2_137     0.000 3395.97 0.0e+00 -46.64
10747       SLCO1B7      12_16     0.000  624.38 0.0e+00  12.26
3556          HJURP      2_137     0.000  373.30 0.0e+00  10.96
8651           MSL2       3_84     0.091   81.97 2.5e-05  10.28
2584        SLCO1B3      12_16     0.000  287.51 0.0e+00   9.93
2586         GOLT1B      12_16     0.000   62.28 0.0e+00   7.53
537           DGAT2      11_42     0.608   46.46 9.6e-05  -7.51
11290  MAPKAPK5-AS1      12_67     0.127   37.01 1.6e-05  -7.21
12467 RP11-219B17.3      15_27     0.985   44.27 1.5e-04   7.18
2541          ALDH2      12_67     0.107   35.10 1.3e-05   7.10
1320        CWF19L1      10_64     0.981   27.59 9.2e-05  -7.09
2536          SH2B3      12_67     0.055   30.20 5.7e-06   6.80
10667         HLA-G       6_24     0.762 1150.05 3.0e-03  -6.69
2170            AHR       7_17     0.052   30.99 5.5e-06  -6.58
4962          EXOC6      10_59     0.099   42.61 1.4e-05  -6.37
8505         HECTD4      12_67     0.209   32.59 2.3e-05   6.33
1191          ERP29      12_67     0.185   29.97 1.9e-05   6.25
10370       TMEM116      12_67     0.185   29.97 1.9e-05  -6.25
9829        ZKSCAN4       6_22     0.095   21.43 6.9e-06  -6.24
10425        AKR1C4       10_6     0.198   31.54 2.1e-05   6.16

Locus 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: 2_137"
        genename region_tag susie_pip     mu2 PVE      z
10567     GIGYF2      2_137         0   54.05   0  -5.12
9340     C2orf82      2_137         0    5.37   0   0.45
620         NGEF      2_137         0   44.21   0   2.52
8006      INPP5D      2_137         0   36.50   0   4.05
879         DGKD      2_137         0  796.64   0  -0.07
1088       USP40      2_137         0 3395.97   0 -46.64
3556       HJURP      2_137         0  373.30   0  10.96
11098 AC006037.2      2_137         0   15.79   0   1.50

Version Author Date
b14741c wesleycrouse 2021-09-06
[1] "Region: 12_16"
      genename region_tag susie_pip    mu2 PVE     z
2584   SLCO1B3      12_16         0 287.51   0  9.93
10747  SLCO1B7      12_16         0 624.38   0 12.26
3400      IAPP      12_16         0  95.90   0  5.16
3399   PYROXD1      12_16         0  11.90   0  0.89
36       RECQL      12_16         0 129.06   0  3.81
2586    GOLT1B      12_16         0  62.28   0  7.53
4482       SPX      12_16         0 150.77   0  4.60
2587      LDHB      12_16         0   8.40   0  0.07
3401     KCNJ8      12_16         0   6.08   0 -0.33
689      ABCC9      12_16         0   7.00   0  2.11
2590     C2CD5      12_16         0  10.33   0 -1.39
5073     ETNK1      12_16         0   9.07   0  0.30

Version Author Date
b14741c wesleycrouse 2021-09-06
[1] "Region: 3_84"
     genename region_tag susie_pip   mu2     PVE     z
796   PPP2R3A       3_84     0.145 13.81 6.9e-06 -2.46
8651     MSL2       3_84     0.091 81.97 2.5e-05 10.28
2795     PCCB       3_84     0.064  5.48 1.2e-06  1.22
3148    STAG1       3_84     0.066  4.63 1.0e-06 -0.03
6584     NCK1       3_84     0.064  7.94 1.7e-06 -2.19

Version Author Date
b14741c wesleycrouse 2021-09-06
[1] "Region: 11_42"
           genename region_tag susie_pip   mu2     PVE     z
7611          XRRA1      11_42     0.120  7.42 3.0e-06 -1.00
3170          SPCS2      11_42     0.095  5.54 1.8e-06  0.73
6901           NEU3      11_42     0.085  4.45 1.3e-06  0.39
4848        SLCO2B1      11_42     0.085  4.44 1.3e-06 -0.22
12001         TPBGL      11_42     0.110  6.57 2.5e-06  0.68
6617          GDPD5      11_42     0.120  9.34 3.8e-06  1.95
8328           MAP6      11_42     0.098  5.57 1.9e-06  0.24
7603         MOGAT2      11_42     0.264 12.77 1.2e-05  0.85
537           DGAT2      11_42     0.608 46.46 9.6e-05 -7.51
10381         UVRAG      11_42     0.095  6.40 2.1e-06  1.62
1082          WNT11      11_42     0.117  7.35 2.9e-06  1.06
11773 RP11-619A14.3      11_42     0.095  5.41 1.8e-06  0.61
4849         THAP12      11_42     0.091  5.17 1.6e-06 -0.66
12265 RP11-111M22.5      11_42     0.099  5.88 2.0e-06  0.80
11766 RP11-111M22.3      11_42     0.085  4.48 1.3e-06  0.31
11751  RP11-672A2.4      11_42     0.092  5.08 1.6e-06  0.53
9350           TSKU      11_42     0.089  4.75 1.4e-06  0.25
905           ACER3      11_42     0.265 14.43 1.3e-05 -2.04
5976          CAPN5      11_42     0.196 11.97 8.0e-06  1.72

Version Author Date
b14741c wesleycrouse 2021-09-06
[1] "Region: 12_67"
          genename region_tag susie_pip   mu2     PVE     z
5112          TCHP      12_67     0.495 15.40 2.6e-05  3.21
5111          GIT2      12_67     0.297 14.18 1.4e-05  3.36
8639      C12orf76      12_67     0.069  6.01 1.4e-06 -0.30
3515         IFT81      12_67     0.074  9.67 2.4e-06  2.50
10093       ANAPC7      12_67     0.072  7.73 1.9e-06  2.16
2531         ARPC3      12_67     0.078  6.76 1.8e-06  0.54
10684      FAM216A      12_67     0.063  7.90 1.7e-06  2.42
2532          GPN3      12_67     0.055  4.79 8.9e-07  1.40
2533         VPS29      12_67     0.055  4.82 9.0e-07 -1.41
10683        TCTN1      12_67     0.069  5.68 1.3e-06  0.02
3517         HVCN1      12_67     0.108 11.91 4.4e-06  2.67
9717        PPP1CC      12_67     0.108 11.85 4.4e-06 -2.65
10375      FAM109A      12_67     0.057  6.81 1.3e-06 -1.46
2536         SH2B3      12_67     0.055 30.20 5.7e-06  6.80
10680        ATXN2      12_67     0.057 17.13 3.3e-06  3.97
2541         ALDH2      12_67     0.107 35.10 1.3e-05  7.10
11290 MAPKAPK5-AS1      12_67     0.127 37.01 1.6e-05 -7.21
1191         ERP29      12_67     0.185 29.97 1.9e-05  6.25
10370      TMEM116      12_67     0.185 29.97 1.9e-05 -6.25
2544         NAA25      12_67     0.155 27.83 1.5e-05 -6.12
8505        HECTD4      12_67     0.209 32.59 2.3e-05  6.33
9084        PTPN11      12_67     0.072  6.64 1.6e-06 -0.78

Version Author Date
b14741c wesleycrouse 2021-09-06

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
326005  rs72834643       6_20     1.000   118.88 4.1e-04   9.73
326026 rs115740542       6_20     1.000   180.46 6.2e-04  12.66
370129  rs12208357      6_103     1.000    80.73 2.8e-04  -6.65
381915 rs542176135       7_17     1.000   152.94 5.2e-04  -8.38
533135   rs6480402      10_46     1.000   469.87 1.6e-03   8.90
603061  rs11045819      12_16     1.000  2552.27 8.7e-03 -14.34
603078   rs4363657      12_16     1.000  1590.67 5.4e-03  43.78
797812  rs59616136      19_14     1.000    49.81 1.7e-04   7.00
897096  rs13031505      2_136     1.000    60.72 2.1e-04  -8.07
900079  rs13402456      2_137     1.000  3952.27 1.3e-02 -44.40
900100  rs28948388      2_137     1.000  2663.80 9.1e-03 -35.74
900242   rs6431631      2_137     1.000  8372.22 2.9e-02  49.75
900289 rs150783152      2_137     1.000  3828.43 1.3e-02  45.10
909053   rs1611236       6_24     1.000 27586.62 9.4e-02  -3.60
326868   rs3130253       6_23     0.999    41.11 1.4e-04  -6.63
508420 rs115478735       9_70     0.999    65.26 2.2e-04  -8.02
381937   rs4721597       7_17     0.996    89.08 3.0e-04   1.94
554431  rs76153913       11_2     0.996    47.07 1.6e-04   6.70
805958 rs113345881      19_32     0.996    35.92 1.2e-04   6.12
629891    rs653178      12_67     0.995   155.80 5.3e-04 -13.12
897033  rs34247311      2_136     0.995    38.53 1.3e-04   4.44
613305   rs7397189      12_36     0.993    35.11 1.2e-04   5.69
817689  rs34507316      20_13     0.993    33.09 1.1e-04   5.60
805956    rs814573      19_32     0.989    39.63 1.3e-04  -6.68
36470    rs2779116       1_78     0.981    68.79 2.3e-04  -8.25
744704   rs2608604      16_53     0.976    55.78 1.9e-04  -6.33
796041 rs141645070      19_10     0.969    30.48 1.0e-04  -5.27
798597   rs3794991      19_15     0.952   131.75 4.3e-04  11.64
805892   rs1551891      19_32     0.946    35.90 1.2e-04   7.88
845995  rs34662558      22_10     0.935    30.92 9.9e-05  -5.19
602776   rs7962574      12_15     0.922    45.95 1.4e-04  -8.40
805961  rs12721109      19_32     0.918    33.02 1.0e-04   6.44
494338   rs9410381       9_45     0.913    74.54 2.3e-04   8.62
233289  rs17238095       4_72     0.910    29.26 9.1e-05   5.18
436846  rs12549737       8_24     0.896    30.05 9.2e-05   5.15
508805  rs34755157       9_71     0.896    29.05 8.9e-05  -5.03
900250  rs35025179      2_137     0.890  5886.79 1.8e-02 -49.78
558884  rs34623292      11_10     0.888    30.00 9.1e-05  -5.06
602781  rs73080739      12_15     0.879    34.16 1.0e-04  -7.24
325844  rs75080831       6_19     0.854    60.06 1.8e-04   7.96
602816  rs10770693      12_15     0.851    59.30 1.7e-04   8.86
294604   rs4566840       5_66     0.840    31.47 9.0e-05  -5.46
557453   rs4910498       11_8     0.833    50.56 1.4e-04   6.71

SNPs 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
b14741c wesleycrouse 2021-09-06
#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
909053   rs1611236       6_24     1.000 27586.62 9.4e-02 -3.60
909076   rs1611248       6_24     0.237 27585.63 2.2e-02 -3.63
909080   rs1611252       6_24     0.140 27585.51 1.3e-02 -3.62
909035 rs111734624       6_24     0.146 27585.47 1.4e-02 -3.62
909032   rs2508055       6_24     0.146 27585.46 1.4e-02 -3.62
909097   rs1611260       6_24     0.138 27585.44 1.3e-02 -3.62
909103   rs1611265       6_24     0.137 27585.42 1.3e-02 -3.62
908971   rs1633033       6_24     0.253 27585.30 2.4e-02 -3.63
908984   rs2844838       6_24     0.129 27584.82 1.2e-02 -3.63
909105   rs1611267       6_24     0.057 27584.79 5.4e-03 -3.62
909107   rs2394171       6_24     0.049 27584.76 4.6e-03 -3.61
909109   rs2893981       6_24     0.049 27584.71 4.6e-03 -3.61
909039   rs1611228       6_24     0.063 27584.64 5.9e-03 -3.62
909028   rs1737020       6_24     0.020 27584.20 1.8e-03 -3.60
909029   rs1737019       6_24     0.020 27584.20 1.8e-03 -3.60
908988   rs1633032       6_24     0.066 27582.93 6.2e-03 -3.62
909022   rs1633020       6_24     0.026 27581.16 2.4e-03 -3.61
909026   rs1633018       6_24     0.009 27580.58 8.2e-04 -3.60
909051   rs1611234       6_24     0.009 27578.97 8.1e-04 -3.60
908911   rs1610726       6_24     0.014 27577.73 1.3e-03 -3.61
908979   rs2844840       6_24     0.034 27575.49 3.2e-03 -3.62
909306   rs3129185       6_24     0.041 27575.16 3.9e-03 -3.63
909321   rs1736999       6_24     0.029 27573.96 2.8e-03 -3.63
909334   rs1633001       6_24     0.055 27572.06 5.2e-03 -3.64
909074   rs1611246       6_24     0.014 27570.41 1.4e-03 -3.62
909510   rs1632980       6_24     0.025 27570.15 2.4e-03 -3.63
909007   rs1614309       6_24     0.016 27562.78 1.5e-03 -3.63
909006   rs1633030       6_24     0.000 27539.27 5.5e-06 -3.61
909119   rs9258382       6_24     0.001 27512.70 9.1e-05 -3.69
909116   rs9258379       6_24     0.000 27465.77 3.9e-09 -3.62
909065   rs1611241       6_24     0.004 27439.03 4.0e-04 -3.84
909010   rs1633028       6_24     0.000 27397.93 5.6e-09 -3.72
909068   rs1611244       6_24     0.000 27291.46 3.4e-09 -3.64
909023   rs2735042       6_24     0.000 27248.96 3.5e-09 -3.61
909104   rs1611266       6_24     0.000 27051.57 1.6e-08 -3.97
909077   rs1611249       6_24     0.000 26937.56 3.0e-08 -4.10
909043   rs1611230       6_24     0.000 26870.84 2.6e-08 -4.08
909092 rs145043018       6_24     0.000 26865.68 2.8e-08 -4.09
909102 rs147376303       6_24     0.000 26865.55 2.8e-08 -4.09
909113   rs9258376       6_24     0.000 26865.21 2.8e-08 -4.09
909120   rs1633016       6_24     0.000 26864.88 2.8e-08 -4.09
908965   rs1633035       6_24     0.000 26864.28 2.8e-08 -4.09
909173   rs1633014       6_24     0.000 26862.45 2.7e-08 -4.09
908998   rs1618670       6_24     0.000 26862.21 2.7e-08 -4.09
909025   rs1633019       6_24     0.000 26859.77 2.4e-08 -4.06
909287   rs1633010       6_24     0.000 26854.25 2.7e-08 -4.09
909411    rs909722       6_24     0.000 26849.65 2.6e-08 -4.08
909443   rs1610713       6_24     0.000 26849.60 2.6e-08 -4.09
909368   rs1736991       6_24     0.000 26849.44 2.7e-08 -4.09
909423   rs1610648       6_24     0.000 26848.15 2.8e-08 -4.10

SNPs 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
909053   rs1611236       6_24     1.000 27586.62 0.09400  -3.60
900242   rs6431631      2_137     1.000  8372.22 0.02900  49.75
908971   rs1633033       6_24     0.253 27585.30 0.02400  -3.63
909076   rs1611248       6_24     0.237 27585.63 0.02200  -3.63
900250  rs35025179      2_137     0.890  5886.79 0.01800 -49.78
909032   rs2508055       6_24     0.146 27585.46 0.01400  -3.62
909035 rs111734624       6_24     0.146 27585.47 0.01400  -3.62
900079  rs13402456      2_137     1.000  3952.27 0.01300 -44.40
900289 rs150783152      2_137     1.000  3828.43 0.01300  45.10
909080   rs1611252       6_24     0.140 27585.51 0.01300  -3.62
909097   rs1611260       6_24     0.138 27585.44 0.01300  -3.62
909103   rs1611265       6_24     0.137 27585.42 0.01300  -3.62
908984   rs2844838       6_24     0.129 27584.82 0.01200  -3.63
900100  rs28948388      2_137     1.000  2663.80 0.00910 -35.74
603061  rs11045819      12_16     1.000  2552.27 0.00870 -14.34
908988   rs1633032       6_24     0.066 27582.93 0.00620  -3.62
602962  rs12366506      12_16     0.714  2465.44 0.00600  23.44
909039   rs1611228       6_24     0.063 27584.64 0.00590  -3.62
603078   rs4363657      12_16     1.000  1590.67 0.00540  43.78
909105   rs1611267       6_24     0.057 27584.79 0.00540  -3.62
909334   rs1633001       6_24     0.055 27572.06 0.00520  -3.64
909107   rs2394171       6_24     0.049 27584.76 0.00460  -3.61
909109   rs2893981       6_24     0.049 27584.71 0.00460  -3.61
602968  rs11045612      12_16     0.532  2463.24 0.00450  23.44
909306   rs3129185       6_24     0.041 27575.16 0.00390  -3.63
908979   rs2844840       6_24     0.034 27575.49 0.00320  -3.62
909321   rs1736999       6_24     0.029 27573.96 0.00280  -3.63
909022   rs1633020       6_24     0.026 27581.16 0.00240  -3.61
909510   rs1632980       6_24     0.025 27570.15 0.00240  -3.63
900253  rs17863805      2_137     0.103  5879.66 0.00210 -49.75
602979  rs73062442      12_16     0.215  2463.49 0.00180  23.39
909028   rs1737020       6_24     0.020 27584.20 0.00180  -3.60
909029   rs1737019       6_24     0.020 27584.20 0.00180  -3.60
602971  rs11513221      12_16     0.198  2463.53 0.00170  23.38
533135   rs6480402      10_46     1.000   469.87 0.00160   8.90
909007   rs1614309       6_24     0.016 27562.78 0.00150  -3.63
909074   rs1611246       6_24     0.014 27570.41 0.00140  -3.62
908911   rs1610726       6_24     0.014 27577.73 0.00130  -3.61
909026   rs1633018       6_24     0.009 27580.58 0.00082  -3.60
909051   rs1611234       6_24     0.009 27578.97 0.00081  -3.60
533140  rs35233497      10_46     0.525   352.95 0.00063  -4.24
326026 rs115740542       6_20     1.000   180.46 0.00062  12.66
533143  rs79086908      10_46     0.474   352.88 0.00057  -4.24
629891    rs653178      12_67     0.995   155.80 0.00053 -13.12
381915 rs542176135       7_17     1.000   152.94 0.00052  -8.38
602872   rs3060461      12_16     0.317   477.25 0.00052 -21.86
914179   rs3132129       6_24     0.094  1434.35 0.00046  -7.09
602876  rs35290079      12_16     0.279   474.92 0.00045 -21.76
798597   rs3794991      19_15     0.952   131.75 0.00043  11.64
326005  rs72834643       6_20     1.000   118.88 0.00041   9.73

SNPs 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
900250  rs35025179      2_137     0.890  5886.79 1.8e-02 -49.78
900242   rs6431631      2_137     1.000  8372.22 2.9e-02  49.75
900253  rs17863805      2_137     0.103  5879.66 2.1e-03 -49.75
900246   rs9287649      2_137     0.001  5867.77 1.0e-05 -49.74
900245   rs4663973      2_137     0.006  5873.92 1.3e-04 -49.71
900244   rs4663335      2_137     0.000  5862.95 1.3e-07 -49.63
899151  rs62192912      2_137     0.000  3896.63 0.0e+00 -49.52
900224   rs1500482      2_137     0.000  8260.81 0.0e+00  49.20
899798  rs72980395      2_137     0.000  6695.23 0.0e+00  48.79
900267  rs12468017      2_137     0.000  5622.20 0.0e+00 -48.71
900247   rs6431632      2_137     0.000  8103.88 0.0e+00  48.44
900260  rs17862885      2_137     0.000  5661.36 0.0e+00 -48.44
900220      rs8330      2_137     0.000  8057.30 0.0e+00  48.40
900249  rs10199525      2_137     0.000  8103.55 0.0e+00  48.32
899649   rs7584554      2_137     0.000  4661.10 0.0e+00  48.31
900262 rs116149023      2_137     0.000  5614.64 0.0e+00 -48.31
900254  rs11903524      2_137     0.000  8101.74 0.0e+00  48.30
900248   rs6431633      2_137     0.000  8100.21 0.0e+00  48.29
900251   rs7583811      2_137     0.000  8100.57 0.0e+00  48.29
900257   rs7596785      2_137     0.000  8097.59 0.0e+00  48.26
899999  rs62192806      2_137     0.000  6434.77 0.0e+00  48.16
900258  rs10929304      2_137     0.000  8076.73 0.0e+00  48.13
899620    rs838716      2_137     0.000  2530.54 0.0e+00  48.10
899968  rs62192778      2_137     0.000  4587.29 0.0e+00 -48.02
899602    rs838711      2_137     0.000  2532.49 0.0e+00  47.95
899956  rs62192777      2_137     0.000  4568.29 0.0e+00 -47.92
899596    rs838709      2_137     0.000  2502.54 0.0e+00  47.89
899595    rs838708      2_137     0.000  2499.08 0.0e+00  47.88
900329   rs6716168      2_137     0.000  7644.44 0.0e+00 -47.87
900221  rs17862880      2_137     0.000  5324.94 0.0e+00 -47.84
900206  rs55891750      2_137     0.000 10773.35 0.0e+00  47.71
900326    rs730673      2_137     0.000  7554.83 0.0e+00 -47.70
899842 rs142848041      2_137     0.000  6644.41 0.0e+00  47.35
900264  rs80334668      2_137     0.000  6138.56 0.0e+00  47.35
900226  rs11563250      2_137     0.000  5271.27 0.0e+00 -47.28
900227  rs28900409      2_137     0.000  5270.77 0.0e+00 -47.26
899819 rs139170694      2_137     0.000  6586.89 0.0e+00  47.16
900109  rs12469671      2_137     0.000  4785.04 0.0e+00 -47.07
899793 rs144648801      2_137     0.000  7414.14 0.0e+00  47.00
900255 rs561956448      2_137     0.000  5211.67 0.0e+00 -47.00
900218  rs10929303      2_137     0.000  7618.29 0.0e+00  46.69
899559   rs6727030      2_137     0.000  2362.38 0.0e+00  46.64
899732  rs76218013      2_137     0.000  5132.55 0.0e+00  46.51
899544   rs1550532      2_137     0.000  2338.27 0.0e+00  46.49
899251  rs62192932      2_137     0.000  3449.63 0.0e+00  46.41
900231   rs7586006      2_137     0.000  5197.26 0.0e+00 -46.39
900229  rs10199882      2_137     0.000  5200.72 0.0e+00 -46.38
899704  rs11689634      2_137     0.000  7211.50 0.0e+00  46.34
899738  rs72980357      2_137     0.000  7107.81 0.0e+00  46.06
899983  rs78198652      2_137     0.000  3468.90 0.0e+00 -45.92

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] cowplot_1.0.0 ggplot2_3.3.3

loaded via a namespace (and not attached):
 [1] Biobase_2.44.0              httr_1.4.1                 
 [3] bit64_4.0.5                 assertthat_0.2.1           
 [5] stats4_3.6.1                blob_1.2.1                 
 [7] BSgenome_1.52.0             GenomeInfoDbData_1.2.1     
 [9] Rsamtools_2.0.0             yaml_2.2.0                 
[11] progress_1.2.2              pillar_1.6.1               
[13] RSQLite_2.2.7               lattice_0.20-38            
[15] glue_1.4.2                  digest_0.6.20              
[17] GenomicRanges_1.36.0        promises_1.0.1             
[19] XVector_0.24.0              colorspace_1.4-1           
[21] htmltools_0.3.6             httpuv_1.5.1               
[23] Matrix_1.2-18               XML_3.98-1.20              
[25] pkgconfig_2.0.3             biomaRt_2.40.1             
[27] zlibbioc_1.30.0             purrr_0.3.4                
[29] scales_1.1.0                whisker_0.3-2              
[31] later_0.8.0                 BiocParallel_1.18.0        
[33] git2r_0.26.1                tibble_3.1.2               
[35] farver_2.1.0                generics_0.0.2             
[37] IRanges_2.18.1              ellipsis_0.3.2             
[39] withr_2.4.1                 cachem_1.0.5               
[41] SummarizedExperiment_1.14.1 GenomicFeatures_1.36.3     
[43] BiocGenerics_0.30.0         magrittr_2.0.1             
[45] crayon_1.4.1                memoise_2.0.0              
[47] evaluate_0.14               fs_1.3.1                   
[49] fansi_0.5.0                 tools_3.6.1                
[51] data.table_1.14.0           prettyunits_1.0.2          
[53] hms_1.1.0                   lifecycle_1.0.0            
[55] matrixStats_0.57.0          stringr_1.4.0              
[57] S4Vectors_0.22.1            munsell_0.5.0              
[59] DelayedArray_0.10.0         AnnotationDbi_1.46.0       
[61] Biostrings_2.52.0           compiler_3.6.1             
[63] GenomeInfoDb_1.20.0         rlang_0.4.11               
[65] grid_3.6.1                  RCurl_1.98-1.1             
[67] VariantAnnotation_1.30.1    labeling_0.3               
[69] bitops_1.0-6                rmarkdown_1.13             
[71] gtable_0.3.0                DBI_1.1.1                  
[73] R6_2.5.0                    GenomicAlignments_1.20.1   
[75] dplyr_1.0.7                 knitr_1.23                 
[77] rtracklayer_1.44.0          utf8_1.2.1                 
[79] fastmap_1.1.0               bit_4.0.4                  
[81] workflowr_1.6.2             rprojroot_2.0.2            
[83] stringi_1.4.3               parallel_3.6.1             
[85] Rcpp_1.0.6                  vctrs_0.3.8                
[87] tidyselect_1.1.0            xfun_0.8