GWAS：使用R，比较GLM和MLM对假阳性的控制差异（复刻Nature genetics 图）

目录

1.数据准备

2.代码

如果想知道横纵坐标设置的原理，移步这篇超级棒的文章！

我们来复刻如下这张2016年发表在Nature genetics上的一篇文章中比较GLM和MLM的QQ plot！

参考文献：

Genetic variation in ZmVPP1 contributes to drought tolerance in maize seedlings

Wang, X., Wang, H., Liu, S. et al. Genetic variation in ZmVPP1 contributes to drought tolerance in maize seedlings. Nat Genet 48, 1233–1241 (2016). https://doi.org/10.1038/ng.3636

Quantile–quantile plot for the GWAS under a general linear model (GLM) and mixed linear model (MLM).

 数据是随便找的，成果图基本上就是这种样式。

1.数据准备

GLM和MLM分析产生的p-value

如果是使用Tassel分析的，可以直接把结果导入下面的代码。

如果不是，直接读入GLM和MLM分析产生的p-value作为下面的GLM和MLM变量的值

2.代码

#1.Tassel分析产生的GLM和MLM数据读入====
GLM <- read.table("GLM_results.txt")
colnames(GLM) <- GLM[1,]
GLM <- GLM[-1,]
GLM$p <- as.numeric(GLM$p) #读入的数据类型均为character，需要转化为numeric data才能进行比较

MLM <- read.table("MLM_results.txt")
colnames(MLM) <- MLM[1,]
MLM <- MLM[-1,]
MLM$p <- as.numeric(MLM$p)

#2.选做：提取性状子集(如果分析包含多个性状再进行这一步，如果没有就忽略这里)====
traits_factor <- as.factor(GLM$Trait)
all_trait <- levels(traits_factor)

trait <- all_trait[i]   #i是要研究的性状在all_trait中的序号
data_glm <- subset(GLM,GLM$Trait==trait)
data_mlm <- subset(MLM,MLM$Trait==trait)
#注意看两个子集的marker数目是否相同
data_mlm <- data_mlm[-1,]

#3.计算均匀分布的分位数的负对数（生成横坐标数据）====
expected_p_value <- c(1:nrow(data_glm))/nrow(data_glm) #计算均匀分布的分位数
minus_log10_exp <- -log10(sort(expected_p_value,decreasing = TRUE))

#4.计算分析结果p的负对数（生成纵坐标数据）====
glm_minus_logp_val <- -log10(sort(data_glm$p,decreasing = TRUE)) #如果没有进行第二步，data_glm就是上面读入的GLM，data_mlm是上面的MLM
#生成绘图表格，group添加分组信息
glm_qq_data <- data.frame(minus_log10_exp=minus_log10_exp,minus_log_p_val=glm_minus_logp_val,group=rep("GLM",nrow(data_glm)))
#同理，处理MLM
mlm_minus_logp_val <- -log10(sort(data_mlm$p,decreasing = TRUE))
mlm_qq_data <- data.frame(minus_log10_exp=minus_log10_exp,minus_log_p_val=mlm_minus_logp_val,group=rep("MLM",nrow(data_mlm)))

#5.最终绘图数据及绘图====
combine_data <- rbind(glm_qq_data,mlm_qq_data) #合并GLM和MLM的p value

library(ggplot2)
p <- ggplot()+
  geom_point(data = combine_data,mapping = aes(minus_log10_exp,minus_log_p_val,color=group))+ 

  geom_abline(intercept = 0,slope = 1,color="#D3D3D3",size=0.75)+  #画y=x直线
  
  scale_color_manual(values = c("black","red"))+  #黑色点为GLM结果，红色点为MLM结果
  
  scale_x_continuous(expand = c(0,0))+
  scale_y_continuous(expand = c(0,0))+
  
  labs(title = trait)+ #添加主标题
  xlab(expression('Expected'~'-log'[10]*'(p)'))+ylab(expression('Observed'~'-log'[10]*'(p)'))+

  theme_classic()+                      #确定基础主题
  
  theme(
    axis.line = element_line(size = 1), #修改坐标轴粗细
    axis.ticks = element_line(size=1),  #修改坐标轴刻度线粗细
    axis.ticks.length = unit(0.25,"cm"),#修改坐标轴刻度线长度
    axis.text = element_text(size = 12,colour = "black"), #修改坐标轴刻度字体
    axis.title = element_text(size = 12), #修改坐标轴标签字体
    plot.title = element_text(size = 12,hjust = 0.5),#修改图题位置和字体大小
    legend.title = element_blank(),      #除去图例标题
    legend.text = element_text(size=12) 
    )

成果图！（主标题打码了。。。）

 如果想改变图例的位置的话，可以直接在AI里面改，或者再在theme里面设置。

如果想知道横纵坐标设置的原理，移步这篇超级棒的文章！

如何理解GWAS中Manhattan plot和QQ plot所传递的信息 - 简书 (jianshu.com)