曼哈頓圖，我自己最常用于的場景就是GWAS分析之后，常使用它展示基因組中與某種表型顯著相關(guān)的SNP位點或基因型信息。一般橫軸是代表染色體，縱軸則表示SNP位點與表型關(guān)聯(lián)的顯著程度，一般-log10(Pvalue)，實線一般表示統(tǒng)計學(xué)上的顯著性cutoff。
曼哈頓圖的實現(xiàn)方法，包括現(xiàn)在的形式也是多種多樣，我們先來測試幾個最簡單的。
第一個：qqman包。
library(qqman)
我們測試數(shù)據(jù)用的是這個包自帶的測試數(shù)據(jù)。
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。如果設(shè)置，那么高于這個Pvalue的SNP將會被標(biāo)記出來。
annotateTop：If TRUE, only annotates the top hit on each chromosome that is below the annotatePval threshold。如果設(shè)置成TRUE，那么只會標(biāo)記每個染色體上最top的SNP（低于annotatePval）。
第二個：ggplot
library(ggplot2)
library(tidyverse)
因為我們位置是根據(jù)染色體變化的，所以我們需要把位置做個轉(zhuǎn)化，轉(zhuǎn)化成X軸上連續(xù)的坐標(biāo)。但是就會出現(xiàn)一個問題，就是X軸上的刻度和坐標(biāo)，我也需要手動設(shè)置，安排好22個染色體的位置。
#計算染色體刻度坐標(biāo)
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函數(shù)是數(shù)據(jù)處理中常用到的函數(shù)?？梢园凑找蟀褦?shù)據(jù)打組聚合，然后對聚合以后的數(shù)據(jù)進行加和、求平均等各種操作。
aggregate(x, by, FUN, ..., simplify = TRUE, drop = TRUE)
當(dāng)然也可以用我們前幾期常用的通道操作。
#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標(biāo)簽加入。
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)
下面我們嘗試用曼哈頓圖畫一個微生物領(lǐng)域常見的圖。
一般情況下，在這種圖中：
散點：代表單個OTU；
散點大小：相對豐度；
散點形狀：顯著富集實心圓點，否則為圓環(huán)；
灰色背景：間隔每個目水平(或強調(diào)是否存在顯著富集OTUs)；
水平線：顯著性水平p = 0.05；
我們也把上述技巧運用以下，來畫下面這個圖，這個是我們測試數(shù)據(jù)的樣子。
otu_stat <- read.table("otu_sign.txt",header=T,sep="\t")
#先按分類門水平排序，這里直接按首字母排序了，大家也可以按照自己感興趣的順序來排。
otu_stat <- otu_stat[order(otu_stat$phylum), ]
#然后和上面一樣，生成連續(xù)的在X軸上對應(yīng)的坐標(biāo)信息
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),? #調(diào)整X軸坐標(biāo)的角度
? ? ? 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給去除了，因為已經(jīng)添加X軸坐標(biāo)了
geom_hline(yintercept = c(-log10(1e-05), -log10(5e-08)), color = c('blue', 'red'), size = 1) +
theme(legend.position = 'top', legend.direction = "horizontal")? #調(diào)整legend的位置和方向
我不喜歡添加矩形框，但是如果想要在每個門類那里添加矩形框，也是可以的。我們就要計算矩形框的起始和結(jié)束位置。
其實很簡單，我們只需要獲得每個group中，otu_sort的最大值和最小值即可。比如，舉個例子，我們獲得Acidobacteria數(shù)據(jù)集。
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循環(huán)既可以實現(xiàn)。
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
