美文网首页生物信息学
GWAS - PRS多基因风险评分计算学习笔记

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

作者: SnowPye | 来源:发表于2020-05-31 16:14 被阅读0次

    一、安装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

    相关文章

      网友评论

        本文标题:GWAS - PRS多基因风险评分计算学习笔记

        本文链接:https://www.haomeiwen.com/subject/bvilzhtx.html