1

我正在尝试自动标记manhattan plot.

对于那些不知道 amanhattan plot是什么的人来说，这真的没关系。

重写旧式 R 代码ggplot2对我来说似乎是一个更大的挑战。

由于我想添加到脚本的唯一功能是能够很好地自动标记一些数据点（使用ggrepel），所以我想也许我可以将标签覆盖在用旧代码绘制的图表上。

我目前的尝试如下：

library(ggplot2);
library(ggrepel);

d=read.table("a.txt",header=T,fill=T, sep=" ");
dmht='';

dmht<-data.frame(chrom=d[,9], txStart = d[,11], "-log10(PValue)" = -log10(d[,5]))
# sort it
o<-order(dmht[,1],dmht[,2]);
dmht<-data.frame(dmht[o,]);
names(dmht)<-c("chrom", "txStart", "-log10(PValue)");
attach(dmht);
chrs<-c('chr1','chr2','chr3','chr4','chr5','chr6','chr7','chr8','chr9','chr10','chr11','chr12','chr13','chr14','chr15','chr16','chr17','chr18','chr19','chrX','chrY');
chrLabels = c('1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14', '15', '16', '17', '18', '19', 'X', 'Y');
chrLen<-rep(0,length(chrs));
chrMin<-rep(0,length(chrs));
chrMax<-rep(0,length(chrs));
totalLen = 0;

for(i in 1:length(chrs)) {
    dchr<-subset(dmht, chrom == chrs[i]);
    chrMin[i]<-min(as.numeric(dchr[,2]),na.rm=T);
    chrMax[i]<-max(as.numeric(dchr[,2]),na.rm=T);
    chrLen[i]<-chrMax[i] - chrMin[i] + 1;
    totalLen = totalLen + chrMax[i];
}


ds=read.table("selected2label.txt",header=T,fill=T,sep="    ");
dsel='';

dsel<-data.frame(ds[,9], ds[,11], -log10(ds[,5]) , ds[,1])
ymin = min(as.numeric(dmht[,3], dsel[,3]), na.rm=T);
ymax = max(as.numeric(dmht[,3], dsel[,3]), na.rm=T);
# chr start positions
chrStart<-rep(0,length(chrs));

for(i in 1:length(chrs)) {
    if(i == 1) {
    chrStart[i] = 1;
    }
    else {
    chrStart[i] = chrStart[i-1] + chrMax[i-1] + 1;
    }
}

#dmht=subset(dmht,dmht[,1] < 24); # remove mitochondra snps
colors <- rep(c("blue", "green", "cyan"),9);

png("result/test/mhtTest.png" , width=1600);
par(las=3, lab=c(length(chrLabels),5,7))

# draw the dots
dchr<-subset(dmht, chrom == chrs[1]);
plot(as.numeric(dchr[,2])+chrStart[1], as.numeric(dchr[,3]), col=colors[1], ylim=c(ymin,ymax),xlim=c(chrMin[1],totalLen),axes=F,ylab="-log10(PValue)", xlab="Chromosome", main="");
for(i in 2:length(chrs)) {
    dchr<-subset(dmht, chrom == chrs[i]);
    points(as.numeric(dchr[,2])+chrStart[i], as.numeric(dchr[,3]), col=colors[i]);
}

axis(side=1,labels=chrLabels,at=chrStart);
axis(side=2);

# draw the quantiles
quants<-quantile(as.numeric(dmht[,3]), p=c(),na.rm=T);

for(q in quants) {
    abline(h=q);
}

# draw the abs values
abss<-c();
for(a in abss) {
    abline(h=a);
}

# sort it
#os<-order(dsel[,1],dsel[,2]);
#dsel<-data.frame(dsel[os,]);
#dsel
#dsel<-data.frame(dsel);
#dsel
colnames(dsel)<-c("chrom", "txStart", "-log10(PValue)" , 'GENE_ID');
detach(dmht);
attach(dsel);

# highlight the selected dots
for(i in 1:length(chrs)) {
    dchr<-subset(dsel, chrom == chrs[i] | paste("chr", chrom, sep="") == chrs[i]);
    if(length(dchr[,2]) > 0) {
    print(dchr)
    # this is the new code
    geom_label_repel(data = dchr, aes(label = GENE_ID, x=as.numeric(txStart)+chrStart[i], y = as.numeric(dchr[,3])), size = 5, box.padding = unit(0.35, "lines"), point.padding = unit(0.5, "lines")) 
    # replacing the line in the old script:
    # text(as.numeric(dchr[,2])+chrStart[i], as.numeric(dchr[,3]), dchr[,4])
    }
}


dev.off();

我替换的唯一行（除了添加libary(ggplot2)& library(ggrepel)）是：

    # this is the new code
    geom_label_repel(data = dchr, aes(label = GENE_ID, x=as.numeric(txStart)+chrStart[i], y = as.numeric(dchr[,3])), size = 5, box.padding = unit(0.35, "lines"), point.padding = unit(0.5, "lines")) 
    # replacing the line in the old script:
    # text(as.numeric(dchr[,2])+chrStart[i], as.numeric(dchr[,3]), dchr[,4])

问题是标签根本没有出现。

您的帮助将不胜感激。

4

2 回答 2

0

如果您尝试将 ggrepel 与 base-R 图一起使用，简短的回答是您不能。ggplot2 包是 ggrepel 的一个附加组件，它基于网格图形系统（就像点阵绘图包一样）。无法将 ggplot 的元素与基本 R 图（您将其描述为旧式代码）混合和匹配。

于 2018-01-04T14:42:51.513 回答
0

这段代码可能有效，也可能无效，我只是盯着你的情节代码。

ggplot(aes(x=as.numeric(dchr[,2]) + chrStart[1], y=as.numeric(dchr[,3]), color = dchr[, 2]) + 
    geom_point() +
    ylim(c(ymin,ymax)) +
    xlim(c(chrMin[1],totalLen)) +
    ylab("-log10(PValue)") +
    xlab("Chromosome")+
    ggtitle("") + 
        geom_label_repel(data = dchr, inherit.aes=FALSE,
        aes(label = GENE_ID, x=as.numeric(txStart)+chrStart[i], y = as.numeric(dchr[,3])), size = 5, box.padding = unit(0.35, "lines"), point.padding = unit(0.5, "lines"))
于 2018-01-05T02:53:21.180 回答