一、前言
好難受凯亮,今天的R崩了、昨天才安裝的包怎么報錯了哄尔、為啥我裝了這個包沒反應(yīng)假消,什么鬼怎么又亂碼了。這玩泥巴(╯‵□′)╯︵┻━┻
不行要佛系生活岭接,好久沒寫簡書了富拗,還是來整理下代碼吧┬—┬ノ('-'ノ)
二、常見的生信分析代碼整理
今天要說的事如題鸣戴,這些都是本人完成的流程中非常常見的R代碼啃沪,以下一一整理并為其注釋,沒事就經(jīng)常翻閱這些代碼瞧瞧窄锅。
不過話說回來创千,R真的好垃圾啊(/‵Д′)/~ ╧╧
老板來了,我永遠(yuǎn)喜歡R語言┬─┬ ノ( ' - 'ノ)
1.系統(tǒng)讀寫相關(guān)
設(shè)置工作路徑
設(shè)置工作路徑,大家習(xí)慣打開R程序追驴,而不是在固定的位置打開R械哟,那么很有可能你就要設(shè)置工作路徑,這樣做的目的是為了節(jié)省語句殿雪。注意一下暇咆,這是/左斜杠,windows下復(fù)制路徑是\斜杠丙曙,新人老容易犯這個錯誤
setwd("D:/R")
讀取文件
讀取txt文件爸业,./為當(dāng)前路徑,與OS文件讀取方法類似亏镰,參數(shù)sep = "\t"為以TAB制表符區(qū)分行列扯旷,comment.char = "#"為#號為注釋內(nèi)容,stringsAsFactors = F 為取消字符串轉(zhuǎn)換為因子拆挥,header = T為有表頭薄霜,fill=TRUE 如果為true,則在行長度不相等的情況下纸兔,將隱式添加空白字段惰瓜。./是指當(dāng)前目錄下的Human_ppi_symbol_2019_7.txt文件,當(dāng)然也可以寫絕對路徑汉矿。
a <- read.table("./Human_ppi_symbol_2019_7.txt",sep = "\t",comment.char = "#", stringsAsFactors = F,header = T, fill=TRUE)
寫入文件
寫入txt文件,不指定路徑則為當(dāng)前路徑下創(chuàng)建,學(xué)會舉一反三崎坊,參數(shù)有很多,這里并不能寫出所有洲拇,因此建議大家用?write.table查閱文檔奈揍,當(dāng)你不懂這個函數(shù)有什么參數(shù)的時候,就強(qiáng)烈建議赋续?+函數(shù)名來查閱文檔男翰,我就是這么學(xué)習(xí)的。
write.table(Gene, 'Gene.txt',sep = '\t', quote=F, row.names=F)
2.數(shù)據(jù)篩選相關(guān)
挑選特定的行和列
很多時候我們需要篩選數(shù)據(jù)纽乱,最常見的比如我希望得到某一列的固定邏輯的數(shù)據(jù)(我這里以大于900的數(shù)據(jù)為例)蛾绎,寫法是這樣的:
篩選大于900分的行,以a[x,y]為基礎(chǔ)鸦列,a$為catch住含有名為“combined_score”的列租冠,>900為篩選邏輯,前提要求數(shù)據(jù)類型為數(shù)字薯嗤,得到數(shù)據(jù)
這句話應(yīng)該這樣解析顽爹,a$combined_score > 900是邏輯值,返回True和Flase骆姐,你也可以理解為返回了對應(yīng)的行數(shù)镜粤,那么我們只要把這個行數(shù)再賦值給a[x,y]中的x捏题,那么得到的就是含有大于900的所有a 的x行。
a = a[a$combined_score > 900,]
取消科學(xué)計數(shù)法
取消科學(xué)計數(shù)法繁仁,有一些基因表達(dá)量太小或者p.Value太小涉馅,會用科學(xué)技術(shù)進(jìn)行計算,擔(dān)心后面篩除的時候會丟失
options(scipen = 400)
篩選同時滿足兩個要求的數(shù)據(jù)
篩選同時滿足P.Value < 0.05 和 logFC > 0.5的行黄虱,因為logFC有負(fù)數(shù)稚矿,可用abs絕對值化,這句話及上句都常用在做出來的DEG(差異表達(dá)基因)中
DEG_sort_p <- DEG_sort_p[DEG_sort_p$P.Value < 0.05 & abs(DEG_sort_p$logFC) > 0.5,]
常用語句
統(tǒng)計行列
dim(變量名)
行名稱更換
row.names(DEG_sort_p) = DEG_sort_p[,1]
把DEG_sort_p中的差異基因提出來捻浦,制作成一個data.frame
diff_Gene = data.frame(Gene=rownames(DEG_sort_p))
這里是蛋白互作的關(guān)鍵合并步驟晤揣,第①步,把基因中含有該蛋白的“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是去掉基因名字重復(fù)的基因昧识,只保留一個,保留有效基因
Gene=unique(data.frame(Gene=Gene[,2]))
判斷蛋白基因是否在全在基因里,返回真假值210122為True
table(Gene$Gene%in%DEG_sort_p$DEG)
unique是去掉基因名字重復(fù)的基因盗扒,只保留一個,保留有效基因
Gene=unique(data.frame(Gene=Gene[,2]))
簡單粗暴的合并行
Gene=rbind(Gene,diff_Gene)
把重復(fù)基因的行去掉
Gene=unique(Gene)
重命名行
colnames(data2)=c("Gene1","Gene2","score")
把data2的第一列變成character格式跪楞,這樣做的意義是有一些包要求輸入進(jìn)入格式是這種數(shù)據(jù)類型,
as.array()侣灶、as.data.frame()甸祭、as.matrix()、as.numeric()
data2$Gene1 <-as.character(data2[,1])
自創(chuàng)代碼:閱讀當(dāng)前目錄下所有文件褥影,如果里面包含GSE名稱的文件夾池户,則提取出,進(jìn)而閱讀里面包含exp_matrix的表達(dá)矩陣并進(jìn)行merge合并凡怎,要求這些GSE都是基于相同平臺校焦。意思是,當(dāng)前目錄下有N個GSE文件夾统倒,這些文件夾里面有各自已經(jīng)完成好的表達(dá)矩陣寨典,把這些表達(dá)矩陣都取出來并進(jìn)行merge合并成為一個大矩陣,適合做數(shù)據(jù)合并用房匆。
dir = list.dirs()
# 篩出路徑中包含GSE的文件夾
GSE_dir = dir[grep("GSE",dir)]
# 循環(huán)獲取目錄凝赛,要求根目錄下必須含有GSE的文件夾,GSE文件夾內(nèi)含有exp_matrix名稱的表達(dá)矩陣坛缕,要求表達(dá)矩陣是相同平臺,輸出變量名為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])
}
}
優(yōu)雅的增加一行捆昏,把行名稱加回去
GSEall_DEG_sort_p$Symbol = rownames(GSEall_DEG_sort_p)
求平均的函數(shù)赚楚,敲重要
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)
}