Blood cells regions (with lower genotype criteria)

workflowr

Ignored files:
    Ignored:    analysis/figure/

Untracked files:
    Untracked:  scripts/bloodcells.steps.md
    Untracked:  scripts/maf.R
    Untracked:  scripts/mafinfo.R

Unstaged changes:
    Modified:   scripts/get_bloodcells_trait_regions.R

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.6, MAF > 0.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_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")
}

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

gwas_HLR_perc = fread('/gpfs/data/stephens-lab/finemap-uk-biobank/data/raw/BloodCells/gwas_maf001_info6/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')

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

gwas_Reticulocyte_perc = fread('/gpfs/data/stephens-lab/finemap-uk-biobank/data/raw/BloodCells/gwas/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')

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_HLR_perc %>% 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()

There are 975 regoins in total, 80 regions with length greater than 2Mb, 70 regions contain more than 10000 SNPs, 292 regions contain more 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')

gwas_Lymphocyte_perc = fread('/gpfs/data/stephens-lab/finemap-uk-biobank/data/raw/BloodCells/gwas_maf001_info6/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')

Session information

sessionInfo()

R version 3.5.0 (2018-04-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server release 6.7 (Santiago)

Matrix products: default
BLAS: /gpfs/apps/haswell/software/gcc-6.2.0/R/3.5.0/lib64/R/lib/libRblas.so
LAPACK: /gpfs/apps/haswell/software/gcc-6.2.0/R/3.5.0/lib64/R/lib/libRlapack.so

locale:
[1] C

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

other attached packages:
[1] knitr_1.21        kableExtra_1.3.1  dplyr_0.8.3       data.table_1.12.0

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.2        highr_0.7         compiler_3.5.0   
 [4] pillar_1.3.1      later_0.7.5       git2r_0.26.1     
 [7] workflowr_1.6.1   tools_3.5.0       digest_0.6.22    
[10] viridisLite_0.3.0 evaluate_0.12     tibble_2.0.1     
[13] pkgconfig_2.0.3   rlang_0.4.1       rstudioapi_0.9.0 
[16] yaml_2.2.0        xfun_0.4          stringr_1.3.1    
[19] httr_1.4.0        xml2_1.2.0        fs_1.3.1         
[22] webshot_0.5.1     rprojroot_1.3-2   tidyselect_0.2.5 
[25] glue_1.3.1        R6_2.4.0          rmarkdown_1.11   
[28] purrr_0.3.0       magrittr_1.5      whisker_0.3-2    
[31] backports_1.1.5   scales_1.0.0      promises_1.0.1   
[34] htmltools_0.3.6   assertthat_0.2.1  rvest_0.3.2      
[37] colorspace_1.4-1  httpuv_1.4.5.1    stringi_1.2.4    
[40] munsell_0.5.0     crayon_1.3.4