R包”gwasrapidd”------快速获取GWAS Catalog数据库的信息

在往期内容中，米老鼠和大家简单介绍过做孟德尔随机化研究使用到的数据库，主要是OpenGWAS， GWAS Catalog 和Phenoscanner这三个。其中，Open GWAS库的数据可以使用“ieugwasr”包来快速获取，具体请参考往期内容。今天我和大家简单介绍一下可以快速获取GWAS Catalog数据库信息的“gwasrapidd”包，该包于近期加入CRAN集。

从获取数据的角度来看，主要使用的有四个函数：get_studies(), get_associations(), get_variants(),和 get_traits()。

1. 使用get_studies()函数

install.packages("gwasrapidd") #安装R包
library(gwasrapidd) #
my_study1 <- get_studies(study_id ='GCST000858') #
slotNames(my_study1) #
#[1] "studies"                 "genotyping_techs"        "platforms"               "ancestries"             
#[5] "ancestral_groups"         "countries_of_origin"      "countries_of_recruitment""publications"
my_study1@studies #查看slot studies的具体信息

返回的结果是一个S4对象，我们可以使用slotNames()函数来获取每个slot的名字，这里我们发现my_studies主要包括8个slot，其中"studies"代表研究的基本信息，"genotyping_techs"代表采取的基因分型技术，"platforms"代表使用的测序平台，"ancestries"和"ancestral_groups"代表人群所属的种族以及不同人种的样本量，"countries_of_origin"和"countries_of_recruitment"代表国别相关信息，"publications"代表发表的文章信息。

这里的get_studies()的主要参数如下：（1）参数study_id：代表GWASCatalog里研究的accession号，前四个字母是“GCST”，可以是向量类型；（2）参数association_id：代表GWAS catalog里的关联信息ID，是一个数字，可以是向量类型；（3）参数variant_id：代表GWASCatalog里的遗传变异信息，一般均为rsid，可以是向量类型；（4）参数efo_id：代表GWAS Catalog里性状的ID号，以“EFO”开头，可以是向量类型；（5）参数pubmed_id：代表研究的PubMedID号；（6）参数full_pvalue_set： 是一个逻辑参数，代表获取是否有完整summary结果的研究（full summary statistics），如果设置为TRUE则代表只 or studies without it (FALSE)；（7）参数efo_trait：是字符串型向量，代表EFO表型描述信息，如“uric acid measurement”。

我们需要注意如果是多参数输入的话，“gwasrapidd”包返回的是多个参数的并集，比如，如果同时输入参数study_id和variant_id，那么返回的结果就是要么包含study_id的，要么包含variant_id的:

my_study2 <- get_studies(study_id ='GCST000858', variant_id = 'rs12752552')
my_study2@studies 

可以看出在这种情况下，返回的结果是或者包含“GCST000858”，或者包含“rs12752552”的研究信息，其等同于如下代码：

s1 <- get_studies(study_id ='GCST000858')
s2 <- get_studies(variant_id ='rs12752552')
my_study2 <- gwasrapidd::union(s1, s2) #这里的union函数来自gwasrapidd包

2. 使用get_associations()函数

my_associations <-get_associations(study_id = my_study1@studies$study_id)
slotNames(my_associations)
#[1] "associations""loci"        "risk_alleles" "genes"        "ensembl_ids"  "entrez_ids" 
as.data.frame(my_associations@associations)

这里get_associations()函数的参数和get_studies()的差不多，单数参数interactive在get_associations()中是比较特殊的，它是一个逻辑型参数，表示是否反应SNP之间的交互作用，默认值为TRUE。最后，该函数会返回6个slot，分别反映关联值大小，位点信息，风险等位基因信息，基因信息，基因的ENSEMBL编码和基因的ENTREZ编码信息，感兴趣的小伙伴可以都是试着查看一下。

3. 使用get_variants()函数

my_variants <- get_variants(study_id =my_study1@studies$study_id)
slotNames(my_variants)
#[1] "variants"         "genomic_contexts""ensembl_ids"     "entrez_ids"     
as.data.frame(my_variants@variants)
# variant_id merged functional_class chromosome_name chromosome_positionchromosome_region    last_update_date
#1 rs7329174      0   intron_variant              13            40983974          13q14.11 2021-10-19 22:22:55
as.data.frame(my_variants@genomic_contexts)

关于get_variants()函数有一个需要注意的参数genomic_range，该参数表示的是指定遗传变异在基因组上的特定位置，它是一个列表型数据，由三组向量构成，分别是染色体号，七点和终点。该函数返回的结果包含4个slot，分别表示遗传变异的信息（不包含GWAS的汇总数据），遗传变异在基因组上的信息，基因的ENSEMBL编码和基因的ENTREZ编码信息。从上图中我们可以看出一个位于内含子区域上的SNP可以对应多个基因，区别就是和不同基因的距离不同，一般我们会选择距离最近的那个基因。

4. 使用get_traits()函数

my_traits <- get_traits(study_id =my_study1@studies$study_id)
slotNames(my_traits)
as.data.frame(my_traits@traits)
#      efo_id                       trait                                 uri
#1 EFO_0002690 systemic lupus erythematosus http://www.ebi.ac.uk/efo/EFO_0002690

函数get_traits()的参数set_operation值得我们关注，它表示对返回的trait的操作，有两个选项，分别是“union”和“intersection”，前者表示取所有的返回的trait，后者表示取交集，默认值是“union”。

虽然“gwasrapidd”包在Github上已存在好几年了，但知道最近才加入到CRAN里，目前来看，该包还是比较稳定和可靠的，有需要的小伙伴赶快学起来吧！