多个性状曼哈顿图

library(qqman)# 创建示例数据 - 改为3个染色体set.seed(123)n_snps<-900n_chr<-3# 创建数据data<- data.frame(SNP=paste0("rs",1:n_snps), CHR=rep(1:n_chr, each=n_snps/n_chr), BP=unlist(lapply(1:n_chr, function(x)sort(sample(1:2e6, n_snps/n_chr)))), P1=runif(n_snps), P2=runif(n_snps), P3=runif(n_snps))# 计算-log10(P)data$logP1<- -log10(data$P1)data$logP2<- -log10(data$P2)data$logP3<- -log10(data$P3)# 第一步：计算染色体边界，增加间距data$x_pos<- NA lastbase<-0ticks<- NULL chr_ends<- c()# 增加染色体间距的倍数space_multiplier<-0.3# 控制间距大小，值越大间距越大for(iin1:n_chr){chr_data<- data[data$CHR==i,]chr_length<- max(chr_data$BP, na.rm=TRUE)- min(chr_data$BP, na.rm=TRUE)if(i==1){data$x_pos[data$CHR==i]<- data$BP[data$CHR==i]}else{# 增加间距：前一个染色体的长度 + 间距lastbase<- lastbase + max(data$BP[data$CHR==(i-1)], na.rm=TRUE)+(chr_length * space_multiplier)data$x_pos[data$CHR==i]<- data$BP[data$CHR==i]+ lastbase}# 记录染色体结束位置chr_ends<- c(chr_ends, max(data$x_pos[data$CHR==i], na.rm=TRUE))# 计算染色体中心位置（用于x轴标签）ticks<- c(ticks,(min(data$x_pos[data$CHR==i], na.rm=TRUE)+ max(data$x_pos[data$CHR==i], na.rm=TRUE))/2)}# 第二步：创建空绘图区域par(mar=c(5,5,4,2))pdf(file="1.pdf")plot(NA, xlim=c(min(data$x_pos, na.rm=TRUE), max(data$x_pos, na.rm=TRUE)), ylim=c(0, max(c(data$logP1, data$logP2, data$logP3), na.rm=TRUE)+1), xlab="Chromosome", ylab=expression(-log[10](italic(p))), xaxt="n",# 不显示默认x轴yaxt="n",# 不显示默认y轴bty="n", main="Three Chromosomes with Increased Spacing", cex.main=1.2)# 第三步：为每个性状单独绘制# 性状1 - 用三角形points(data$x_pos, data$logP1, col=adjustcolor("blue", alpha.f=0.7), pch=17, cex=0.8)# 性状2 - 用圆形points(data$x_pos, data$logP2, col=adjustcolor("red", alpha.f=0.7), pch=19, cex=0.8)# 性状3 - 用方形points(data$x_pos, data$logP3, col=adjustcolor("green", alpha.f=0.7), pch=15, cex=0.8)# 添加染色体分隔线（更明显）if(n_chr>1){for(iin1:(n_chr-1)){# 染色体间的分隔线abline(v=chr_ends[i]+(min(data$x_pos[data$CHR==(i+1)], na.rm=TRUE)- chr_ends[i])/2, col="gray60", lty=2, lwd=1.5)}}# 添加x轴axis(1, at=ticks, labels=paste("Chr",1:n_chr), cex.axis=1, font=2)# 添加y轴y_ticks<- pretty(c(0, max(c(data$logP1, data$logP2, data$logP3), na.rm=TRUE)))axis(2, at=y_ticks, las=1, cex.axis=0.9)# 添加网格线abline(h=y_ticks, col="gray90", lty=3, lwd=0.5)# 添加显著线abline(h=-log10(5e-8), col="red", lty=2, lwd=2)text(x=max(data$x_pos, na.rm=TRUE)*0.05, y=-log10(5e-8)+0.2, labels="5e-8", col="red", cex=0.8)# 添加图例legend("topright", legend=c("Trait 1","Trait 2","Trait 3"), col=c("blue","red","green"), pch=c(17,19,15), bty="o", box.col="gray",bg="white", cex=0.9, pt.cex=1.2)# 添加边框box(col="gray60", lwd=1)dev.off()