接触R语言已经一年多了,期间都是断断续续的学习。书买了不少,看完的没多少。毕竟,Google在手,一搜就来。
for循环
为什么要写for循环呢?因为我一堆数据需要算很多的回归方程,如果一个一个计算的话,那我估计两三天就没了。碰巧的是,我刚好看到宏基因组推送的代码
),里面有几个我看起来及其复杂的if循环,我就想我是不是可以用for循环来解决我的问题呢?说干就干。想到这个主意的时候已经十点半了,合上电脑,赶回宿舍,匆匆洗漱完,打开电脑撸代码。不算顺利吧,做到十二点五十多完成了for循环的基本构架,能够批量计算回归方程和绘图了。满意的爬上床睡觉了。突然想到不行啊,我还是得把回归方程抄一遍啊,整不成。第二天早上赶到实验室,匆匆改了代码,期间失败很多次,总是无法批量把回归方程写进矩阵(数据框)中。突然一下子就想明白是某个变量的问题了。一改,成了。哈哈。贴上代码,留个念想,鼓励自己。
#载入绘图相关的包
library(ggplot2)
#数据载入
raw_data <- read.table("raw_data.txt", header = T, sep = "\t")
#创建矩阵
matrix <- matrix(0,ncol(raw_data),4)
#for循环
for (i in 2:ncol(raw_data)) {
conc <- raw_data$conc
amino_acid <- raw_data[i]
data <- data.frame(x = conc, y = amino_acid)
colnames(data) <- c("x","y")
amino_acid_name <- colnames(raw_data[i])
fit <- lm(data$y ~ data$x)
slope <- round(fit$coefficients[2],4)
intercept <- round(fit$coefficients[1],4)
info <- summary(fit)
maintitle <- paste(amino_acid_name,"\n","y=",slope,"x+",intercept, "\n","R=", round(sqrt(info$r.squared),4))
regression_equation <- paste("y=",slope,"x+",intercept)
xlab <- paste(amino_acid_name," Concentration(mM/L)")
ylab <- ("Peak Area")
ggplot(data = data, mapping = aes(x = data$x, y = data$y))+
geom_point(position = "identity")+
geom_smooth(method = lm, se = FALSE,color = "black")+
annotate("text", x = max(raw_data[1])-0.025, y = max(raw_data[i])+0.05,label = maintitle)+
theme_bw()+
theme(panel.grid = element_blank(),
plot.title = element_text(hjust = 0.5),
axis.line = element_line(linetype = "solid"),
#axis.ticks = element_line(size = 1,colour = "black"),
axis.ticks = element_blank(),
plot.background = element_rect(linetype = "solid"),
panel.border = element_rect(colour = "black", size = 1),
axis.text = element_text(colour = "black"),
axis.text.x = element_text(colour = "black"),
axis.text.y = element_text(colour = "black")) +
labs(x = xlab, y = ylab, title = " ")
FileName <- paste(amino_acid_name,"_regression_equation", ".pdf", sep = "_")
ggsave(FileName, width = 8, height = 8)
matrix[i-1,1] <- amino_acid_name
matrix[i-1,2] <- regression_equation
matrix[i-1,3] <- paste("r^2=", round(info$r.squared,4))
matrix[i-1,4] <- paste("R=", round(sqrt(info$r.squared),4))
if (i == ncol(raw_data)) {
write.csv(matrix,file = "regrression equation.csv")
}
}
这个代码看起来是臃肿的,但是是我的第一个for循环,我喜欢,哈哈哈,慢慢改进。
样品分布展示
出去办点事情,地铁上收到老师的一张图片:
然后,想让我实现图中的样品分布地图展示。
一看,我*,这图肯定是R实现的,那一定会有现成的R包。瞬间想到之前中科院遗传发育所白洋老师课题组在Nature Biotechnology上的文章【2】,那里面是一定有代码和原始数据的。一找,果然,找到了。
火速下载,哈哈,拿到原始数据,开始上手撸代码:
data <- read.table("2.txt",header = T, sep = "\t")
mp <- NULL
mapworld <- borders("world",colour = "gray50", fill = "cornsilk")
mp <- ggplot() + mapworld + ylim(-90,90) + xlim(-180,180)
map <- mp + geom_point(aes(x = data$Longitude, y = data$Latitude, color = data$Subspecies))+
scale_size(range = c(2, 9))+
ggtitle("lixiang")+
theme(plot.subtitle = element_text(vjust = 1),
plot.caption = element_text(vjust = 1),
panel.grid.major = element_line(linetype = "blank"),
panel.grid.minor = element_line(linetype = "blank"),
panel.background = element_rect(fill = NA))
map
图是什么亚子的呢?咻咻咻→
xiangxiang的丑图
那他们文章里的图片是啥样的呢?
嚯嚯嚯→
NBT文章 Figure 1a【2】
差点。明天改进。
批量化处理是计算机的优势,也应该成为解决问题的一种思考方式。
多看,多想,多记,多练,多试。
下班,回寝。
参考文献
【1】Wang H, Xu X, Vieira FG, et al. The Power of Inbreeding: NGS-Based GWAS of Rice Reveals Convergent Evolution during Rice Domestication. Mol Plant. 2016;9(7):975-85
【2】 Zhang J, Liu YX, Zhang N, et al. NRT1.1B is associated with root microbiota composition and nitrogen use in field-grown rice. Nat Biotechnol. 2019;37(6):676-684
网友评论