Java自学者论坛

 找回密码
 立即注册

手机号码,快捷登录

恭喜Java自学者论坛(https://www.javazxz.com)已经为数万Java学习者服务超过8年了!积累会员资料超过10000G+
成为本站VIP会员,下载本站10000G+会员资源,会员资料板块,购买链接:点击进入购买VIP会员

JAVA高级面试进阶训练营视频教程

Java架构师系统进阶VIP课程

分布式高可用全栈开发微服务教程Go语言视频零基础入门到精通Java架构师3期(课件+源码)
Java开发全终端实战租房项目视频教程SpringBoot2.X入门到高级使用教程大数据培训第六期全套视频教程深度学习(CNN RNN GAN)算法原理Java亿级流量电商系统视频教程
互联网架构师视频教程年薪50万Spark2.0从入门到精通年薪50万!人工智能学习路线教程年薪50万大数据入门到精通学习路线年薪50万机器学习入门到精通教程
仿小米商城类app和小程序视频教程深度学习数据分析基础到实战最新黑马javaEE2.1就业课程从 0到JVM实战高手教程MySQL入门到精通教程
查看: 38|回复: 0

GWAS:拒绝假阳性之case和control数量比例严重失衡的解决方案(SAIGE模型的应用)

[复制链接]
  • TA的每日心情
    奋斗
    昨天 21:54
  • 签到天数: 95 天

    [LV.6]常住居民II

    1378

    主题

    1435

    帖子

    5万

    积分

    管理员

    Rank: 9Rank: 9Rank: 9

    积分
    50644
    发表于 2021-6-13 15:47:55 | 显示全部楼层 |阅读模式

    一、为什么要校正case和control数量比例不平衡情况

    试问作为生信届人员,最怕的是什么,当然是统计结果不靠谱。统计结果不靠谱包括两方面:一个是假阴性,一个是假阳性。假阴性可以理解为白天鹅被误当成丑小鸭了,假阳性可以理解为一大堆青蛙,你不知道哪个才是你的真命天子。假阴性就罢了,最多让你错过发现真理的机会,但万一假阳性呢,你拿着一个看似完美的结果吭哧吭哧做实验验证,一年半载的周期下来,什么结果都验证不出来,岂不是坑了做实验的人。因此,我们就要在源头上,把这个不靠谱的统计结果杜绝出去。

    上一篇文章什么!GWAS研究中case和control的比例是有讲究的?就讲到GWAS分析中,如果case和control数量比例失衡的话,会产出非常多的假阳性结果,而用SAIGE模型做GWAS分析可以校正这种数量比例不平衡的情况。下面具体讲讲怎么应用SAIGE模型。

     

    二、怎么校正:SAIGE的下载和安装

    1、下载SAIGE

    此操作在Linux上进行,系统要求R-3.5.1, gcc >= 5.5.0, cmake 3.8.1

    wget https://github.com/weizhouUMICH/SAIGE/blob/master/SAIGE_0.35.8.1_R_x86_64-pc-linux-gnu.tar.gz
    

      

     

    2、安装SAIGE

    R CMD INSTALL SAIGE_XX_R_x86_64-pc-linux-gnu.tar.gz
    

      

     

    3、安装SAIGE所依赖的其他R包:Rcpp, RcppArmadillo, RcppParallel, data.table, SPAtest, RcppEigen, Matrix, methods, optparse

    以下两个方法二选一:

    如果是用conda的话,则用以下命令:

    conda install -n r-env r-Rcpp #r-env是指conda下的R环境
    conda install -n r-env r-RcppArmadillo
    conda install -n r-env r-RcppParallel
    conda install -n r-env r-SPAtest
    conda install -n r-env r-RcppEigen
    conda install -n r-env r-optparse
    

      

    也可以进入R,用install.packages( )安装:

    install.packages("Rcpp")
    install.packages("RcppArmadillo")
    install.packages("RcppParallel")
    install.packages("SPAtest")
    install.packages("RcppEigen")
    install.packages("optparse")
    

      

     

    三、怎么校正:SAIGE的分析、解读

    1、第一步,计算NULL GLMM

    1)计算NULL GLMM的命令:

    Rscript step1_fitNULLGLMM.R \  
    --plinkFile=./input/plinkforGRM_1000samples_10kMarkers \  
    --phenoFile=./input/pheno_1000samples.txt \  
    --phenoCol=y \   
    --covarColList=x1,x2 \  
    --sampleIDColinphenoFile=IID \  
    --traitType=binary \  
    --outputPrefix=./output/example \  
    --nThreads=4 \  
    --LOCO=TRUE
    

      

    plinkFile为plink的输入文件(bed, bim, fam格式)

    phenoFile文件格式如下,第一列代表研究的表型,第二列及第N-1列代表协变量,最后一列IID为个体的ID号:

     

    --phenoCol=y # 指定你要研究的表型列名,在本次例子中,指定y的表型分析。

    --covarColList=x1,x2 #指定加入的协变量

    --sampleIDColinphenoFile=IID #指定样本的ID

    --traitType=binary  #指定研究的表型的类型,binary指二分类,即case和control

    --outputPrefix #生成文件的输出路径

     

    2)输出文件的结果解读:

    这个步骤会生成三个文件,分别为:example.rda、example.varianceRatio.txt、example_30markers.SAIGE.results.txt

    第一个文件:example.rda,是一个model file

    可以用R打开:

    load("./output/example.rda")
    names(modglmm)
    modglmm$theta
    

      

    第二个文件:example_30markers.SAIGE.results.txt,随意选取位点的关联分析结果

     

    第三个文件:example.varianceRatio.txt

     

    2、计算每个marker的SPA得分

    1)计算每个marker的SNP得分命令:

    Rscript step2_SPAtests.R \  
    --dosageFile=./input/dosage_10markers.txt \  
    --dosageFileNrowSkip=1 \  
    --dosageFileNcolSkip=6 \  
    --dosageFilecolnamesSkip=CHR,SNP,CM,POS,EFFECT_ALLELE,ALT_ALLELE \  
    --minMAF=0.0001 \  
    --sampleFile=./input/sampleIDindosage.txt \  
    --GMMATmodelFile=./output/example.rda \  
    --varianceRatioFile=./output/example.varianceRatio.txt \  
    --SAIGEOutputFile=./output/example.plainDosage.SAIGE.txt \  
    --numLinesOutput=2 \        
    --IsOutputAFinCaseCtrl=TRUE
    

      

    dosage_10markers.txt的文件格式如下,类似于plink的ped格式,前六列分别为:CHR, SNP, CM, POS, COUNTED, ALT, 后面为个体的ID号:

     

    sampleIDindosage.txt文件为样本ID名:

     

    example.rda、example.varianceRatio.txt为第一步生成的两个文件。

     

    2)输出文件的结果解读:

    生成example.plainDosage.SAIGE.txt文件,其内容如下:

     

    其中,P值(红框)即为我们校正case和control数量比例不平衡以后得到的GWAS结果,p.value.NA为不校正case和control数量不平衡的结果。

    参数说明:

    CHR: chromosome

    POS: genome position 

    SNPID: variant ID

    Allele1: Ref allele

    Allele2: Alt allele

    AC_Allele2: allele count of Alt allele

    AF_Allele2: allele frequency of Alt allele

    N: sample size

    BETA: effect size

    SE: standard error of BETA

    Tstat: score statistic

    p.value: p value with SPA applied

    p.value.NA: p value when SPA is not applied

    Is.SPA.converge: whether SPA is converged or not

    varT: estimated variance of score statistic with sample related incorporated

    varTstar: variance of score statistic without sample related incorporated

    AF.Cases: allele frequency of allele 2 in cases (only for binary traits and if --IsOutputAFinCaseCtrl=TRUE)

    AF.Controls: allele frequency of allele 2 in controls (only for binary traits and if --IsOutputAFinCaseCtrl=TRUE)

     

     至此,校正GWAS分析中case和control数量比例严重失衡的工作就完成了。导师再也不用担心你给出一堆假阳性结果了。

     

    哎...今天够累的,签到来了1...
    回复

    使用道具 举报

    您需要登录后才可以回帖 登录 | 立即注册

    本版积分规则

    QQ|手机版|小黑屋|Java自学者论坛 ( 声明:本站文章及资料整理自互联网,用于Java自学者交流学习使用,对资料版权不负任何法律责任,若有侵权请及时联系客服屏蔽删除 )

    GMT+8, 2021-8-3 09:42 , Processed in 0.105883 second(s), 29 queries .

    Powered by Discuz! X3.4

    Copyright © 2001-2020, Tencent Cloud.

    快速回复 返回顶部 返回列表