R語言學(xué)習(xí)之線圖的繪制詳解

更新時(shí)間：2022年03月01日 16:01:39 作者：陳同

線圖是反映趨勢變化的一種方式，其輸入數(shù)據(jù)一般也是一個(gè)矩陣。本文將利用R語言繪制單線圖、多線圖以及橫軸文本線圖，感興趣的可以了解一下

線圖

線圖是反映趨勢變化的一種方式，其輸入數(shù)據(jù)一般也是一個(gè)矩陣。

單線圖

假設(shè)有這么一個(gè)矩陣，第一列為轉(zhuǎn)錄起始位點(diǎn)及其上下游5 kb的區(qū)域，第二列為H3K27ac修飾在這些區(qū)域的豐度，想繪制一張線圖展示。

profile="Pos;H3K27ac
-5000;8.7
-4000;8.4
-3000;8.3
-2000;7.2
-1000;3.6
0;3.6
1000;7.1
2000;8.2
3000;8.4
4000;8.5
5000;8.5"

讀入數(shù)據(jù) (經(jīng)過前面幾篇的聯(lián)系，這應(yīng)該都很熟了)

profile_text <- read.table(text=profile, header=T, row.names=1, quote="",sep=";")
profile_text

 ? ? ?H3K27ac
-5000 ? ? 8.7
-4000 ? ? 8.4
-3000 ? ? 8.3
-2000 ? ? 7.2
-1000 ? ? 3.6
0 ? ? ? ? 3.6
1000 ? ? ?7.1
2000 ? ? ?8.2
3000 ? ? ?8.4
4000 ? ? ?8.5
5000 ? ? ?8.5

# 在melt時(shí)保留位置信息
# melt格式是ggplot2畫圖最喜歡的格式
# 好好體會下這個(gè)格式，雖然多占用了不少空間，但是確實(shí)很方便

# 這里可以用 `xvariable`，也可以是其它字符串，但需要保證后面與這里的一致
# 因?yàn)檫@一列是要在X軸顯示，所以起名為`xvariable`。
profile_text$xvariable = rownames(profile_text)
library(ggplot2)
library(reshape2)
data_m <- melt(profile_text, id.vars=c("xvariable"))
data_m

 ? xvariable variable value
1 ? ? ?-5000 ?H3K27ac ? 8.7
2 ? ? ?-4000 ?H3K27ac ? 8.4
3 ? ? ?-3000 ?H3K27ac ? 8.3
4 ? ? ?-2000 ?H3K27ac ? 7.2
5 ? ? ?-1000 ?H3K27ac ? 3.6
6 ? ? ? ? ?0 ?H3K27ac ? 3.6
7 ? ? ? 1000 ?H3K27ac ? 7.1
8 ? ? ? 2000 ?H3K27ac ? 8.2
9 ? ? ? 3000 ?H3K27ac ? 8.4
10 ? ? ?4000 ?H3K27ac ? 8.5
11 ? ? ?5000 ?H3K27ac ? 8.5

然后開始畫圖，與上面畫heatmap一樣。

# variable和value為矩陣melt后的兩列的名字，內(nèi)部變量, variable代表了點(diǎn)線的屬性，value代表對應(yīng)的值。
p <- ggplot(data_m, aes(x=xvariable, y=value),color=variable) + geom_line()
p
# 圖會存儲在當(dāng)前目錄的Rplots.pdf文件中，如果用Rstudio，可以不運(yùn)行dev.off()
dev.off()

滿心期待一個(gè)倒鐘形曲線，結(jié)果，

什么也沒有。

仔細(xì)看，出來一段提示

geom_path: Each group consists of only one observation. 
Do you need to adjust the group aesthetic?

原來默認(rèn)ggplot2把每個(gè)點(diǎn)都視作了一個(gè)分組，什么都沒畫出來。而data_m中的數(shù)據(jù)都來源于一個(gè)分組H3K27ac，分組的名字為variable，修改下腳本，看看效果。

p <- ggplot(data_m, aes(x=xvariable, y=value,color=variable,group=variable)) + 
 ? ? geom_line() + theme(legend.position=c(0.1,0.9))
p
dev.off()

圖出來了，一條線，看一眼沒問題；再仔細(xì)看，不對了，怎么還不是倒鐘形，原來橫坐標(biāo)錯(cuò)位了。

檢查下數(shù)據(jù)格式

summary(data_m)

 ?xvariable ? ? ?variable ? ? ? ?
 Length:11 ? ? ? H3K27ac:11 ? ??
 Class :character ? ? ? ? ? ??
 Mode ?:character ? ? ? ?

問題來了，xvariable雖然看上去數(shù)字，但存儲的實(shí)際是字符串 (因?yàn)槭亲鳛樾忻肿x取的)，需要轉(zhuǎn)換為數(shù)字。

data_m$xvariable <- as.numeric(data_m$xvariable)

#再檢驗(yàn)下
is.numeric(data_m$xvariable)

[1] TRUE

好了，繼續(xù)畫圖。

# 注意斷行時(shí)，加號在行尾，不能放在行首
p <- ggplot(data_m, aes(x=xvariable, y=value,color=variable,group=variable)) +
 ? ? geom_line() + theme(legend.position=c(0.1,0.8))
p
dev.off()

圖終于出來了，調(diào)了下legend的位置，看上去有點(diǎn)意思了。

有點(diǎn)難看，如果平滑下，會不會好一些，stat_smooth可以對繪制的線進(jìn)行局部擬合。在不影響變化趨勢的情況下，可以使用 (但慎用)。

p <- ggplot(data_m, aes(x=xvariable, y=value,color=variable,group=variable)) + 
 ? ? geom_line() + stat_smooth(method="auto", se=FALSE) + 
 ? ? theme(legend.position=c(0.1,0.8))
p
dev.off()

從圖中看，趨勢還是一致的，線條更優(yōu)美了。另外一個(gè)方式是增加區(qū)間的數(shù)量，線也會好些，而且更真實(shí)。

stat_smooth和geom_line各繪制了一條線，只保留一條就好。

p <- ggplot(data_m, aes(x=xvariable, y=value,color=variable,group=variable)) + 
 ? ? stat_smooth(method="auto", se=FALSE) + theme(legend.position=c(0.1,0.8))
p
dev.off()

好了，終于完成了單條線圖的繪制。

多線圖

那么再來一個(gè)多線圖的例子吧，只要給之前的數(shù)據(jù)矩陣多加幾列就好了。

profile = "Pos;h3k27ac;ctcf;enhancer;h3k4me3;polII
-5000;8.7;10.7;11.7;10;8.3
-4000;8.4;10.8;11.8;9.8;7.8
-3000;8.3;10.5;12.2;9.4;7
-2000;7.2;10.9;12.7;8.4;4.8
-1000;3.6;8.5;12.8;4.8;1.3
0;3.6;8.5;13.4;5.2;1.5
1000;7.1;10.9;12.4;8.1;4.9
2000;8.2;10.7;12.4;9.5;7.7
3000;8.4;10.4;12;9.8;7.9
4000;8.5;10.6;11.7;9.7;8.2
5000;8.5;10.6;11.7;10;8.2"

profile_text <- read.table(text=profile, header=T, row.names=1, quote="",sep=";")

profile_text$xvariable = rownames(profile_text)
data_m <- melt(profile_text, id.vars=c("xvariable"))
data_m$xvariable <- as.numeric(data_m$xvariable)

# 這里group=variable，而不是group=1 (如果上面你用的是1的話)
# variable和value為矩陣melt后的兩列的名字，內(nèi)部變量, variable代表了點(diǎn)線的屬性，value代表對應(yīng)的值。
p <- ggplot(data_m, aes(x=xvariable, y=value,color=variable,group=variable)) + 
 ? ? stat_smooth(method="auto", se=FALSE) + theme(legend.position=c(0.85,0.2))
p
dev.off()

橫軸文本線圖

如果橫軸是文本，又該怎么調(diào)整順序呢？還記得之前熱圖旁的行或列的順序調(diào)整嗎？重新設(shè)置變量的factor水平就可以控制其順序。

profile = "Pos;h3k27ac;ctcf;enhancer;h3k4me3;polII
-5000;8.7;10.7;11.7;10;8.3
-4000;8.4;10.8;11.8;9.8;7.8
-3000;8.3;10.5;12.2;9.4;7
-2000;7.2;10.9;12.7;8.4;4.8
-1000;3.6;8.5;12.8;4.8;1.3
0;3.6;8.5;13.4;5.2;1.5
1000;7.1;10.9;12.4;8.1;4.9
2000;8.2;10.7;12.4;9.5;7.7
3000;8.4;10.4;12;9.8;7.9
4000;8.5;10.6;11.7;9.7;8.2
5000;8.5;10.6;11.7;10;8.2"

profile_text <- read.table(text=profile, header=T, row.names=1, quote="",sep=";")

profile_text_rownames <- row.names(profile_text)

profile_text$xvariable = rownames(profile_text)
data_m <- melt(profile_text, id.vars=c("xvariable"))

# 就是這一句，會經(jīng)常用到
data_m$xvariable <- factor(data_m$xvariable, levels=profile_text_rownames, ordered=T)

# geom_line設(shè)置線的粗細(xì)和透明度
p <- ggplot(data_m, aes(x=xvariable, y=value,color=variable,group=variable)) + 
 ? ? geom_line(size=1, alpha=0.9) + theme(legend.position=c(0.85,0.2)) +
 ? ? theme(axis.text.x=element_text(angle=45,hjust=1, vjust=1))

# stat_smooth
#p <- ggplot(data_m, aes(x=xvariable, y=value,color=variable,group=variable)) + 
# ? ? stat_smooth(method="auto", se=FALSE) + theme(legend.position=c(0.85,0.2)) +
# ? ? theme(axis.text.x=element_text(angle=45,hjust=1, vjust=1))

p
dev.off()

比較下位置信息做為數(shù)字(前面的線圖)和位置信息橫軸的差別。當(dāng)為數(shù)值時(shí)，ggplot2會選擇合適的幾個(gè)刻度做標(biāo)記，當(dāng)為文本時(shí)，會全部標(biāo)記。另外文本橫軸，smooth效果不明顯 (下面第2張圖)。