GWAS分析

首先准备输入文件（vcf文件和表型文件）

基因型推断

beagle -Xmx4g  -Djava.io.tmpdir=./TMP  gt=./yuanguoxin.vcf   out=yuanguoxin.impute impute=true window=10 nthreads=2
#用beagle软件做基因型推断

格式转换

plink --vcf ../data/Real_snp.vcf --maf 0.05 --geno 0.1  --recode12  --output-missing-genotype 0 --transpose --out snp_filter   --set-missing-var-ids @:#  --allow-extra-chr
##转换为tped格式，emmax软件所必须的

会生成 tfam、tped、map文件
根据tfam文件ID将你的表型文件排序

表型文件
tfam文件

die "perl $0 <tfam> <pheno.table> >pheno.emmax.table\n" unless @ARGV==2;
my $tfam = shift;
my $trait = shift;

my %trait;
open IN, $trait || die $!;
while(<IN>){
        chomp;
        next if (/^\s*$/);
        my @t = split /\s+/;
        next if ( scalar(@t) != 2);
        if($t[1] =~ /^[0-9.]+$/){
                $trait{$t[0]} = $t[1];
#               print "phe:$_\n";
        }
}
close IN;

open FM, $tfam || die $!;
while(<FM>){
        chomp;
        next if (/^\s*$/);
        my @t = split /\s+/;
        if(exists $trait{$t[0]}){
                print "$t[0]\t$t[0]\t$trait{$t[0]}\n";
        }
        else{
                print "$t[0]\t$t[0]\tNA\n";
        }
}
close FM;

生成亲缘关系矩阵

/home/data/t0202008/software/emmax-kin-intel64 -v -d 10  -o ./pop.kinship  snp_filter

GWAS

/home/data/t0202008/software/emmax-intel64 -v -d 10 -t snp_filter  -p trait.sort.txt -k pop.kinship -o emmax.out 1> emmax.log 2>emmax.err

GWAS结果

画图