一、前言
好难受,今天的R崩了、昨天才安装的包怎么报错了、为啥我装了这个包没反应,什么鬼怎么又乱码了。这玩泥巴(╯‵□′)╯︵┻━┻
不行要佛系生活,好久没写简书了,还是来整理下代码吧┬—┬ノ('-'ノ)
二、常见的生信分析代码整理
今天要说的事如题,这些都是本人完成的流程中非常常见的R代码,以下一一整理并为其注释,没事就经常翻阅这些代码瞧瞧。
不过话说回来,R真的好垃圾啊(/‵Д′)/~ ╧╧
老板来了,我永远喜欢R语言┬─┬ ノ( ' - 'ノ)
1.系统读写相关
设置工作路径
设置工作路径,大家习惯打开R程序,而不是在固定的位置打开R,那么很有可能你就要设置工作路径,这样做的目的是为了节省语句。注意一下,这是/左斜杠,windows下复制路径是\斜杠,新人老容易犯这个错误
setwd("D:/R")
读取文件
读取txt文件,./为当前路径,与OS文件读取方法类似,参数sep = "\t"为以TAB制表符区分行列,comment.char = "#"为#号为注释内容,stringsAsFactors = F 为取消字符串转换为因子,header = T为有表头,fill=TRUE 如果为true,则在行长度不相等的情况下,将隐式添加空白字段。./是指当前目录下的Human_ppi_symbol_2019_7.txt文件,当然也可以写绝对路径。
a <- read.table("./Human_ppi_symbol_2019_7.txt",sep = "\t",comment.char = "#", stringsAsFactors = F,header = T, fill=TRUE)
写入文件
写入txt文件,不指定路径则为当前路径下创建,学会举一反三,参数有很多,这里并不能写出所有,因此建议大家用?write.table查阅文档,当你不懂这个函数有什么参数的时候,就强烈建议?+函数名来查阅文档,我就是这么学习的。
write.table(Gene, 'Gene.txt',sep = '\t', quote=F, row.names=F)
2.数据筛选相关
挑选特定的行和列
很多时候我们需要筛选数据,最常见的比如我希望得到某一列的固定逻辑的数据(我这里以大于900的数据为例),写法是这样的:
筛选大于900分的行,以a[x,y]为基础,a$为catch住含有名为“combined_score”的列,>900为筛选逻辑,前提要求数据类型为数字,得到数据
这句话应该这样解析,a$combined_score > 900是逻辑值,返回True和Flase,你也可以理解为返回了对应的行数,那么我们只要把这个行数再赋值给a[x,y]中的x,那么得到的就是含有大于900的所有a 的x行。
a = a[a$combined_score > 900,]
取消科学计数法
取消科学计数法,有一些基因表达量太小或者p.Value太小,会用科学技术进行计算,担心后面筛除的时候会丢失
options(scipen = 400)
筛选同时满足两个要求的数据
筛选同时满足P.Value < 0.05 和 logFC > 0.5的行,因为logFC有负数,可用abs绝对值化,这句话及上句都常用在做出来的DEG(差异表达基因)中
DEG_sort_p <- DEG_sort_p[DEG_sort_p$P.Value < 0.05 & abs(DEG_sort_p$logFC) > 0.5,]
常用语句
统计行列
dim(变量名)
行名称更换
row.names(DEG_sort_p) = DEG_sort_p[,1]
把DEG_sort_p中的差异基因提出来,制作成一个data.frame
diff_Gene = data.frame(Gene=rownames(DEG_sort_p))
这里是蛋白互作的关键合并步骤,第①步,把基因中含有该蛋白的“protein1“merge得到
Gene1=merge(diff_Gene, String_hsa_PPIs_Symbol1, by.x = "Gene",by.y = "protein1",all=FALSE)
判断蛋白基因是否在全在基因里,返回真假值为True
table(Gene$Gene%in%DEG_sort_p$DEG)
unique是去掉基因名字重复的基因,只保留一个,保留有效基因
Gene=unique(data.frame(Gene=Gene[,2]))
判断蛋白基因是否在全在基因里,返回真假值210122为True
table(Gene$Gene%in%DEG_sort_p$DEG)
unique是去掉基因名字重复的基因,只保留一个,保留有效基因
Gene=unique(data.frame(Gene=Gene[,2]))
简单粗暴的合并行
Gene=rbind(Gene,diff_Gene)
把重复基因的行去掉
Gene=unique(Gene)
重命名行
colnames(data2)=c("Gene1","Gene2","score")
把data2的第一列变成character格式,这样做的意义是有一些包要求输入进入格式是这种数据类型,
as.array()、as.data.frame()、as.matrix()、as.numeric()
data2$Gene1 <-as.character(data2[,1])
自创代码:阅读当前目录下所有文件,如果里面包含GSE名称的文件夹,则提取出,进而阅读里面包含exp_matrix的表达矩阵并进行merge合并,要求这些GSE都是基于相同平台。意思是,当前目录下有N个GSE文件夹,这些文件夹里面有各自已经完成好的表达矩阵,把这些表达矩阵都取出来并进行merge合并成为一个大矩阵,适合做数据合并用。
dir = list.dirs()
# 筛出路径中包含GSE的文件夹
GSE_dir = dir[grep("GSE",dir)]
# 循环获取目录,要求根目录下必须含有GSE的文件夹,GSE文件夹内含有exp_matrix名称的表达矩阵,要求表达矩阵是相同平台,输出变量名为GSE_all
for (index in 1:length(GSE_dir)){
if (index == 1){
GSE_all <- read.table(paste(GSE_dir[index],"/exp_matrix.txt",sep=""), sep = "\t",stringsAsFactors = F, header = T, fill=TRUE, quote = "")
GSE_name = strsplit(GSE_dir[index],"./", fixed = FALSE, perl = FALSE, useBytes = FALSE)
print(GSE_name[[1]][2])
}
if (index >= 2){
GSE_data <- read.table(paste(GSE_dir[index],"/exp_matrix.txt",sep=""), sep = "\t",stringsAsFactors = F, header = T, fill=TRUE, quote = "")
GSE_all = merge(GSE_all,GSE_data,by.x = "Symbol",by.y = "Symbol")
GSE_name = strsplit(GSE_dir[index],"./", fixed = FALSE, perl = FALSE, useBytes = FALSE)
print(GSE_name[[1]][2])
}
}
优雅的增加一行,把行名称加回去
GSEall_DEG_sort_p$Symbol = rownames(GSEall_DEG_sort_p)
求平均的函数,敲重要
meanfun <- function(x) {
x1 <- data.frame(unique(x[,1]))
colnames(x1) <- c("Symbol")
for (i in 2:ncol(x)){
x2 <- data.frame(tapply(x[,i],x[,1],mean))
x2[,2] <- rownames(x2)
colnames(x2) <- c(colnames(x)[i], "Symbol")
x1 <- merge(x1,x2,by.x = "Symbol",by.y = "Symbol",all=FALSE)
}
return(x1)
}