GWAS - PRS多基因风险评分计算学习笔记

一、安装PRSice(mac版)

经试验我觉得直接从git hub中下载对应的安装包是最快的:https://github.com/choishingwan/PRSice,下载之后解压,解压文件如图所示

  • PRSice_mac为软件的运行脚本
  • PRSice.R是执行脚本的封装
  • TOY开头的是数据集,分为BASE和TARGET两部分
    *接下来可以先用这个数据集测试一下

跟plink的安装一样,打开terminal,cd到解压的文件夹,./PRSice_mac如果显示下面的画面表现为安装成功:



之后下载R包

Rscript PRSice.R --dir .

Rscript 表示用R脚本来调用该软件
dir参数为指定R包ggplot2安装的路径
因为结果展示会调用ggplot2画图进行可视化,如果你的R中已经安装了这个包,这个参数可以不要!!!!(我的R中有,所以后文中的此条命令我均不需要运行)
**因为直接运行该命令总容易报错,因此也可以通过R studio来安装这个包:点击面板中的package,选择ggplot2点击install,如果太慢的话选择清华的镜像就会快一点

二、文件准备

上面安装成功的图中可以看到PRSice的注释信息,需要准备base和target两个文件

  • Base Dataset:以空格分割的SNP关联分析文件,为.assoc结尾的。关联分析的教程见前面。
    必须要有的几列是:SNP, A1, OR or BETA(对于分类型表型提供OR值), P,其他几列不是必须要的,如下图:


    BASE

    如果表头没有以上述名称命名,必须用参数--chr, --A1, --A2, --stat, --snp, --bp, --se, --pvalue指明表头。
    参数可以这样设置:--snp SNP --chr CHR --bp BP --A1 A1 --A2 A2 --stat OR --se SE --pvalue P
    或者:--snp 0 --chr 1 --bp 2 --A1 3 --A2 4 --stat 5 --se 6 --pvalue 7 --index

  • Target Dataset 是plink产生的二进制文件.bed, .bim, 和 .fam file,要经过质控QC,过程翻看之前的教程。

三、PRSice运行

#二元性状用的是OR值,命令如下
Rscript PRSice.R --dir . \ 
  --prsice ./PRSice \
  --base TOY_BASE_GWAS.assoc \
  --target TOY_TARGET_DATA \
  --thread 1 \
  --stat OR \
  --binary-target T

#数量性状用的是BETA值,命令如下
Rscript PRSice.R --dir . \
  --prsice ./PRSice \
  --base TOY_BASE_GWAS.assoc \
  --target TOY_TARGET_DATA \
  --thread 1 \
  --stat BETA \
  --beta \
  --binary-target F
  • --dir参数指定可执行文件的路径;(安装过ggplot2的可以直接忽略掉这一条命令)
  • base参数指定base data(验证样本)的关联分析结果,包含每个遗传变异的效应值和 p 值。
  • target参数指定target data(测试样本)的分型结果,指 target 样本二进制 plink 格式的基因型数据前缀,也支持 BGEN 格式
  • thread参数指定线程数
  • stat指定base data关联分析结果中的效应值,有OR和BETA两种取值
  • binary-target参数指定target data中分型结果是否为二分类的表型,T=ture
    ·

  • 如果下载下来的数据第六行并不是按我们的需求写上是表型的,为了后面的关联分析,要手工定义第六行是实验组,即需要把.fam 文件中的第6行,全部变成2。操作:
#第一步: mv操作,重命名.fam文件为tmp文件,这样可以在上tmp上面修改
mv xxxx.fam tmp
#第二步:awk操作,把第六行全部变成2,然后这个东西再重新写入.fam文件
awk '{print $1,$2,$3,$4,$5,$6=2}' tmp > xxxx.fam
  • 在上述情况下,或者如果fam中没有表型文件,需要另外再创建一个表型文件
    表型文件一般表型文件为txt格式,表型文件有三列,分别为:
    FID(Family ID)、IID(Individual ID)、Pheno(Phenotype)
    在命令语句中加入【--pheno xxx.txt
  • 用bar-level设定你想计算的阈值,或者PRSice会默认计算upper和lower间一定间隔的p值下的所有PRS,命令语句中加入【--all-score】命令可以得到这所有的PRS结果

四、运行结果解释

运行完毕后,PRSice文件夹中应该有如下的文件,两个图和几个文本文件


  • PRSice_BARPLOT_.png 显示了不同阈值得到的多基因评分的对应 R2 分布。柱状图最高的点表示模型最优,如该图显示的是P值阈值为0.4463时,模型最优,R方=0.0520082(P=4.7*10-18)这些信息也可以在.summary文件中查到 。


    BARPLOT
  • PRSice_HIGH-RES_PLOT_ .png条形图显示了许多不同的 p值阈值,以及对应 R2 的 p 值(黑点),绿线为趋势线。y 轴最高点代表最优解best-fit PRS
    这里最佳的P值阈值为绿色最高点处,此时P值的阈值为0.4463


    high-res

    这两个图均表明,许多影响 base 样本性状的 SNP 可用于预测 target 样品中的性状。 注意,这两个性状可以相同也可以不同。 如果使用相同的性状,则预测值与该性状的遗传性(以及 base 样本的样本量)有关。 如果分析不同的性状,则预测值与两个性状之间的遗传重叠有关。无论哪种方式,多基因风险评分分析通常都表明,宽松 p 值阈值的模型通常比更严格阈值的模型有更好的预测能力,这表明许多统计学上无关紧要的 SNP 在多基因性状上具有预测价值。


  • PRSice.prsice文件包含:在给定不同阈值的P值以后,符合要求的SNP数量(Num_SNP),SNP的解释度(R2),回归系数
  • PRSice.best文件包含FID,IID,是否回归,PRS值。这个文件计算的是每个个体最优的PRS值。
  • PRSice.summary文件,包含表型,P的阈值,PRS的解释方差,所有变量的解释方差,协变量的解释方差,回归系数,P值,该阈值下的SNP数量。 这个文件给出的是该表型下最优的模型。
  • Clumping :识别并选择每个LD block中最显著的 SNP(即 p 值最低)以进行进一步分析。这样可以减少 SNP 之间的相关性,同时保留具有最强统计证据的 SNP。

附1:OR值与log OR值的转换

下载PGC数据库关于精神分裂症的数据https://www.med.unc.edu/pgc/data-index/
下载完成为.gz文件,解压语句:

tar -xzvf file.tar.gz      #解压xxxtar.gz文件
gunzip FileName.gz    #解压xxx.gz文件

我下载的为第二种文件,解压完成后得到一个文稿,可以用txt打开
打开之后检查一下文稿里面是OR值 或者 log OR值(计算PRS需要用OR),可以用R进行转换:

dat <- read table("xxxx.txt",header=T)
dat$OR <- exp(dat$LogOR)
write.table(dat,"xxxtransformed.txt",quote=F, row.names=F)
q()
  • 文本文件可以打开txt查看

附2:万能查看文件语句

用txt可以打开任何一个文件查看,什么后缀的都可以但是!!!!!!排列不整齐!!!!进行后续筛选非常不方便!!!!!!
这个时候sort -o可以帮你

sort -o happy.txt sad.assoc 

上面的命令就可以把关联分析得到的sad.assoc文件 另存为 happy.txt,然后就可以用EXCEL打开了,尽情查看吧。

参考资料

http://www.360doc.com/content/19/1224/13/68068867_881784568.shtml
https://choishingwan.github.io/PRSice/step_by_step/
https://www.jianshu.com/p/636048889b2a
https://www.jianshu.com/p/656573127d11
https://www.cnblogs.com/chenwenyan/p/10686136.html
https://www.jianshu.com/c/7df7c15887bd

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