GWAS的曼哈顿图和QQ图diamante

科技   科技   2024-10-22 18:41   河南  

大家好,我是邓飞,对于GWAS分析结果,第一个要看的是曼哈顿图,看看有没有显著性的点,没有显著性的点,项目白做了!第二个要看的是QQ图,比较翘就非常理想。下面介绍一下常用的可视化方法,包括:qqman和cmplot两个包。

相关软件,比如gapitrMVP,都会自动出图,而GEMMAGCTA则是需要后期自己作图。

无论是软件自动出图,还是需要自己作图,学习根据GWAS结果手动作图都是必须的。

我们一般使用qqman作图和cmplot两个包画GWAS的QQ图和曼哈顿图,后者颜色更漂亮。

这篇博客,介绍一下这两个包如何画GWAS的结果可视化图。

第一个是qqman, 因为这个软件函数很方便。

「安装qqman软件」

# Install the stable release from CRAn
install.packages("qqman")
# 直接从Github上安装
devtools::install_github("stephenturner/qqman")
# 从Github上安装最新的开发版
devtools::install_github("stephenturner/qqman",ref="dev")

如果你只知道qqman,而不知道cmplot的话,这说明你还不太专业。

「安装cmplot代码」

install.packages("CMplot")

看一下cmplot包的出图:

另外,cmplot也可以出这样的图:

示例数据

使用qqman中的gwasResults数据集:

library(qqman)
data("gwasResults")

dat = gwasResults
head(dat)

str(dat)
table(dat$CHR)

第一列是SNP的ID,第二列是染色体,第三列是物理位置,第四列是P值。

共有16470个SNP数据,共有22条染色体。

qqman作图

「QQ图绘制」

这里,只需要一列P值即可。

qq(dat$P)

「曼哈顿图」如果数据结构如上所示,直接调用数据即可:

manhattan(dat)

当然,更通用的是指定染色体、物理位置、P值:

manhattan(dat,
            chr = "CHR"# 染色体
            bp = "BP"# 物理位置
            p = "P"# P值
            )

这里,还可以设置阈值线,默认是5和8:

suggestiveline Where to draw a "suggestive" line. Default -log10(1e-5). Set to FALSE to disable.

genomewideline Where to draw a "genome-wide sigificant" line. Default -log10(5e-8). Set to FALSE to disable.

cmplot作图

我们用同样的数据,使用cmplot作图。

「qq图绘制」

CMplot(dat,plot.type = "q",threshold = 0.05)

对比一下cmplotqqman的QQ图:可以看到,cmplot的QQ图更好看,而且还有置信区间!666

「曼哈顿图:」代码:

CMplot(dat,plot.type = "m",threshold = c(0.01,0.05)/nrow(dat),threshold.col=c('grey','black'),
       threshold.lty = c(1,2),threshold.lwd = c(1,1), amplify = T,
       signal.cex = c(1,1), signal.pch = c(20,20),signal.col = c("red","orange"))

比较一下两个软件的曼哈顿图:可以看到,cmplot的作图更好看,不同染色体不同的颜色,不同阈值的结果不同的颜色。

高级作图

「SNP密度图」

CMplot(dat,plot.type = "d",bin.size = 1e6, col = c("darkgreen","yellow","red"))

「环形曼哈顿图:」

CMplot(dat,plot.type="c",r=0.5,threshold=c(0.01,0.05)/nrow(pig60K),cex = 0.5, 
       threshold.col = c("red","orange"), threshold.lty = c(1,2),amplify = T, cir.chr.h = 2,
       signal.cex = c(2,2), signal.pch = c(19,20), signal.col=c("red","green"),outward=TRUE)

注意,这个图,很多参数都是默认的,具体还可以再设置美观一下,但是轮廓图是有了!

合并密度图和圆形曼哈顿图:

CMplot(dat,plot.type="c",r=0.4,col=c("grey30","grey60"),chr.labels=paste("Chr",c(1:22),sep=""),
       threshold=c(1e-6,1e-4),cir.chr.h=1.5,amplify=TRUE,threshold.lty=c(1,2),threshold.col=c("red","blue"),signal.line=1,signal.col=c("red","green"),chr.den.col=c("darkgreen","yellow","red"),
       bin.size=1e6,outward=FALSE,file="jpg",memo="",dpi=300,file.output=TRUE,verbose=TRUE)

分割线



大家好,我是邓飞,一个持续分享的农业数据分析师,这里我将自己公众号的干货内容挑重点罗列一下,方便大家阅读和使用。


1,GWAS学习教程(快来领取 | 飞哥的GWAS分析教程),这个pdf是我将公众号的内容进行了汇总,更方便从头学习GWAS分析,里面配套了数据、代码和讲解,属于干货推荐的Number 1。


2,农学人如何入门数据分析资料汇总(飞哥汇总 | 入门数据分析资源推荐),里面推荐了免费的教程,包括编程、统计和专业书籍。


3,数量遗传学电子书下载(数量遗传学,分享几本书的电子版


4,R语言电子书线上书籍推荐(学习R语言这几本电子书就够了!


育种数据分析之放飞自我
本公众号主要介绍动植物育种数据分析中的相关问题, 算法及程序代码.
 最新文章