GWAS（3）——曼哈顿图和QQ图，PCA plot

曼哈顿图

GAPIT作图

###GAPIT 网站：http://www.zzlab.net/GAPIT/
###源码：http://www.zzlab.net/GAPIT/gapit_functions.txt
###R代码 
###加载包 
library(multtest) 
library(gplots) 
library(LDheatmap) 
library(genetics) 
library(ape) 
library(EMMREML) 
library(compiler) 
source("http://zzlab.net/GAPIT/gapit_functions.txt") 
source("http://zzlab.net/GAPIT/emma.txt") 
setwd("你的路径") 
myG <- read.table("trait.txt",head=TRUE) 
myY <- read.table("5out.hmp.txt",head=F) 
myGAPIT <- GAPIT( 
 Y = myY, 
 G = myG, 
 PCA.total = 5 ###展示 PCA 的主成分数量 
) 

这样就画出曼哈顿图和 QQ 图了。
只是画出的图没有下面的 GEMMA 方法效果好。

GEMMA分析

在分析前需要计算 PCA，并且准备表型文件。

###准备表型文件，需要将 fam 文件中的第六列（一般是-9 表示缺失）替换掉，转换成手头的表型数据。 

library(data.table) 
phe <- fread(“phe.txt”,header=T,sep=”\t”) 
colnames(phe) <- c(“V1”,colnames(phe)[-1]) 
fam <- fread(“5out.fam”) 
fam <- fam[,-6] 
fam.o <- merge(x=fam,y=phe,by=”V1”,all.x=T) 
fwrite(x=fam.o,file=”5out.fam.o”,sep=” “,col.names=F,na=-9,quote=F) 

也可以这样进行，前提是表型数据按照 fam 中基因型数据整理对应。

library(data.table) 
fam <- fread(“5out.fam”) 
fam <- fam[,-6] 
write.table(fam,file=”fam.csv”,row.names=F,col.names=F) 
#得到 fam.csv 后，打开文件手动粘贴表型文件数值到 fam.csv 中保存 
a<-read.table(“fam.csv”,header=F,sep=”\t”) 
write.table(a,file=”5out.fam”, row.names=F,col.names=F) 

###以上步骤需要检验一下是否适用于其他情况。 

接着，进行 kinship 矩阵运算，一般命令为

###5out 是二进制文件前缀，-gk 2 指生成 kinship 矩阵时进行 scale，-o 输出文件，后面可以
指定输出文件前缀。 
gemma –bfile 5out –gk 2 –o kin

进行 GWAS 分析

#-k 后面是 kinship 文件的全称 
#-c 后面是协变量文件，根据 pca 文件得来，需要将 pca 文件中后缀为.eigenvec 的文件修
改一下：第一第二列基因 ID 删除，在第一列添加 1，表示截距 
#-n 后面的数字代表了文件的第几列，注意的是，“-n 1”表示 fam 文件第六列，以此类推。
gemma –bfile 5out –k kin.sXX.txt –c cov.txt –n [num] –lmm [num] –o

多性状时可以进行循环命令 ：

KIN=/你的 kinship 文件路径/kinship.sXX.txt 
for i in {1..性状数} 
do 
echo $i 
gemma –bfile 5out –k $KIN –lmm –n $i –o T$i-GWA 
done

这样就输出了一连串的 assoc.txt 文件。
接下来画图可以用 qqman，rMVP 等进行画图。

rMVP代码

library(rMVP)
MVP.Data(fileHMP="gwas.MVP.hmp.txt",
         filePhe="gwas.MVP.pheno.txt",
         sep.hmp="\t",
         sep.phe="\t",
         SNP.effect="Add",
         fileKin=FALSE,
         filePC=FALSE,
         out="hmp"
)

genotype <- attach.big.matrix("hmp.geno.desc")
phenotype <- read.table("hmp.phe",head=TRUE)

map <- read.table("hmp.geno.map" , head = TRUE)

for(i in 3:ncol(phenotype)){
  phe <- phenotype[, c(1,5)]
  head <- colnames(phe)[2]
  phena <- phe %>% 
    filter(is.na(phe[[head]]))
  pheno <- na.omit(phe)
  pheno[,2] <- pheno[sample(1:nrow(pheno), replace=F),2]
  phe <- rbind(pheno,phena)
  imMVP <- MVP(
    phe=phe,
    geno=genotype,
    map=map,
    nPC.GLM=5,
    #    nPC.MLM=3,
    # nPC.FarmCPU=3,
    priority="speed",
    vc.method="BRENT",
    maxLoop=10,
    method.bin="static",
    threshold=0,
    method=c("GLM")
    # permutation.threshold=TRUE,
    # permutation.rep=1000
    
  )
  gc()
}

QQ图

rMVP

library(rMVP)
library(dplyr)
df <- read.csv("X10.csv",header=T)
df <- df[c("Marker","Chr","Pos","p")]
names(df) <- c("rs","chr","pos","P")
df[df$chr==11,]$pos <- df[df$chr==11,]$pos*1000
df <- na.omit(df)
# MVP.Report(df,plot.type="m",LOG10=TRUE,threshold=0.05/nrow(df),threshold.lty=2,threshold.lwd=1,threshold.col="grey",amplify=TRUE,col=c("dodgerblue4","deepskyblue"),signal.col="red",cex=0.7,chr.den.col=NULL,file="pdf",memo="",dpi=600)
# MVP.Report(df,plot.type="q",conf.int.col=NULL,box=TRUE,file="pdf",dpi=600)
MVP.Report(df,plot.type = "m",file.type = "jpg")
MVP.Report(df,plot.type="q",file.type = "jpg")

qqman

install.packages("qqman")
library(qqman)
qq(gwasResults$P, main = "Q-Q plot of GWAS p-values", xlim = c(0, 7), ylim = c(0,    12), pch = 18, col = "blue4", cex = 1.5, las = 1)

PCA plot

rm(list = ls())
library(ggplot2)
library(RColorBrewer)
data <- read.table("SNP.eigenvec",sep = "\t",header = T)
ggplot(data,aes(x=pc1,y=pc2,color=species))+ geom_point()
#######
tiff(filename = "pca.tif", width = 30, height = 16, units = "cm", compression = "lzw", res = 600)
ggplot(data,aes(x=pc1,y=pc2,color=species))+ 
  geom_point()+ 
  theme_bw() +
  theme(panel.border=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank(),axis.line= element_line(colour = "black"))
dev.off()