#基因型质控
plink --chr-set 33 --bfile genotype864_QC --maf 0.05 --geno 0.1 --mind 0.1 --hwe 0.000001 --make-bed --out genotype864_QC1

#主成分分析
plink --chr-set 33 --bfile genotype864_QC1 --pca 10 --out genotype864_QC1

#主成分作图
setwd("E:/Duck/F2/")
pca <- read.table("genotype864_QC1.eigenvec")
eigenval <- read.table("genotype864_QC1.eigenval")
colnames(pca) <- c("FID","IID",paste0("PC", 1:(ncol(pca)-2)))
pve <- data.frame(PC = 1:10, pve = eigenval$V1/sum(eigenval$V1)*100)
pdf("PCA_explain.pdf",height = 7,width = 7)
ggplot(pve, aes(PC, pve)) + geom_bar(stat = "identity")+ ylab("Percentage variance explained") + theme_classic()
dev.off()
pdf("population_structure_PCA.pdf",height = 7,width = 7)
ggplot(pca, aes(PC1, PC2)) + geom_point(size = 2)+coord_equal()+theme_classic()+  xlab(paste0("PC1 (", signif(pve$pve[1], 3), "%)"))+ ylab(paste0("PC2 (", signif(pve$pve[2], 3), "%)"))
dev.off()

####一般线性模型GWAS
 plink --chr-set 33 --allow-no-sex --bfile genotype864_QC1 --linear --pheno phenotype_hiblup.txt --mpheno 2 --covar Covariate864.txt --covar-number 1-4 --out trait1

####GWAS作图
library(qqman)
library(data.table)
setwd("E:/Duck/F2")
gwasResults <- fread("trait1.assoc.linear", head=TRUE)
head(gwasResults)

tiff(filename ="trait1.qqplot.tiff",width = 8, height = 6, units = "in",res=300,compression = "jpeg",bg = "white")
qq(gwasResults$P, main = "Q-Q plot of GWAS p-values")
dev.off()


tiff(filename = paste(Trait,".manhattan.tiff",sep = ""),width = 8, height = 6, units = "in",res=300,compression = "jpeg",bg = "white")
manhattan(gwasResults1, chr="CHR", bp="BP", snp="SNP", p="P",cex=0.5,suggestiveline =F,genomewideline =-log10(0.00000005),col = c("blue","lightblue"))
dev.off()



####混合线性模型GWAS
#构建G矩阵
gcta64 --autosome-num 33 --bfile genotype864_QC1 --make-grm --make-grm-alg 0 --out genotype864_QC1 --thread-num 10
#GWAS运行
gcta64 --mlma --bfile genotype864_QC1 --grm genotype864_QC1 --pheno phenotype_hiblup.txt --mpheno 2 --covar dcovariate.txt --qcovar qcovariate.txt --out Trait1_gcta --thread-num 4


####GWAS作图
library(qqman)
library(data.table)
setwd("E:/Duck/F2")
gwasResults <- fread("Trait1_gcta.mlma", head=TRUE)
head(gwasResults)

tiff(filename ="trait1_gcta.qqplot.tiff",width = 8, height = 6, units = "in",res=300,compression = "jpeg",bg = "white")
qq(gwasResults$p, main = "Q-Q plot of GWAS p-values")
dev.off()

tiff(filename ="trait1_gcta.manhattan.tiff",width = 8, height = 6, units = "in",res=300,compression = "jpeg",bg = "white")
manhattan(gwasResults, chr="Chr", bp="bp", snp="SNP", p="p",cex=0.5,suggestiveline =F,genomewideline =-log10(0.00000005),col = c("blue","lightblue"))
dev.off()





#####计算遗传力
hiblup --single-trait --pheno phenotype_hiblup.txt --pheno-pos 4 --pedigree pedigree_all2.txt --dcovar 3 --out trait1





#####GBLUP
hiblup --single-trait --bfile genotype864_QC --pheno phenotype_Ref1.txt --pheno-pos 4 --dcovar 3 --add --out trait1.Ref1 


##可靠性估计
val=read.table("phenotype_Val1.txt",header = T)
predict=read.table("trait1.Ref1.rand",header=T)
mergefile=merge(val,predict, by.x = "ID", by.y = "ID")
head(mergefile)
Result=cor(mergefile["che"],mergefile["GA"],use = "complete.obs")
cat(Result,"\n")


#BLUP育种值估计
hiblup--single-trait --pedigree pedigree_all2.txt --pheno phenotype_Ref1.txt --pheno-pos 4 --dcovar 3 --out trait1.Ref1.allBLUP 

val=read.table("phenotype_Val1.txt",header = T)
predict=read.table("trait1.Ref1.allBLUP.rand",header=T)
mergefile=merge(val,predict, by.x = "ID", by.y = "ID")
head(mergefile)
Result=cor(mergefile["che"],mergefile["PA"],use = "complete.obs")
cat(Result,"\n")


##ssGBLUP育种值估计
hiblup --single-trait --pedigree pedigree_all2.txt --bfile genotype864_QC --pheno phenotype_Ref1.txt --pheno-pos 4 --dcovar 3 --out trait1.Ref1.ssGBLUP
val=read.table("phenotype_Val1.txt",header = T)
predict=read.table("trait1.Ref1.ssGBLUP.rand",header=T)
mergefile=merge(val,predict, by.x = "ID", by.y = "ID")
head(mergefile)
Result=cor(mergefile["che"],mergefile["HA"],use = "complete.obs")
cat(Result,"\n")