Using haploR, an R package for querying HaploReg and RegulomeDB

Arguments

• query: Query (a vector of rsIDs).
• file: A text file (one refSNP ID per line).
• study: A particular study. See function getHaploRegStudyList(…). Default: NULL.
• ldThresh: LD threshold, r2 (select NA to show only query variants). Default: 0.8.
• ldPop: 1000G Phase 1 population for LD calculation. Can be: AFR (Africa), AMR (America), ASN (Asia). Default: EUR (Europe).
• epi: Source for epigenomes. Possible values: vanilla for ChromHMM (Core 15-state model); imputed for ChromHMM (25-state model using 12 imputed marks); methyl for H3K4me1/H3K4me3 peaks; acetyl for H3K27ac/H3K9ac peaks. Default: vanilla.
• cons: Mammalian conservation algorithm. Possible values: gerp for GERP (http://mendel.stanford.edu/SidowLab/downloads/gerp/), siphy for SiPhy-omega, both for both. Default: siphy.
• genetypes: Show position relative to. Possible values: gencode for Gencode genes; refseq for RefSeq genes; both for both. Default: gencode.
• url: HaploReg url address. Default: https://pubs.broadinstitute.org/mammals/haploreg/haploreg.php
• timeout: A timeout parameter for curl. Default: 10
• encoding: Set the encoding for correct retrieval web-page content. Default: UTF-8
• verbose: Verbosing output. Default: FALSE.

Output

A data frame (table) wrapped into a tibble object contains data extracted from HaploReg site. The colums (33 in total at the time of writing this vignette) are specific to HaploReg output. Below we describe the columns:

• chr: Chromosome, type: numeric
• pos_hg38: Position on the human genome, type: numeric.
• r2: Linkage disequilibrium. Type: numeric.
• D’: Linkage disequilibrium, alternative definition. Type: numeric.
• is_query_snp: Indicator shows query SNP, 0 - not query SNP, 1 - query SNP. Type: numeric.
• rsID: refSNP ID. Type: character.
• ref: Reference allele. Type: character.
• alt: Alternative allele. Type: character.
• AFR: r2 calculated for Africa. Type: numeric.
• AMR: r2 calculated for America. Type: numeric.
• ASN: r2 calculated for Asia. Type: numeric.
• EUR: r2 calculated for Europe. Type: numeric.
• GERP_cons: GERP scores. Type: numeric.
• SiPhy_cons: SiPhy scores. Type: numeric.
• Chromatin_States: Chromatin states: reference epigenome identifiers (EID) of chromatin-associated proteins and histone modifications in that region. Type: character.
• Chromatin_States_Imputed: Chromatin states based on imputed data. Type: character.
• Chromatin_Marks: Chromatin marks Type: character.
• DNAse: DNAse. Type: chararcter.
• Proteins: A list of protein names. Type: character.
• eQTL: Expression Quantitative Trait Loci. Type: character.
• gwas: GWAS study name. Type: character.
• grasp: GRASP study name: character.
• Motifs: Motif names. Type: character.
• GENCODE_id: GENCODE transcript ID. Type: character.
• GENCODE_name: GENCODE gene name. Type: character.
• GENCODE_direction: GENCODE direction (transcription toward 3’ or 5’ end of the DNA sequence). Type: numeric.
• GENCODE_distance: GENCODE distance. Type: numeric.
• RefSeq_id: NCBI Reference Sequence Accession number. Type: character.
• RefSeq_name: NCBI Reference Sequence name. Type: character.
• RefSeq_direction: NCBI Reference Sequence direction (transcription toward 3’ or 5’ end of the DNA sequence). Type: numeric.
• RefSeq_distance: NCBI Reference Sequence distance. Type: numeric.
• dbSNP_functional_annotation Annotated proteins associated with the SNP. Type: numeric.
• query_snp_rsid: Query SNP rs ID. Type: character.
• Promoter_histone_marks: Promoter histone marks. Type: factor.
• Enhancer_histone_marks: Enhancer histone marks. Type: factor.

Number of rows is not constant, at least equal or more than the number of query SNPs, and depends on r2 parameter choosen in a query (default 0.8). This means that the program outputs not only query SNPs, but also those SNPs that have r2 ≥ 0.8 with the query SNPs.

One or several genetic variants

library(haploR)
x <- queryHaploreg(query=c("rs10048158","rs4791078"))
x

## # A tibble: 33 × 35
##    chr     pos_hg38    r2  `D'` is_query_snp rsID  ref   alt     AFR   AMR   ASN
##    <chr>   <chr>    <dbl> <dbl>        <dbl> <chr> <chr> <chr> <dbl> <dbl> <dbl>
##  1 17      66240200  1     1               1 rs10… T     C      0.93  0.66  0.92
##  2 17      66240200  0.82  0.93            0 rs10… T     C      0.93  0.66  0.92
##  3 Array17 66231972  0.82  0.99            0 rs11… G     T      0.88  0.63  0.9 
##  4 Array17 66248387  0.99  1               0 rs12… T     C      0.93  0.66  0.92
##  5 Array17 66248387  0.81  0.92            0 rs12… T     C      0.93  0.66  0.92
##  6 17      66214285  0.98  0.99            0 rs19… G     C      0.86  0.67  0.91
##  7 17      66214285  0.83  0.93            0 rs19… G     C      0.86  0.67  0.91
##  8 Array17 66240359  0.83  1               0 rs20… A     G      0.64  0.6   0.91
##  9 Array17 66219777  0.82  0.99            0 rs22… A     G      0.89  0.63  0.89
## 10 17      66219453  0.83  0.93            0 rs22… G     A      0.91  0.67  0.91
## # ℹ 23 more rows
## # ℹ 24 more variables: EUR <dbl>, GERP_cons <chr>, SiPhy_cons <chr>,
## #   Chromatin_States <chr>, Chromatin_States_Imputed <chr>,
## #   Chromatin_Marks <chr>, DNAse <chr>, Proteins <chr>, eQTL <chr>, gwas <chr>,
## #   grasp <chr>, Motifs <chr>, GENCODE_id <chr>, GENCODE_name <chr>,
## #   GENCODE_direction <chr>, GENCODE_distance <chr>, RefSeq_id <chr>,
## #   RefSeq_name <chr>, RefSeq_direction <chr>, RefSeq_distance <chr>, …

Here query is a vector with names of genetic variants.

We then can create a subset from the results, for example, to choose only SNPs with r2 > 0.9:

if(length(x)!=0) {
  subset.high.LD <- x[as.numeric(x$r2) > 0.9, c("rsID", "r2", "chr", "pos_hg38", "is_query_snp", "ref", "alt")]
  subset.high.LD
}

## # A tibble: 13 × 7
##    rsID          r2 chr     pos_hg38 is_query_snp ref   alt  
##    <chr>      <dbl> <chr>   <chr>           <dbl> <chr> <chr>
##  1 rs10048158  1    17      66240200            1 T     C    
##  2 rs12603947  0.99 Array17 66248387            0 T     C    
##  3 rs1971682   0.98 17      66214285            0 G     C    
##  4 rs2215415   0.99 17      66219453            0 G     A    
##  5 rs3744317   0.99 Array17 66220526            0 G     A    
##  6 rs4366742   0.99 Array17 66216124            0 T     C    
##  7 rs4790914   1    17      66213160            0 C     G    
##  8 rs4791078   1    17      66213896            1 A     C    
##  9 rs4791079   1    17      66213422            0 T     G    
## 10 rs71160546  0.94 Array17 66230111            0 GA    G    
## 11 rs7342920   0.99 Array17 66248527            0 T     G    
## 12 rs8178827   0.96 Array17 66227121            0 C     T    
## 13 rs9895261   1    17      66244318            0 A     G

We can then save the subset.high.LD into an Excel workbook:

require(openxlsx)
write.xlsx(x=subset.high.LD, file="subset.high.LD.xlsx")

This was an example of gathering post-gwas information directly from the online tool. haploR has an additional advantage because it downloads the full information for query retrieved by HaploReg. For example, if you go online and submit these two SNPs to HaploReg (https://pubs.broadinstitute.org/mammals/haploreg/haploreg.php), you will see that some cells of columns “Motifs changed” and “Selected eQTL hits” are hidded (only number of hits are given). haploR retrives this information in a form of a data frame which can be saved into Excel file.

if(length(x)!=0) {
    x[, c("Motifs", "rsID")]
    x[, c("eQTL", "rsID")]
}

Uploading file with variants

If you have a file with your SNPs you would like to analyze, you can supply it as an input as follows:

library(haploR)
x <- queryHaploreg(file=system.file("extdata/snps.txt", package = "haploR"))
x

File “snps.txt” is a text file which contains one rs-ID per line:

rs10048158
rs4791078

Using existing studies

Sometimes one would like to explore results from already performed study. In this case you should first explore the existing studies from HaploReg web site (https://pubs.broadinstitute.org/mammals/haploreg/haploreg.php) and then use one of them as an input parameter. See example below:

library(haploR)
# Getting a list of existing studies:
studies <- getStudyList()
# Let us look at the first element:
if(!is.null(studies)) {
    studies[[1]]
    # Let us look at the second element:
    studies[[2]]
    # Query Hploreg to explore results from 
    # this study:
    x <- queryHaploreg(study=studies[[1]])
    x
}

Extended view of SNP

If you want to see more information about SNP of interest, you can use the getExtendedView: Parameters: snp - A SNP of interest. url - A url to the HaploReg. Default: “https://pubs.broadinstitute.org/mammals/haploreg/detail_v4.1.php?query=&id=”

Return A list of tables t1, t2, …, etc, depending on information contained in HaploReg database.

Example:

library(haploR)
tables <- getExtendedView(snp="rs10048158")
tables

Arguments

• query: Query (a vector of rsIDs).
• genomeAssembly : Genome assembly built: can be GRCh37 or GRCh38
• limit : It controls how many variants will be queried and returned for a large region. It can be a number (1000 by default) or “all”. Please note that large number or “all” may get yourself hurt because you could get timeout or may even crash the server.
• timeout: A ‘timeout’ parameter for ‘curl’. Default: 100

Output

A data frame (table) OR a list with the following items: - guery_coordinates - features - regulome_score - variants - nearby_snps - assembly

Example

library(haploR)

# With RsIDs:
x <- queryRegulome(c("rs4791078","rs10048158"))
x

# With region:
y <- queryRegulome("chr1:39492461-39492462")
y