目录

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 valuelibrary(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)