【R画图学习13.3】散点图---曼哈顿图

曼哈顿图，我自己最常用于的场景就是GWAS分析之后，常使用它展示基因组中与某种表型显著相关的SNP位点或基因型信息。一般横轴是代表染色体，纵轴则表示SNP位点与表型关联的显著程度，一般-log10(Pvalue)，实线一般表示统计学上的显著性cutoff。

曼哈顿图的实现方法，包括现在的形式也是多种多样，我们先来测试几个最简单的。

第一个：qqman包。

library(qqman)

我们测试数据用的是这个包自带的测试数据。

data(gwasResults)

head(gwasResults)

manhattan(gwasResults, col = rainbow(22), suggestiveline = -log10(1e-05), genomewideline = -log10(5e-08), annotatePval = 5e-08, annotateTop = FALSE)

suggestiveline：一般是显著性cutoff

genomewideline：高置信显著性

annotatePval：If set, SNPs below this p-value will be annotated on the plot。如果设置，那么高于这个Pvalue的SNP将会被标记出来。

annotateTop：If TRUE, only annotates the top hit on each chromosome that is below the annotatePval threshold。如果设置成TRUE，那么只会标记每个染色体上最top的SNP（低于annotatePval）。

第二个：ggplot

library(ggplot2)

library(tidyverse)

因为我们位置是根据染色体变化的，所以我们需要把位置做个转化，转化成X轴上连续的坐标。但是就会出现一个问题，就是X轴上的刻度和坐标，我也需要手动设置，安排好22个染色体的位置。

#计算染色体刻度坐标

gwasResults$SNP1 <- seq(1, nrow(gwasResults), 1)

gwasResults$CHR <- factor(gwasResults$CHR, levels = unique(gwasResults$CHR))

chr <- aggregate(gwasResults$SNP1, by = list(gwasResults$CHR), FUN = median)

aggregate函数是数据处理中常用到的函数。可以按照要求把数据打组聚合，然后对聚合以后的数据进行加和、求平均等各种操作。

aggregate(x, by, FUN, ..., simplify = TRUE, drop = TRUE)

当然也可以用我们前几期常用的通道操作。

#data_new <- gwasResults %>% group_by(CHR) %>% mutate(mean = mean(SNP1))

#data_new2 <- unique(data_new[,c(2,6)])

colnames(chr) <- c("label","location")

下面我们就可以用ggplot来画图了。

ggplot(gwasResults, aes(SNP1, -log10(P)))+

geom_point(aes(color = CHR), show.legend = FALSE) +

scale_color_manual(values = rainbow(22)) +

geom_hline(yintercept = c(-log10(1e-05), -log10(5e-08)), color = c('blue', 'red'), size = 0.35) +

scale_x_continuous(breaks = chr$location, labels = chr$label, expand = c(0, 0)) +

scale_y_continuous(breaks=seq(1,9,2),labels=as.character(seq(1,9,2)),expand = c(0, 0), limits = c(0, 9)) +

theme(panel.grid = element_blank(), axis.line = element_line(colour = 'black'), panel.background = element_rect(fill = 'transparent')) +

labs(x = 'Chromosome', y=expression(-log[10](P)))+

#annotate('rect', xmin = 0, xmax = max(gwasResults$SNP1), ymin = -log10(1e-05), ymax = -log10(5e-08), fill = 'gray98') +  #可以加个矩形框，但是我觉得没啥用

theme(text = element_text(size = 20))

把超过显著性的Pvalue的SNP标签加入。

top <- gwasResults%>% filter(P<=5e-08)

ggplot(gwasResults, aes(SNP1, -log10(P)))+

geom_point(aes(color = CHR), show.legend = FALSE) +

scale_color_manual(values = rainbow(22)) +

geom_hline(yintercept = c(-log10(1e-05), -log10(5e-08)), color = c('blue', 'red'), size = 0.35) +

scale_x_continuous(breaks = chr$location, labels = chr$label, expand = c(0, 0)) +

scale_y_continuous(breaks=seq(1,9,2),labels=as.character(seq(1,9,2)),expand = c(0, 0), limits = c(0, 9)) +

theme(panel.grid = element_blank(), axis.line = element_line(colour = 'black'), panel.background = element_rect(fill = 'transparent')) +

labs(x = 'Chromosome', y = expression(''~-log[10]~'(P)'))+

#annotate('rect', xmin = 0, xmax = max(gwasResults$SNP1), ymin = -log10(1e-05), ymax = -log10(5e-08), fill = 'gray98') +

theme(text = element_text(size = 20))+

geom_label_repel(data=top,aes(x=SNP1, y=-log10(P), label = SNP),size = 4,box.padding = unit(0.5, 'lines'),show.legend = FALSE)

下面我们尝试用曼哈顿图画一个微生物领域常见的图。

一般情况下，在这种图中：

散点：代表单个OTU；

散点大小：相对丰度；

散点形状：显著富集实心圆点，否则为圆环；

灰色背景：间隔每个目水平(或强调是否存在显著富集OTUs)；

水平线：显著性水平p = 0.05；

我们也把上述技巧运用以下，来画下面这个图，这个是我们测试数据的样子。

otu_stat <- read.table("otu_sign.txt",header=T,sep="\t")

#先按分类门水平排序，这里直接按首字母排序了，大家也可以按照自己感兴趣的顺序来排。

otu_stat <- otu_stat[order(otu_stat$phylum), ]

#然后和上面一样，生成连续的在X轴上对应的坐标信息

otu_stat$otu_sort <- 1:nrow(otu_stat)

我们先生成一个最基本的图。

ggplot(otu_stat, aes(otu_sort, -log10(p_value))) +

geom_point(aes(size = abundance, color = phylum, shape = enrich)) +

scale_size(range = c(1, 5))+

theme_bw()+

scale_shape_manual(values=c("Enriched" = 17, "Depleted" = 25, "Non-signficant" = 20))+ #指定enrich类别的shape

theme(panel.grid = element_blank(), axis.line = element_line(colour = 'black'), panel.background = element_rect(fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) +

labs(x="OTUs", y = expression(~-log[10](P)))

下面我们像上面那个曼哈顿图一样手动修改X轴的刻度和注释。

chr <- aggregate(otu_stat$otu_sort, by = list(otu_stat$phylum ), FUN = median)

colnames(chr) <- c("label","location")

计算每个门类的均值位置，也就是label放的位置。

ggplot(otu_stat, aes(otu_sort, -log10(p_value))) +

geom_point(aes(size = abundance, color = phylum, shape = enrich)) +

scale_x_continuous(breaks = chr$location, labels = chr$label, expand = c(0, 0)) +  #手动添加刻度注释，在每个门类的均值位置

scale_size(range = c(1, 5))+

theme_bw()+

scale_shape_manual(values=c("Enriched" = 17, "Depleted" = 25, "Non-signficant" = 20))+

theme(panel.grid = element_blank(), axis.line = element_line(colour = 'black'),

      panel.background = element_rect(fill = 'transparent'),

      axis.text.x = element_text(angle = 45, hjust = 1),  #调整X轴坐标的角度

      axis.text = element_text(face = "bold"),

      legend.key = element_rect(fill = 'transparent')) +

theme(text = element_text(size = 20))+

labs(x="OTUs", y = expression(~-log[10](P)))

ggplot(otu_stat, aes(otu_sort, -log10(p_value))) +

geom_point(aes(size = abundance, color = phylum, shape = enrich)) +

scale_x_continuous(breaks = chr$location, labels = chr$label, expand = c(0, 0)) +

scale_size(range = c(1, 5))+

theme_bw()+

scale_shape_manual(values=c("Enriched" = 17, "Depleted" = 25, "Non-signficant" = 20))+

theme(panel.grid = element_blank(), axis.line = element_line(colour = 'black'),

      panel.background = element_rect(fill = 'transparent'),

      axis.text.x = element_text(angle = 45, hjust = 1),

      axis.text= element_text(face = "bold",size=12),

      axis.title=element_text(face = "bold",size=15),

      legend.title = element_text(size = 12),

      legend.key = element_rect(fill = 'transparent')) +

labs(x = NULL, y = expression(''~-log[10]~'(P)'), size = 'relative abundance (%)', shape = 'significantly enriched') +

guides(color = 'none') +  #把color显示的legend给去除了，因为已经添加X轴坐标了

geom_hline(yintercept = c(-log10(1e-05), -log10(5e-08)), color = c('blue', 'red'), size = 1) +

theme(legend.position = 'top', legend.direction = "horizontal")  #调整legend的位置和方向

我不喜欢添加矩形框，但是如果想要在每个门类那里添加矩形框，也是可以的。我们就要计算矩形框的起始和结束位置。

其实很简单，我们只需要获得每个group中，otu_sort的最大值和最小值即可。比如，举个例子，我们获得Acidobacteria数据集。

x<- otu_stat %>% filter(phylum=="Acidobacteria")

然后利用min(x$otu_sort)，max(x$otu_sort)即可得到最大或者最小值。

p<-ggplot(otu_stat, aes(otu_sort, -log10(p_value))) +

geom_point(aes(size = abundance, color = phylum, shape = enrich)) +

scale_x_continuous(breaks = chr$location, labels = chr$label, expand = c(0, 0)) +

scale_size(range = c(1, 5))+

theme_bw()+

scale_shape_manual(values=c("Enriched" = 17, "Depleted" = 25, "Non-signficant" = 20))+

theme(panel.grid = element_blank(), axis.line = element_line(colour = 'black'),

      panel.background = element_rect(fill = 'transparent'),

      axis.text.x = element_text(angle = 45, hjust = 1),

      axis.text= element_text(face = "bold",size=12),

      axis.title=element_text(face = "bold",size=15),

      legend.title = element_text(size = 12),

      legend.key = element_rect(fill = 'transparent')) +

labs(x = NULL, y = expression(''~-log[10]~'(P)'), size = 'relative abundance (%)', shape = 'significantly enriched') +

guides(color = 'none') +

geom_hline(yintercept = c(-log10(1e-05), -log10(5e-08)), color = c('blue', 'red'), size = 1) +

theme(legend.position = 'top', legend.direction = "horizontal")

把前面的图保存在p变量中。

p+

annotate('rect', xmin =min(x$otu_sort), xmax=max(x$otu_sort), ymin = -Inf, ymax = Inf, alpha=0.5,

fill ="gray85")

这样我们就添加了一个灰色框。

接下来，我们写个for循环既可以实现。

group <- unique(otu_stat$phylum)

for(i in 1:length(group))

{

  x<- otu_stat %>% filter(phylum==group[i])

  p<- p+annotate('rect', xmin =min(x$otu_sort), xmax=max(x$otu_sort), ymin = -Inf, ymax = Inf, alpha=0.5,

  fill = ifelse(i %% 2 == 0, 'gray95', 'gray85'))

}

p