if(!require("qqman"))install.packages("qqman")
if(!require("CMplot"))install.packages("CMplot")
library(qqman)
#library(CMplot)
library(tidyverse)
setwd('F:/R/tmp/')
#输入文件格式如下：
#  SNP   CHR   BP   P
# 其中SNP列为SNP名称，这一列可以没有
# CHR，BP分别为SNP所在的染色体和位置
# P表示该SNP的关联显著性p值
data<-readxl::read_xlsx("rMVP.xlsx",sheet = 4)
#use search() to determine whether all the packages have loaded
# 取bonferroni矫正阈值
fdr=0.01/3352885
data1<-data %>% select(,1:4)
#CMplot(data1,plot.type = "m",multracks = TRUE,threshold = fdr,dpi = 300,file = "jpg")
#data(package="CMplot")#查看R包带的数据集
#view("pig60k")查看R包示例数据
#计算染色体长度
chr_len <-data1 %>% group_by(CHR) %>%
  summarise(chr_len=max(BP))
#计算每条chr的初始位置
chr_pos<- chr_len %>%
  mutate(total=cumsum(chr_len)-chr_len) %>%
  select(-chr_len)
#计算累计SNP的位置
Snp_pos<- chr_pos %>%
  left_join(data1, .,by="CHR")%>%
  arrange(CHR,BP) %>%
  mutate(BPcum=BP+total)
head(Snp_pos,2)

X_axis <-  Snp_pos %>% group_by(CHR) %>% summarize(center=( max(BPcum) +min(BPcum) ) / 2 )
p <- ggplot(Snp_pos, aes(x=BPcum, y=-log10(P))) +
  #设置点的大小，透明度
  geom_point(aes(color=as.factor(SNP)), alpha=0.8, size=1.3) +
  #设置颜色
  #scale_color_manual(values = rep(c("grey", "skyblue"), 22 )) +
  #设定X轴
  scale_x_continuous( label = X_axis$CHR, breaks= X_axis$center ) +
  #去除绘图区和X轴之间的gap
  scale_y_continuous(expand = c(0, 0) )+
  labs(x="Chromosome",y="-log10(Pvalue)",title = "")+
  theme(legend.title=element_blank())+
  theme(legend.text = element_text(face = "bold",family = "serif"))+
  theme(axis.text.x = element_text(size = 10,colour = "black",face = "bold",family = "serif"))+
  theme(axis.text.y = element_text(size = 10,colour = "black",face = "bold",family = "serif"))+
  theme(axis.title.x = element_text(size = 12,face = "bold",family = "serif"))+
  theme(axis.title.y = element_text(size = 12,face = "bold",family = "serif"))
ggsave(p,filename = "different_condition.png")

#部分内容来自https://blog.csdn.net/ddxygq/article/details/103555955