在往期内容中,米老鼠和大家简单介绍过做孟德尔随机化研究使用到的数据库,主要是OpenGWAS, GWAS Catalog 和Phenoscanner这三个。其中,Open GWAS库的数据可以使用“ieugwasr”包来快速获取,具体请参考往期内容。今天我和大家简单介绍一下可以快速获取GWAS Catalog数据库信息的“gwasrapidd”包,该包于近期加入CRAN集。
从获取数据的角度来看,主要使用的有四个函数:get_studies(), get_associations(), get_variants(),和 get_traits()。
1. 使用get_studies()函数
代码语言:javascript复制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的:
代码语言:javascript复制my_study2 <- get_studies(study_id ='GCST000858', variant_id = 'rs12752552')
my_study2@studies
可以看出在这种情况下,返回的结果是或者包含“GCST000858”,或者包含“rs12752552”的研究信息,其等同于如下代码:
代码语言:javascript复制s1 <- get_studies(study_id ='GCST000858')
s2 <- get_studies(variant_id ='rs12752552')
my_study2 <- gwasrapidd::union(s1, s2) #这里的union函数来自gwasrapidd包
2. 使用get_associations()函数
代码语言:javascript复制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()函数
代码语言:javascript复制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()函数
代码语言:javascript复制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里,目前来看,该包还是比较稳定和可靠的,有需要的小伙伴赶快学起来吧!