## Warning: replacing previous import 'utils::findMatches' by
## 'S4Vectors::findMatches' when loading 'AnnotationDbi'
NHGRI maintains and routinely updates a database of selected genome-wide association studies. This document describes R/Bioconductor facilities for working with contents of this database.
Once the package has been installed, use to obtain interactive access to all the facilities. After executing this command, use to obtain an overview. The current version of this vignette can always be accessed at www.bioconductor.org, or by suitably navigating the web pages generated with .
## # A tibble: 447,939 × 38
## DATE ADDED TO CATALO…¹ PUBMEDID `FIRST AUTHOR` DATE JOURNAL LINK STUDY
## <date> <dbl> <chr> <date> <chr> <chr> <chr>
## 1 2008-06-16 15761122 Klein RJ 2005-03-10 Science www.… Comp…
## 2 2008-06-16 16252231 Maraganore DM 2005-09-09 Am J H… www.… High…
## 3 2008-06-16 16648850 Arking DE 2006-04-30 Nat Ge… www.… A co…
## 4 2008-06-16 17052657 Fung HC 2006-09-28 Lancet… www.… Geno…
## 5 2008-06-16 17052657 Fung HC 2006-09-28 Lancet… www.… Geno…
## 6 2008-06-16 17052657 Fung HC 2006-09-28 Lancet… www.… Geno…
## 7 2008-06-16 17053108 Dewan A 2006-10-19 Science www.… HTRA…
## 8 2008-06-16 17068223 Duerr RH 2006-10-26 Science www.… A ge…
## 9 2008-06-16 17158188 Bierut LJ 2006-12-07 Hum Mo… www.… Nove…
## 10 2008-06-16 17158188 Bierut LJ 2006-12-07 Hum Mo… www.… Nove…
## # ℹ 447,929 more rows
## # ℹ abbreviated name: ¹`DATE ADDED TO CATALOG`
## # ℹ 31 more variables: `DISEASE/TRAIT` <chr>, `INITIAL SAMPLE SIZE` <chr>,
## # `REPLICATION SAMPLE SIZE` <chr>, REGION <chr>, CHR_ID <chr>, CHR_POS <dbl>,
## # `REPORTED GENE(S)` <chr>, MAPPED_GENE <chr>, UPSTREAM_GENE_ID <chr>,
## # DOWNSTREAM_GENE_ID <chr>, SNP_GENE_IDS <chr>, UPSTREAM_GENE_DISTANCE <dbl>,
## # DOWNSTREAM_GENE_DISTANCE <dbl>, `STRONGEST SNP-RISK ALLELE` <chr>, …
We can produce a GRanges in two forms. By default
we get an mcols that has a small set of columns. Note
that records that lack a CHR_POS
value are omitted.
Records that have complicated CHR_POS
values, including
semicolons or " x " notation are kept, but only the first
position is retained. The CHR_ID
field may have
complicated character values, these are not normalized,
and are simply used as seqnames
“as is”.
## dropping 20963 records that have NA for CHR_POS
## GRanges object with 426976 ranges and 4 metadata columns:
## seqnames ranges strand | PUBMEDID DATE
## <Rle> <IRanges> <Rle> | <numeric> <Date>
## [1] 6 31600692 * | 25194280 2014-09-07
## [2] 11 61803910 * | 25194280 2014-09-07
## [3] 6 92812416 * | 25194280 2014-09-07
## [4] 5 83523300 * | 25194280 2014-09-07
## [5] 16 9344835 * | 25194280 2014-09-07
## ... ... ... ... . ... ...
## [426972] 11 61856709 * | 35393526 2022-04-07
## [426973] 1 242773726 * | 35393526 2022-04-07
## [426974] 11 61834531 * | 35393526 2022-04-07
## [426975] 7 72673613 * | 35393526 2022-04-07
## [426976] 11 61843278 * | 35393526 2022-04-07
## DISEASE/TRAIT SNPS
## <character> <character>
## [1] Laryngeal squamous c.. rs2857595
## [2] Laryngeal squamous c.. rs174549
## [3] Laryngeal squamous c.. rs9445023
## [4] Laryngeal squamous c.. rs310518
## [5] Laryngeal squamous c.. rs40129
## ... ... ...
## [426972] Triglycerides (58:10.. rs97384
## [426973] Triglycerides (58:10.. rs116817123
## [426974] Triglycerides (58:10.. rs174575
## [426975] Triglycerides (58:5).. rs6460025
## [426976] Triglycerides (58:9).. rs174584
## -------
## seqinfo: 24 sequences from GRCh38 genome; no seqlengths
We can set the seqinfo as follows, retaining only records that employ standard chromosomes.
library(GenomeInfoDb)
data(si.hs.38)
gg = keepStandardChromosomes(gg, pruning="coarse")
seqlevels(gg) = seqlevels(si.hs.38)
seqinfo(gg) = si.hs.38
gg
## GRanges object with 426976 ranges and 4 metadata columns:
## seqnames ranges strand | PUBMEDID DATE
## <Rle> <IRanges> <Rle> | <numeric> <Date>
## [1] 6 31600692 * | 25194280 2014-09-07
## [2] 11 61803910 * | 25194280 2014-09-07
## [3] 6 92812416 * | 25194280 2014-09-07
## [4] 5 83523300 * | 25194280 2014-09-07
## [5] 16 9344835 * | 25194280 2014-09-07
## ... ... ... ... . ... ...
## [426972] 11 61856709 * | 35393526 2022-04-07
## [426973] 1 242773726 * | 35393526 2022-04-07
## [426974] 11 61834531 * | 35393526 2022-04-07
## [426975] 7 72673613 * | 35393526 2022-04-07
## [426976] 11 61843278 * | 35393526 2022-04-07
## DISEASE/TRAIT SNPS
## <character> <character>
## [1] Laryngeal squamous c.. rs2857595
## [2] Laryngeal squamous c.. rs174549
## [3] Laryngeal squamous c.. rs9445023
## [4] Laryngeal squamous c.. rs310518
## [5] Laryngeal squamous c.. rs40129
## ... ... ...
## [426972] Triglycerides (58:10.. rs97384
## [426973] Triglycerides (58:10.. rs116817123
## [426974] Triglycerides (58:10.. rs174575
## [426975] Triglycerides (58:5).. rs6460025
## [426976] Triglycerides (58:9).. rs174584
## -------
## seqinfo: 24 sequences from GRCh38 genome
We use BiocFileCache to manage downloaded TSV from EBI’s download site. The file is provided without compression, so prepare for 200+MB download if you are not working from a cache. There is no etag set, so you have to check for updates on your own.
## function (url = "http://www.ebi.ac.uk/gwas/api/search/downloads/alternative",
## cache = BiocFileCache::BiocFileCache(), refresh = FALSE,
## ...)
## NULL
This is converted to a manageable extension of GRanges using
process_gwas_dataframe
.
## function (df)
## NULL
Available functions are:
## gwascat loaded. Use makeCurrentGwascat() to extract current image.
## from EBI. The data folder of this package has some legacy extracts.
## [1] "as_GRanges" "bindcadd_snv" "getRsids"
## [4] "getTraits" "get_cached_gwascat" "gwascat_from_AHub"
## [7] "gwcat_snapshot" "gwcex2gviz" "ldtagr"
## [10] "locs4trait" "makeCurrentGwascat" "obo2graphNEL"
## [13] "process_gwas_dataframe" "riskyAlleleCount" "subsetByChromosome"
## [16] "subsetByTraits" "topTraits" "traitsManh"
An extended GRanges instance with a sample of 50000 SNP-disease associations reported
on 30 April 2020 is
obtained as follows, with addresses based on the GRCh38 genome build.
We use gwtrunc
to refer to it in the sequel.
To determine the most frequently occurring traits in this sample:
##
## Blood protein levels
## 1941
## Heel bone mineral density
## 1309
## Body mass index
## 1283
## Height
## 1238
## Metabolite levels
## 691
## Systolic blood pressure
## 654
## Schizophrenia
## 643
## Educational attainment (years of education)
## 642
## Post bronchodilator FEV1/FVC ratio
## 479
## Type 2 diabetes
## 466
For a given trait, obtain a GRanges with all recorded associations; here only three associations are shown:
## gwasloc instance with 3 records and 38 attributes per record.
## Extracted: 2020-04-30 23:24:51
## metadata()$badpos includes records for which no unique locus was given.
## Genome: GRCh38
## Excerpt:
## GRanges object with 3 ranges and 3 metadata columns:
## seqnames ranges strand | DISEASE/TRAIT SNPS P-VALUE
## <Rle> <IRanges> <Rle> | <character> <character> <numeric>
## [1] 7 21567734 * | LDL cholesterol rs12670798 6e-09
## [2] 5 75355259 * | LDL cholesterol rs3846662 2e-11
## [3] 2 43837951 * | LDL cholesterol rs6756629 3e-10
## -------
## seqinfo: 24 sequences from GRCh38 genome
A basic Manhattan plot is easily constructed with the ggbio package facilities. Here we confine attention to chromosomes 4:6. First, we create a version of the catalog with \(-log_{10} p\) truncated at a maximum value of 25.
requireNamespace("S4Vectors")
mcols = S4Vectors::mcols
mlpv = mcols(gwtrunc)$PVALUE_MLOG
mlpv = ifelse(mlpv > 25, 25, mlpv)
S4Vectors::mcols(gwtrunc)$PVALUE_MLOG = mlpv
library(GenomeInfoDb)
# seqlevelsStyle(gwtrunc) = "UCSC" # no more!
seqlevels(gwtrunc) = paste0("chr", seqlevels(gwtrunc))
gwlit = gwtrunc[ which(as.character(seqnames(gwtrunc)) %in% c("chr4", "chr5", "chr6")) ]
library(ggbio)
mlpv = mcols(gwlit)$PVALUE_MLOG
mlpv = ifelse(mlpv > 25, 25, mlpv)
S4Vectors::mcols(gwlit)$PVALUE_MLOG = mlpv
## Warning: The `facets` argument of `facet_grid()` is deprecated as of ggplot2 2.2.0.
## ℹ Please use the `rows` argument instead.
## ℹ The deprecated feature was likely used in the ggbio package.
## Please report the issue at <https://github.com/lawremi/ggbio/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.