美文网首页群体遗传学genome生信研究生必备
全基因组关联分析(GWAS)protocol

全基因组关联分析(GWAS)protocol

作者: 只看不写_nathan | 来源:发表于2022-09-07 10:29 被阅读0次

    全基因组关联分析(Genome-wide association study,GWAS)是应用基因组中数以百万级别的单核苷酸多样性(Single Nucleotide Ploymorphism,SNP)作为分子遗传标记,进行全基因组水平上表型与marker的相关性分析,从而发现影响复杂性状的基因突变的一种策略。


    使用的软件为:

    • plink
    • EMMAX
    • R
      平台为Linux

    00软件安装

    Plink作为一款老牌软件,它的文件格式是多个计算软件所需要的。Plink自身也可以进行全基因组的关联分析,但是速度较慢不推荐。我们本次使用Plink的目的是利用它对SNP数据进行提取和过滤。Plink的安装也十分简单,可以直接用conda安装:
    conda install -y plink
    EMMAX是2010年Kang开发的EMMA升级版,它的计算速度特别快,并且在棉花,大豆,水稻甚至椰子等复杂性状探究中广泛应用。每一个SNP对复杂性状的解释率都很低,每个组件的方差在运算中都只计算一次,从而达到简化计算的目的。
    软件下载地址:
    EMMAX的下载地址非常友好,不需要翻墙。

    EMMAX下载
    我用的是老版本,应该没什么区别,看介绍只有编译上的差别,然后我们安装
    tar xvf emmax-interl-binary-20120210.tar.gz
    #然后加到环境变量
    emmax-interl64
    

    R语言大家应该熟知,不需要过度介绍

    01phenotype和genotype文件准备

    phenotype文件

    emmax识别的phenotype文件需要一个表型一个单独的文本,我给大家截图展示一下: phenotype phenotype格式

    文件的格式也非常简单,只有三列。前两列为品种编号和genotype中对应,第三列为表型。

    genotype文件

    plink格式的SNP文件需要用plink进行抽取和转置:
    plink --noweb --bfile /public/home/hguo/public/data/SNP/GWAS_Analysis/Rice_GWAS/gwas_all_chr.filter --keep keep_individual.txt --make-bed --maf 0.05 --geno 0.1 --out gwas_all_chr.filter
    --bfile 二进制plink文件
    --keep 需要提取的品种集
    --make-bed 生成二进制的plink文件
    --maf 过滤最大等基因频率
    --geno 过滤缺失率
    --out 输出的文件名

    02关联分析

    亲缘关系矩阵

    emmax有子命令可以计算亲缘关系矩阵,也可以用plink去计算。
    emmax-kin-intel64 emmax_in -v -h -s -d 10 new_all_chr.filter
    最后的new_all_chr.filter是genotype名

    群体结构

    群体结构对于群体分化不明显的数据,我觉得应该是可以用加的。群体结构我们使用Admixture软件进行计算,这个很多网文都有就不过多介绍,用conda是可以直接安装的。计算之后得到最小的CV值,然后把群体结构整理成协变量的格式。这里我使用的数据是亲缘关系较近的群体,就没有加群体结构。

    计算关联性

    计算关联性之前需要将我们得到的plink文件转置,emmax接受的是转置后的tped文件,代码如下:

    plink --bfile gwas_all_chr.filter --recode12 --output-missing-genotype 0 --transpose --out new_all_chr.filter
    

    接着是计算关联分析的p值。这一步很简单,直接上代码:

    for line in $(cat name.txt)
    do
     emmax -v -d 10 -t new_all_chr.filter -p pheno/$line -k new_all_chr.filter.hIBS.kinf  -o result/$line
    done
    

    解释一下,name.txt是我的表型文件的文件名,为了循环跑每一个sample
    -t 是tped的文件名
    -p 是phenotype文件名
    -k 是刚才得到的kinship文件,如果要用群体结构的话可以加一个协方差的参数
    -o 是输出文件的文件名
    -v 我记得是输出当前进度的设置
    -d 我记得是给予的线程数(ps:这两个参数记不太清了,反正没影响

    03结果文件解析

    跑完上面的程序后,我们会得到三个文件分别是log,reml和ps文件。其中ps文件使我们需要的结果文件。目前手上的项目都已经结束,所以搬运一张ps文件的格式: ps文件格式

    这个文件四列分别为SNP_ID,回归系数,标准差和p值。
    根据这个文件,我们把数据整理成SNP,CHR,BP(snp所在染色体的物理距离)和P(关联分析的P值)四列,然后用于下一步画图。

    04结果画图

    比较简单的曼哈顿图就是使用R包qqman,只要把文件格式整理成我上述描述的那样可以直接使用下面的代码:

    temp<-list.files(pattern = "*.txt")
    for (i in 1:length(temp)) {
    library("qqman")
    mydata<-read.table(temp[i],header=TRUE,sep="\t")
    mydata$SNP<-as.character(mydata$SNP)
    mydata$CHR<-as.numeric(mydata$CHR)
    mydata1<-na.omit(mydata)
    
    outfile<-paste("MM/",temp[i],".png",sep = "")
    png(outfile,height = 900,width = 1200,pointsize = 12)
    manhattan(mydata1,main="Manhattan Plot",cex=1.5,cex.axis=1,col=c("red","blue"),suggestiveline =F,
    genomewideline=F,,chrlabs = c("1:12"))
    abline(h=-log10(1.33E-6), col="orange",lty=2,lwd=2)
    abline(h=-log10(6.66E-8), col="purple",lty=2,lwd=2)
    dev.off()
    }
    

    解释一下,*txt是当前文件夹下所有的txt文件,这个地方我已经把ps文件进行抽提整理成上述所说的那四列的形式,在这里是为了批量画图。

    最后我们看一下画出的曼哈顿图的情况: 最终结果
    这里解释一下,为什么0-1的位置没有东西。这些点都是snp,0-1之间的snp的p值为0.1-1,也就是完全不显著的snp,为了节省画图时间和图片空间,我在ps到txt的文件整理过程就给过滤了。第二点,为什么有两个阈值线,这里分别代表显著和极显著。

    好了,这就是GWAS工作的全部内容,有补充的小伙伴请在评论区留言。谢谢~

    相关文章

      网友评论

        本文标题:全基因组关联分析(GWAS)protocol

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