Skip to content

Instantly share code, notes, and snippets.

@alswl
Last active December 21, 2015 06:09
Show Gist options
  • Save alswl/6262039 to your computer and use it in GitHub Desktop.
Save alswl/6262039 to your computer and use it in GitHub Desktop.
# 15. 等级丰度曲线
# 根据 name 返回 color, lty
get_line_type = function(name) {
if (length(grep('-1$', name)) > 0) {
color <- 1
} else if (length(grep('-2$', name)) > 0) {
color <- 2
} else if (length(grep('-3$', name)) > 0) {
color <- 3
}
if (length(grep('^0-', name)) > 0) {
lty <- 1
} else if (length(grep('^2-', name)) > 0) {
lty <- 2
} else if (length(grep('^3-', name)) > 0) {
lty <- 3
} else if (length(grep('^8-', name)) > 0) {
lty <- 4
} else if (length(grep('^13-', name)) > 0) {
lty <- 5
} else if (length(grep('^15-', name)) > 0) {
lty <- 6
}
c(color, lty)
}
# 1. 准备数据
table <- read.table(file="1.csv", sep="\t", head=T, check.names=F) # 打开 CSV
rownames(table) <- table[, 1] # 将第一行指定为表格 Header
t_data <- table[, -1] # 剩下的是数据
abud <- lapply(1 : ncol(t_data), # 排序,去除不存在的
function(x) {
sort(t_data[t_data[, x] > 0, x],
decreasing=T) * 100 / sum(t_data[, x])
})
names(abud) <- colnames(t_data) # 设置名称
tips <- sapply(abud, length) # 获取长度
# 2. 准备画板
pdf("Rank-Abudance-Curves.pdf", width=20, height=15)
par(mar=c(10, 10, 10, 10))
plot(x=c(0, 0), y=c(1, 1), font.lab=2, font.main=2, cex.lab=2,
cex.main=2, main="Rank-abundance distribution curve", xlab="OTURank",
ylab="RelativeAbundance", xlim=c(1, max(tips)), ylim=c(-3, 2), xaxs="i",
yaxs="i", xaxt="s", yaxt="n", las=1, cex.axis=2) # 背景框架
axis(side=2, pos=c(0, 0), at=-3:2, labels=as.character(10^(-3:2)), cex.axis=2,
las=1) # 坐标
legend("topright", legend=names(tips),
col=sapply(lapply(names(tips), get_line_type), function(x) x[[1]]),
lty=sapply(lapply(names(tips), get_line_type), function(x) x[[2]]),
cex=2, bty="n",
ncol=2)
# 3. 画图
for(i in 1 : length(abud)) { # 画图
name <- names(tips)[i]
if(length(abud[[i]])> 0) {
x = 1 : length(abud[[i]])
y = abud[[i]]
y = log10(y)
points(x, y, type="l",
col=get_line_type(name)[1],
lty=get_line_type(name)[2])
}
}
# 4. 保存 PDF
dev.off()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment