r15.R

# 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()