PheWAS(全表型组关联分析)----PheWAS的R包测试及分析(二)

由于上次写完教程一后还有许多疑问还没解决,在师兄的建议下看了一下R包中核心函数的源码(到底在做什么回归?),大概解开了心中的部分疑惑。最近在看很多问题时候,总会有些急躁,反而忽略了很多核心的东西,还是要好好的静下来反思一下的。那么如何做PheWAS呢?下面我先以R包测试数据为例进行介绍,系统分析如何借助这个包跑自己的数据。

1.R包下载

#相关R包的安装(GitHub上由)
install.packages("devtools")
install.packages(c("dplyr","tidyr","ggplot2","MASS","meta","ggrepel","DT"))
devtools::install_github("PheWAS/PheWAS")

2.R测试(看看,我的重点在第三大点)

vignette("PheWAS-package")  #打开该R包的帮助文档
library(PheWAS)                        #加载包
set.seed(1)                                 #设置随机种子
ex=generateExample()             #产生测试数据
#提取list中的三个部分(这三个文件的解释我后面会讲)
id.vocab.code.count=ex$id.vocab.code.count
genotypes=ex$genotypes        #基因型文件
id.sex=ex$id.sex                        #性别,作为协变量
#创建phecode表-翻译代码
phenotypes=createPhenotypes(id.vocab.code.count, aggregate.fun=sum, id.sex=id.sex)
#运行PheWAS
results=phewas(phenotypes,genotypes,cores=1)      
#绘图
phewasManhattan(results, annotate.angle=0,title="My Example PheWAS Manhattan Plot")

3.个人解读

其实我用这个R包刚开始也很困惑,我要做的是小鼠的PheWAS(全表型组关联分析),而这个包是为人类疾病准备的(需要有ICD9/10CM),是不是这个R包不能用呢?人类的疾病有特殊的编码方式(ICD9/10CM),一个编码代表一种疾病的类型。ICD10CM编码规则较为复杂,我简单介绍下ICD9CM。以糖尿病为例,糖尿病有很多种:Type1,Type2.... ICD9CM用前三个数字表示这是糖尿病,后面两个数字表示糖尿病的类型。但我仔细看了后发现其实这些根本不用管,我可以直接自己生成表型和基因型文件,下面我会做出说明。

3.1从list提取出的三个文件介绍

id.vocab.code.count=ex$id.vocab.code.count      
#id.vocab.code.count 解读
第一列ID:代表的是个体编号(在人中就是不同的个体)。
第二列vocabulary_id:代表采用哪种编码规则(ICD9CM/ICD10CM)。
第三列code:相应的数字和字母代表一种表型。
第四列count:计数。
file
genotypes=ex$genotypes 
#基因型文件解读
第一列ID:代表个体的编号。
第二列rsExample:代表一个SNP在所有个体中的基因型。例如:aa代表0,Aa代表1,AA代表2。
file
id.sex=ex$id.sex 
#性别文件解读
#性别作为一个协变量,不要求一定要有。
第一列ID:代表个体的编号。
第二列sex:代表性别。
file

3.2 phenotypes,genotypes文件介绍

其实仔细看R包的测试就会发现前面做那么多其实就是为了产生这两个文件,这意味着可以跳过前面的步骤直接自己写脚本产生。由于genotypes文件前面已经介绍过下面就介绍下phenotypes文件。

#真正跑的时候是用phenotypes和genotypes,下面解读这两个文件
results=phewas(phenotypes,genotypes,cores=1) 
#如何看自己电脑的核心数呢?(该包支持并行计算,可有效的利用计算资源,加快计算速度)
library(parallel)        
detectCores()                          #看核心数,然后可以调整上面的cores,cores<detectCores()
#phenotypes文件介绍
第一列ID:个体的编号。
第二列到最后的列名:表示疾病的类型(ICD9/10CM,一种对应一种疾病,你如果要跑自己的数据可以自己对你的表型进行标号)
TRUE、FALSE、NA:对于人疾病来说TRUE代表患病,FALSE代表未患病,NA代表该表型缺失。
TRUE、FASLE:相当于是离散型随机变量,看过R包核心函数phewas就知道对离散型变量该包做的是逻辑斯谛回归(处理二分类)。
当然如果你的数据是连续型的(表型有具体数值),该包会对你的连续型变量做线性回归。
file

本文由博客一文多发平台 OpenWrite 发布!

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 213,417评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,921评论 3 387
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 158,850评论 0 349
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,945评论 1 285
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,069评论 6 385
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,188评论 1 291
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,239评论 3 412
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,994评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,409评论 1 304
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,735评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,898评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,578评论 4 336
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,205评论 3 317
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,916评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,156评论 1 267
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,722评论 2 363
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,781评论 2 351

推荐阅读更多精彩内容