Last updated: 2020-11-23

Checks: 7 0

Knit directory: finemap-uk-biobank/

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.

R Markdown file: up-to-date

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.

Environment: empty

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.

Seed: set.seed(20191114)

The command set.seed(20191114) 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.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: 7342408

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 7342408. 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: 


Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    analysis/.Rhistory
    Ignored:    data/.DS_Store
    Ignored:    scripts/.DS_Store

Untracked files:
    Untracked:  data/bloodcells1.csv
    Untracked:  data/height.chr3.matrix
    Untracked:  data/susie_ss_input_sex.rds

Unstaged changes:
    Modified:   analysis/compare_result_more.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/bloodcells_regions.Rmd) and HTML (docs/bloodcells_regions.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 7342408 zouyuxin 2020-11-23 wflow_publish("analysis/bloodcells_regions.Rmd")
html 2bb79f6 zouyuxin 2020-11-22 Build site.
Rmd e3da525 zouyuxin 2020-11-22 wflow_publish("analysis/bloodcells_regions.Rmd")
html 61d98d4 zouyuxin 2020-11-22 Build site.
Rmd 885829c zouyuxin 2020-11-22 wflow_publish("analysis/bloodcells_regions.Rmd")
html 9ab9598 zouyuxin 2020-11-22 Build site.
Rmd 8b1d17e zouyuxin 2020-11-22 wflow_publish("analysis/bloodcells_regions.Rmd")
html 2b09b42 zouyuxin 2020-11-20 Build site.
Rmd 4044ce5 zouyuxin 2020-11-20 wflow_publish("analysis/bloodcells_regions.Rmd")
html 50a6e41 zouyuxin 2020-11-20 Build site.
Rmd 3bfe683 zouyuxin 2020-11-20 wflow_publish("analysis/bloodcells_regions.Rmd")
html 717d6b1 zouyuxin 2020-11-20 Build site.
Rmd c2021d0 zouyuxin 2020-11-20 wflow_publish("analysis/bloodcells_regions.Rmd")

There are 248,980 individuals of white British ancestries with 16 blood cells phenotypes. The script to prepare the phenotypes and covariates is get_bloodcells. The filtering steps are also described here

For genotype data, variants with imputation score (INFO) > 0.9, MAF > 1% are included in association studies.

The script to run GWAS is GWAS

For each phenotype, regions for fine-mapping are defined by greedily starting with the most significantly associated SNP, including SNPs within a window of 500kb centered at the SNP, until we include all significant SNPs (p < 5e-8). We merge ovelapping regions. We exclude HLA region (chr6: 25Mb - 36Mb). The steps are

  1. Find the most significantly associated SNP.

  2. Choose region +- 250kb around the SNP.

  3. Find the next most significantly associated SNP ouside the selected regions.

  4. Choose region +- 250kb around the SNP.

  5. Merge regions if they overlap. ...

When we select region across traits, we include all regions from each pheotype and merge overlapping regions. This produces some very large regions with more than 10000 SNPs.

library(data.table)
library(dplyr)
library(kableExtra)
library(knitr)
pheno_names = c("WBC_count", "RBC_count", "Haemoglobin", "MCV", "RDW", "Platelet_count",
                "Plateletcrit", "PDW", "Lymphocyte_perc", "Monocyte_perc",
                "Neutrophill_perc", "Eosinophill_perc", "Basophill_perc",
                "Reticulocyte_perc", "MSCV", "HLR_perc")
trait_regions = list()
for(trait in pheno_names){
  # region = fread(paste0('/gpfs/data/stephens-lab/finemap-uk-biobank/data/raw/BloodCells/regions/', trait, '_regions'))
  region = fread(paste0('~/Desktop/ukb-bloodcells/regions_raw/', trait, '_regions'))
  region = region %>% arrange(desc(logp))
  region_r = c()
  for(i in 1:22){
    region.chr = region %>% filter(CHR == i) %>% arrange(start)
    if(nrow(region.chr) == 0){
      next
    }
    tmp = region.chr %>% group_by(g = cumsum(cummax(lag(end, default = first(end))) < start)) %>%
      summarise(start = first(start), end = max(end), logp = max(logp),.groups = 'drop') %>% 
      mutate(length = end - start) %>%
      mutate(CHR = i) %>% select(CHR, start, end, length, logp)
    region_r = rbind(region_r, tmp)
  }
  trait_regions[[trait]] = region_r
}

Summary of region length for each phenotype before selecting regions across traits:

for(trait in pheno_names){
  tb = rbind(summary(trait_regions[[trait]]$length), summary(trait_regions[[trait]]$logp))
  rownames(tb) = c('region_length', 'region_max_log10p')
  tb = round(tb,3) %>% kbl(caption = paste0(trait, ': ', nrow(trait_regions[[trait]]), ' regions')) %>% kable_styling()
  cat(tb)
  cat("\n")
}
WBC_count: 280 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 378854.000 5.000e+05 500000.000 610315.389 500000.000 2301820.000
region_max_log10p 7.327 8.531e+00 10.807 15.286 14.601 213.266
RBC_count: 310 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 5.000e+05 5.000e+05 500000.000 653220.6 755004.000 2940124
region_max_log10p 7.307e+00 8.779e+00 11.363 Inf 18.728 Inf
Haemoglobin: 250 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 301703.000 5.000e+05 500000.000 611650.000 500000.000 3251867.000
region_max_log10p 7.302 8.249e+00 10.408 17.123 16.425 261.853
MCV: 326 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 5.000e+05 5.000e+05 500000.000 670089.4 772562.750 2920824
region_max_log10p 7.312e+00 9.576e+00 13.151 Inf 21.175 Inf
RDW: 275 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 275738.000 5.000e+05 500000.000 650422 500000.00 4089364
region_max_log10p 7.309 9.108e+00 12.833 Inf 24.55 Inf
Platelet_count: 397 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 270883.000 5.000e+05 500000.000 678599.6 766174.000 3646000
region_max_log10p 7.324 9.018e+00 12.161 Inf 21.916 Inf
Plateletcrit: 382 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 270883.000 5.000e+05 500000.000 640938.8 754987.500 3593140
region_max_log10p 7.317 8.826e+00 11.431 Inf 19.293 Inf
PDW: 299 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 268602.000 5.000e+05 500000.000 650832 751519.000 4239970
region_max_log10p 7.318 8.975e+00 12.738 Inf 21.975 Inf
Lymphocyte_perc: 232 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 378854.000 5.000e+05 500000.000 621718.961 753623.00 2239954.000
region_max_log10p 7.325 9.048e+00 11.351 16.049 16.54 152.602
Monocyte_perc: 262 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 310142.000 5.00e+05 500000.000 659465.2 751945.000 3864804
region_max_log10p 7.361 8.85e+00 12.112 Inf 21.705 Inf
Neutrophill_perc: 218 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 368127.000 5.000e+05 500000.000 612052.165 500000.000 2208916.000
region_max_log10p 7.303 8.403e+00 10.917 15.978 16.034 193.111
Eosinophill_perc: 277 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 397874.000 5.000e+05 500000.000 656940.123 750663.000 3994646.000
region_max_log10p 7.328 9.148e+00 13.378 20.676 20.883 212.968
Basophill_perc: 85 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 5.000e+05 5.000e+05 500000.000 561068.212 500000.000 1748896.000
region_max_log10p 7.367e+00 8.567e+00 11.392 17.352 16.739 134.511
Reticulocyte_perc: 237 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 303607.000 5.000e+05 500000.000 653925.3 753414.000 5547602
region_max_log10p 7.302 8.534e+00 12.249 Inf 20.093 Inf
MSCV: 285 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 272793.000 5.000e+05 500000.000 656840.8 754786.000 3254290
region_max_log10p 7.328 9.262e+00 12.434 Inf 21.282 Inf
HLR_perc: 246 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 381586.000 5.000e+05 500000.000 674728.8 768513.250 5973226
region_max_log10p 7.304 8.853e+00 12.103 Inf 19.311 Inf

For HLR_perc, the maximum region is at CHR 3 from 46234573 to 52207799, which includes 9572 SNPs.

# gwas_HLR_perc = fread('/gpfs/data/stephens-lab/finemap-uk-biobank/data/raw/BloodCells/gwas/bloodcells_gwas_HLR_perc')
gwas_HLR_perc = fread('~/Desktop/ukb-bloodcells/bloodcells_gwas_HLR_perc')
colnames(gwas_HLR_perc)[1] = 'CHR'
gwas_HLR_perc$P = as.numeric(gwas_HLR_perc$P)
gwas_HLR_perc = gwas_HLR_perc %>% select(CHR, POS, T_STAT, P) %>% mutate(logp = -log10(P))
gwas_HLR_perc.sub = gwas_HLR_perc %>% filter(CHR == 3, POS >= 46234573, POS <=52207799)
plot(gwas_HLR_perc.sub$POS, gwas_HLR_perc.sub$logp, xlab='CHR 3 POS', ylab='-log10(p)', main='HLR_perc')

Version Author Date
50a6e41 zouyuxin 2020-11-20
717d6b1 zouyuxin 2020-11-20

For Reticulocyte_perc, the maximum region is at CHR 3 from 48155661 to 53703263, which includes 8888 SNPs.

# gwas_Reticulocyte_perc = fread('/gpfs/data/stephens-lab/finemap-uk-biobank/data/raw/BloodCells/gwas/bloodcells_gwas_Reticulocyte_perc')
gwas_Reticulocyte_perc = fread('~/Desktop/ukb-bloodcells/bloodcells_gwas_Reticulocyte_perc')
colnames(gwas_Reticulocyte_perc)[1] = 'CHR'
gwas_Reticulocyte_perc$P = as.numeric(gwas_Reticulocyte_perc$P)
gwas_Reticulocyte_perc = gwas_Reticulocyte_perc %>% select(CHR, POS, T_STAT, P) %>% mutate(logp = -log10(P))
gwas_Reticulocyte_perc.sub = gwas_Reticulocyte_perc %>% filter(CHR == 3, POS >= 48155661, POS <=53703263)
plot(gwas_Reticulocyte_perc.sub$POS, gwas_Reticulocyte_perc.sub$logp, xlab='CHR 3 POS', ylab='-log10(p)', main='Reticulocyte_perc')

Version Author Date
50a6e41 zouyuxin 2020-11-20
717d6b1 zouyuxin 2020-11-20

For PDW, the maximum region is at CHR 8 from 7838230 to 12078200, which includes 16605 SNPs.

# gwas_PDW = fread('/gpfs/data/stephens-lab/finemap-uk-biobank/data/raw/BloodCells/gwas/bloodcells_gwas_PDW')
gwas_PDW = fread('~/Desktop/ukb-bloodcells/bloodcells_gwas_PDW')
colnames(gwas_PDW)[1] = 'CHR'
gwas_PDW$P = as.numeric(gwas_PDW$P)
gwas_PDW = gwas_PDW %>% select(CHR, POS, T_STAT, P) %>% mutate(logp = -log10(P))
gwas_PDW.sub = gwas_PDW %>% filter(CHR == 8, POS >= 7838230, POS <=12078200)
plot(gwas_PDW.sub$POS, gwas_PDW.sub$logp, xlab='CHR 3 POS', ylab='-log10(p)', main='PDW')

Version Author Date
50a6e41 zouyuxin 2020-11-20
717d6b1 zouyuxin 2020-11-20

Select regions across phenotype:

tb = bind_rows(trait_regions, .id = "column_label")
res.final = c()
for(i in 1:22){
  tb.chr = tb %>% filter(CHR == i) %>% arrange(start)
  if(nrow(tb.chr) == 0){
    next
  }
  tmp = tb.chr %>% group_by(g = cumsum(cummax(lag(end, default = first(end))) < start)) %>%
    summarise(start = first(start), end = max(end), logp = max(logp), .groups = 'drop') %>% 
    mutate(length = end - start) %>%
    mutate(CHR = i) %>% select(CHR, start, end, length, logp)
  res.final = rbind(res.final, tmp)
}
snpsnum = c()
for(i in 1:nrow(res.final)){
  snpsnum = c(snpsnum, gwas_PDW %>% filter(CHR == res.final$CHR[i],
                                           POS >= res.final$start[i], 
                                           POS <= res.final$end[i]) %>% nrow )
}
res.final$snpsnum = snpsnum

Summary of regions:

tb = rbind(summary(res.final$length), summary(res.final$snpsnum),  summary(res.final$logp))
rownames(tb) = c('region_length', 'region_num_snps', 'region_max_log10p')
round(tb,3)  %>% kbl(caption = paste0(nrow(res.final), ' regions')) %>% kable_styling()
972 regions
Min. 1st Qu. Median Mean 3rd Qu. Max.
region_length 307855.000 500000.00 658060.500 932085.792 1057863.250 8729501
region_num_snps 21.000 1512.50 2015.000 2641.288 3146.250 21219
region_max_log10p 7.328 9.94 15.255 Inf 31.907 Inf

There are 972 regoins in total, 68 regions with length greater than 2Mb, 84 regions contain greater than 5000 SNPs.

The region with maximum length and maximum number of SNPs:

gwas_HLR_perc.max = gwas_HLR_perc %>% filter(CHR == 17, POS >= 39754910, POS <=48484411)
plot(gwas_HLR_perc.max$POS, gwas_HLR_perc.max$logp, xlab='CHR 17 POS', ylab='-log10(p)', main='HLR_perc')

Version Author Date
9ab9598 zouyuxin 2020-11-22 
gwas_Reticulocyte_perc.max = gwas_Reticulocyte_perc %>% filter(CHR == 17, POS >= 39754910, POS <=48484411)
plot(gwas_Reticulocyte_perc.max$POS, gwas_Reticulocyte_perc.max$logp, xlab='CHR 17 POS', ylab='-log10(p)', main='Reticulocyte_perc')

Version Author Date
9ab9598 zouyuxin 2020-11-22 
gwas_PDW.max = gwas_PDW %>% filter(CHR == 17, POS >= 39754910, POS <=48484411)
plot(gwas_PDW.max$POS, gwas_PDW.max$logp, xlab='CHR 17 POS', ylab='-log10(p)', main='PDW')

Version Author Date
9ab9598 zouyuxin 2020-11-22 
# gwas_Lymphocyte_perc = fread('/gpfs/data/stephens-lab/finemap-uk-biobank/data/raw/BloodCells/gwas/bloodcells_gwas_Lymphocyte_perc')
gwas_Lymphocyte_perc = fread('~/Desktop/ukb-bloodcells/bloodcells_gwas_Lymphocyte_perc')
colnames(gwas_Lymphocyte_perc)[1] = 'CHR'
gwas_Lymphocyte_perc$P = as.numeric(gwas_Lymphocyte_perc$P)
gwas_Lymphocyte_perc = gwas_Lymphocyte_perc %>% select(CHR, POS, T_STAT, P) %>% mutate(logp = -log10(P))
gwas_Lymphocyte_perc.max = gwas_Lymphocyte_perc %>% filter(CHR == 17, POS >= 39754910, POS <=48484411)
plot(gwas_Lymphocyte_perc.max$POS, gwas_Lymphocyte_perc.max$logp, xlab='CHR 17 POS', ylab='-log10(p)', main='Lymphocyte_perc')

Version Author Date
9ab9598 zouyuxin 2020-11-22 

sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.16

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] knitr_1.30        kableExtra_1.3.1  dplyr_1.0.2       data.table_1.13.2
[5] workflowr_1.6.2  

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.5        highr_0.8         pillar_1.4.6      compiler_3.6.3   
 [5] later_1.1.0.1     git2r_0.27.1      tools_3.6.3       digest_0.6.27    
 [9] viridisLite_0.3.0 evaluate_0.14     lifecycle_0.2.0   tibble_3.0.4     
[13] pkgconfig_2.0.3   rlang_0.4.8       rstudioapi_0.11   yaml_2.2.1       
[17] xfun_0.19         xml2_1.3.2        stringr_1.4.0     httr_1.4.2       
[21] generics_0.1.0    fs_1.5.0          vctrs_0.3.4       webshot_0.5.2    
[25] rprojroot_1.3-2   tidyselect_1.1.0  glue_1.4.2        R6_2.5.0         
[29] rmarkdown_2.5     purrr_0.3.4       magrittr_1.5      whisker_0.4      
[33] scales_1.1.1      backports_1.2.0   promises_1.1.1    ellipsis_0.3.1   
[37] htmltools_0.5.0   rvest_0.3.6       colorspace_1.4-1  httpuv_1.5.4     
[41] stringi_1.5.3     munsell_0.5.0     crayon_1.3.4