曼哈顿图，我自己最常用于的场景就是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
网友评论