CMplot & rMVP | 全基因组曼哈顿图和QQ图轻松可视化！

1.CMplot

1.1 CMplot介绍

CMplot：https://github.com/YinLiLin/CMplot

这是一个做全基因组对SNP可视化神器了，尹立林教授写的R包。主打两个功能，曼哈顿图（Manhattan plots）和QQ图（Quantile-Quantile plots）

• 曼哈顿图：用于展示 GWAS 分析的结果，其中每个点代表一个单核苷酸多态性（SNP），点的位置表示 SNP 在基因组中的位置，而点的高度则表示该 SNP 与研究性状的关联强度（通常用 -log10(p-value) 表示）。
• QQ图：用于检查数据的分布是否符合期望的分布（通常是正态分布）。在遗传学研究中，QQ图可以用来检测数据中是否存在过度的分散或聚集，这可能是由群体结构、样本大小或其他因素引起的。

1.2 CMplot-DEMO

安装R包，及加载数据

install.packages("CMplot")
library("CMplot")

# github最新版本
source("https://raw.githubusercontent.com/YinLiLin/CMplot/master/R/CMplot.r")

常见的SNP密度图：

CMplot(pig60K,plot.type="d",bin.size=1e6,chr.den.col=c("#7CC767", "#088247", "black"),file="pdf",file.name="",dpi=300,
       main="illumilla_60K",file.output=TRUE,verbose=TRUE,width=9,height=6)

环状曼哈顿图：

CMplot(pig60K,type="p",plot.type="c",chr.labels=paste("Chr",c(1:18,"X","Y"),sep=""),r=0.4,cir.axis=TRUE,
       outward=FALSE,cir.axis.col="black",cir.chr.h=1.3,chr.den.col="black",file="pdf",
       file.name="",dpi=300,file.output=TRUE,verbose=TRUE,width=10,height=10)

矩阵曼哈顿图：

CMplot(pig60K, plot.type="m", LOG10=TRUE, ylim=NULL, threshold=c(1e-6,1e-4),threshold.lty=c(1,2),
       threshold.lwd=c(1,1), threshold.col=c("black","grey"), amplify=TRUE,bin.size=1e6,
       chr.den.col=c("#7CC767", "#088247", "black"),signal.col=c("#D20A13","#223D6C"),signal.cex=c(1.5,1.5),
       signal.pch=c(19,19),file="jpg",file.name="",dpi=300,file.output=TRUE,verbose=TRUE,
       width=14,height=6)

单条染色体可视化：

CMplot(pig60K[pig60K$Chromosome==13, ], plot.type="m",LOG10=TRUE,col=c("grey60"),highlight=SNPs,
       highlight.col="darkgreen",highlight.cex=1,highlight.pch=19,file="jpg",file.name="", 
       threshold=c(1e-6,1e-4),threshold.lty=c(1,2),threshold.lwd=c(1,2), width=9,height=6,
       threshold.col=c("red","blue"),amplify=FALSE,dpi=300,file.output=TRUE,verbose=TRUE)

标出指定基因/SNP：

SNPs <- pig60K[pig60K[,5] < (0.05 / nrow(pig60K)), 1]
genes <- paste("GENE", 1:length(SNPs), sep="_")
set.seed(666666)
CMplot(pig60K[,c(1:3,5)], plot.type="m",LOG10=TRUE,col=c("grey30","grey60"),highlight=SNPs,
         highlight.col=rep(c("yellow","red"),length=length(SNPs)),highlight.cex=1, highlight.text=genes,      
         highlight.text.col=rep("darkgreen",length(SNPs)),threshold=0.05/nrow(pig60K),threshold.lty=2,   
         amplify=FALSE,file="jpg",file.name="",dpi=300,file.output=TRUE,verbose=TRUE,width=14,height=6)

多重曼哈顿图：

CMplot(pig60K, plot.type="m",multraits=TRUE,threshold=c(1e-6,1e-4),threshold.lty=c(1,2), 
        threshold.lwd=c(1,1), threshold.col=c("black","grey"), amplify=TRUE,bin.size=1e6,
        chr.den.col=c("darkgreen", "yellow", "red"), signal.col=c("red","green"),
        signal.cex=1, file="jpg",file.name="",dpi=300,file.output=TRUE,verbose=TRUE,
        points.alpha=100,legend.ncol=1, legend.pos="left")

QQ图：

CMplot(pig60K,plot.type="q",box=FALSE,file="jpg",file.name="",dpi=300,
    conf.int=TRUE,conf.int.col=NULL,threshold.col="red",threshold.lty=2,
    file.output=TRUE,verbose=TRUE,width=5,height=5)

多组QQ图：

CMplot(pig60K,plot.type="q",col=c("dodgerblue1", "olivedrab3", "darkgoldenrod1"),multraits=TRUE,
        threshold=1e-6,ylab.pos=2,signal.pch=c(19,6,4),signal.cex=1.2,signal.col="red",
        conf.int=TRUE,box=FALSE,axis.cex=1,file="jpg",file.name="",dpi=300,file.output=TRUE,
        verbose=TRUE,ylim=c(0,8),width=5,height=5)

1.3 CMplot参数

• Pmap：数据框，至少四列。第一列是SNP的名称，第二列是SNP的染色体，第三列是SNP的位置，剩余的列是每个性状的P值（每个性状占据一列）。
• col：向量或矩阵，如果是向量，则每个圈使用相同的颜色，表示相同染色体使用相同颜色；如果是矩阵，则每行是一个圈（性状），列是用户想要使用的不同颜色，缺失值可以用NA替换。
• pch：点的类型，与中的"pch"相同。
• band：染色体之间的间隔，默认为1。
• cir.band：圈之间的间隔，默认为1。
• H：每个圈的高度，每个圈代表一个性状，默认为1。
• ylim：Y轴的范围。
• cex.axis：X轴数字和标签的大小。
• plot.type：要绘制的图表类型，可以是"c"、“m”、“q"或"b”。
• multracks：是否将曼哈顿图绘制在多个轨道上。
• cex：点的大小。
• r：圆的半径。
• xlab：X轴标签。
• ylab：Y轴标签。
• xaxs：X轴刻度样式。
• yaxs：Y轴刻度样式。
• outward：是否将点从内向外绘制。
• threshold：显著性水平。
• threshold.col：阈值水平线的颜色。
• threshold.lwd：阈值水平线的宽度。
• threshold.lty：阈值水平线的类型。
• amplify：是否突出显示显著点。
• signal.cex：显著点的大小。
• signal.pch：显著点的形状。
• signal.col：显著点的颜色。
• cir.chr：是否绘制染色体边界。
• cir.chr.h：边界的宽度。
• cir.chr.col：边界的颜色。
• chr.labels：圆形曼哈顿图中染色体的标签。
• cir.legend：是否添加每个圈的图例。
• cir.legend.cex：图例中数字的大小。
• cir.legend.col：图例轴的颜色。
• LOG10：是否将P值转换为log10(P值)。
• box：是否在曼哈顿图周围绘制边框。
• conf.int：是否在Q-Q图上绘制置信区间。
• conf.int.col：置信区间的颜色。
• plot0：是否绘制染色体编号为0的染色体。
• file.output：是否输出图表结果。
• file：输出的文件格式，可以是"jpg"、“pdf"或"tiff”。
• dpi：.jpg和.tiff文件的像素密度，默认为300。

2.rMVP

2.1 rMVP介绍

Github：https://github.com/xiaolei-lab/rMVP

这也是一个做全基因组可视化的R包，同课题组开发，CMplot被集成进来了，特点是并行和算法优化。

• rMVP 支持包括常见的单变量线性/混合线性模型（如 MLM、GLM）以及更先进的GWAS方法，，比如 FarmCPU、SUPER等。

• 高效的内存和计算管理，针对大规模数据集设计，rMVP 优化了内存使用和计算速度，可以有效处理包含数百万个标记的大型基因组数据集。

• 多种基因组选择算法，提供了包括最佳线性无偏预测（BLUP）、贝叶斯方法等多种基因组选择算法，帮助研究人员评估和选择优良性状。

2.2 rMVP-DEMO

曼哈顿图和QQ图总的来说大同小异，多了个PCA和柱状图，毕竟亮点还是在处理大数据并行。

MVP.Hist(phe=phenotype, file.type="jpg", breakNum=18, dpi=300)

pca <- attach.big.matrix("mvp.pc.desc")[, 1:3]
#pca <- prcomp(t(as.matrix(genotype)))$x[, 1:3]
MVP.PCAplot(PCA=pca, Ncluster=3, class=NULL, col=c("red", "green", "yellow"), file.type="jpg")

2.3 rMVP参数

MVP 包中包含三种模型：广义线性模型（GLM）、混合线性模型（MLM）和 FarmCPU。

• phe, 表型数据
• geno, 基因型数据
• map, 坐标数据
• K, 关系矩阵
• CV.GLM, 在 GLM 中添加的协变量
• CV.MLM, 在 MLM 中添加的协变量
• CV.FarmCPU, 在 FarmCPU 中添加的协变量
• nPC.GLM, 在 GLM 中添加的主成分的前几列数目
• nPC.MLM, 在 MLM 中添加的主成分的前几列数目
• nPC.FarmCPU, 在 FarmCPU 中添加的主成分的前几列数目
• priority, 在计算基因组关系矩阵时选择“速度”或“内存”
• ncpus, 用于并行计算的 CPU 数目，如果未设置，默认使用所有 CPU
• vc.method, 方差成分分析的方法，有三种方法可用：“BRENT”、“EMMA”和“HE”
• maxLoop, 仅适用于 FarmCPU 的参数，FarmCPU 允许的最大迭代次数
• method.bin, 仅适用于 FarmCPU 的参数，有两个选项：“static”和“FaST-LMM”
• permutation.threshold, 如果为 TRUE，将在曼哈顿图中使用置换阈值。表型数据会被置换以打破与基因型的关系。该实验会重复多次。记录所有实验的最小 P 值向量，并建议使用该向量的 95% 分位数作为显著阈值
• permutation.rep, 置换重复次数，仅在 permutation.threshold 为 TRUE 时使用
• threshold, 0.05/标记大小，曼哈顿图上的截止线
• method, 关联测试的模型，MVP 中可用的三种模型是“GLM”、“MLM”和“FarmCPU”，可选择一种、两种或三种模型进行关联测试
• file.output, 布尔值或字符串向量。如果为 TRUE，输出所有类型的文件。如果为 FALSE，则不输出文件。对于字符串向量，可用值为 c(“pmap”, “pmap.signal”, “plot”, “log”)。其中，pmap 代表所有 SNP P 值文件，pmap.signal 代表显著 SNP 文件，Plot 代表可视化结果，log 代表日志文件。