Last updated: 2020-11-22

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.


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(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.

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 885829c. 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 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)
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:

library(kableExtra)

Attaching package: 'kableExtra'
The following object is masked from 'package:dplyr':

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

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() %>% print
<table class="table" style="margin-left: auto; margin-right: auto;">
<caption>972 regions</caption>
 <thead>
  <tr>
   <th style="text-align:left;">   </th>
   <th style="text-align:right;"> Min. </th>
   <th style="text-align:right;"> 1st Qu. </th>
   <th style="text-align:right;"> Median </th>
   <th style="text-align:right;"> Mean </th>
   <th style="text-align:right;"> 3rd Qu. </th>
   <th style="text-align:right;"> Max. </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> region_length </td>
   <td style="text-align:right;"> 307855.000 </td>
   <td style="text-align:right;"> 500000.00 </td>
   <td style="text-align:right;"> 658060.500 </td>
   <td style="text-align:right;"> 932085.792 </td>
   <td style="text-align:right;"> 1057863.250 </td>
   <td style="text-align:right;"> 8729501 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> region_num_snps </td>
   <td style="text-align:right;"> 21.000 </td>
   <td style="text-align:right;"> 1512.50 </td>
   <td style="text-align:right;"> 2015.000 </td>
   <td style="text-align:right;"> 2641.288 </td>
   <td style="text-align:right;"> 3146.250 </td>
   <td style="text-align:right;"> 21219 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> region_max_log10p </td>
   <td style="text-align:right;"> 7.328 </td>
   <td style="text-align:right;"> 9.94 </td>
   <td style="text-align:right;"> 15.255 </td>
   <td style="text-align:right;"> Inf </td>
   <td style="text-align:right;"> 31.907 </td>
   <td style="text-align:right;"> Inf </td>
  </tr>
</tbody>
</table>

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] kableExtra_1.3.1  dplyr_1.0.2       data.table_1.13.2 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] knitr_1.30        generics_0.1.0    fs_1.5.0          vctrs_0.3.4      
[25] webshot_0.5.2     rprojroot_1.3-2   tidyselect_1.1.0  glue_1.4.2       
[29] R6_2.5.0          rmarkdown_2.5     purrr_0.3.4       magrittr_1.5     
[33] whisker_0.4       scales_1.1.1      backports_1.2.0   promises_1.1.1   
[37] ellipsis_0.3.1    htmltools_0.5.0   rvest_0.3.6       colorspace_1.4-1 
[41] httpuv_1.5.4      stringi_1.5.3     munsell_0.5.0     crayon_1.3.4