1.导入gz压缩文件
library(R.utils)
gzfile ='D:/data by multiviolin/GSE66507_human_blastocyst_rnaseq_counts.txt.gz'
gunzip(gzfile, remove = F)
gf = gzfile(gzfile,'rt')
data = read.table(gf, header = F, sep = '\t', quote = '', comment.char = '#')
2.数值型转换
tdf2=as.data.frame(lapply(tdf,as.numeric))
3.保留不为0的行或列
X[which(rowSums(X) > 0),]
4、循环导出gz文件
file=list.files()
n=length(file)
for (i in file) {gzfile =i
gunzip(gzfile, remove = F)}
5.循环处理txt文件
5.1 取文件名的名字
name<-substr(basename(i),start = 22,stop = 25)
5.2 给数据框增加一列内容
b<-transform(b,group=name)
5.3 取文件夹下的所有文件
file=list.files('all late')
合集
setwd('C:/Users/l1047/Desktop/mice')
file=list.files('all late')
n=length(file)
setwd('all late')
for (i in file) {
name<-substr(basename(i),start = 22,stop = 25)
b<-read.table(file=i,header = T)
b<-b[,-c(2,4,5,6)]
b<-transform(b,group=name)
write.table(b,i,row.names = FALSE,sep = "\t")
}
###改变列名
用colnames<- c (" "," ")
###导入txt文件时,由于第一个字符是#而导致第一行丢失
在read.table中加入参数comment.char=""即可
导入txt文件时报错显示‘line 8 did not have 9 element’
在read.table中加入参数 sep='\t' 即可
##用for循环合并多个txt文件
setwd('C:/Users/l1047/Desktop/mice')
allfile<-list.files('all late')
n=length(allfile)
newdata<-c()
setwd('all late')
for (i in allfile){
rb1<-read.table(i)
name<-substr(basename(i),start = 22,stop = 25)
colnames(rb1)<-rb1[1,]
rb1<-rb1[-1,]
colnames(rb1)<-c("gene",name,"group")
rb1<-rb1[,-3] #对第一个文件处理完毕
newdata<-cbind(newdata,rb1)
}
##删掉某列中重复的行
newdata<-newdata[!duplicated(newdata$gene), ] #删掉行
###删掉重复的列
newdata<-select(newdata,-starts_with("gene")) #删掉以“gene”开头的列
#保留不为0的列,即去掉全是0的列
exp2<-exp[,which(colSums(exp)>0)]
给数据框里的列中的值改名
newsp6[,3][newsp6[,3]=='epiblast']='EPI'
删掉数据框里某一列的含有某个值的所有行
sp6<-newsp6[-which(newsp6$lineage=='not applicable'),]
给数据框增加一列
b<-transform(b,group=name)
提取数据框某列每一行的前两个字符
name<-substr(human2$lineage,1,2)
#substr(data,start, stop)
更改数据框列的顺序
library(dplyr)
human4<-select(human3,SP6,lineage,group)
#select(data,x,y,z,everything()),everything是其他都不变
ggplot画图
library(ggplot2)
nocell$Sp6<-as.numeric(nocell$Sp6)#先把数据变成numeric,否则报错
ggplot(data=nocell,aes(x=lineage,y=Sp6,fill=group))+
geom_violin(scale = "width")+geom_jitter( position = position_jitterdodge(0.35))
更换字符串中的字符
chartr(old,new, X)#X是文件名
##ggplot颜色
##配色
install.packages('RColorBrewer')
library(RColorBrewer)
help("RColorBrewer")
brewer.pal.info #查看色板信息
#查看所有色板。
display.brewer.all()
#查看单色渐变色板。
display.brewer.all(type='seq')
#查看双色渐变色板。
display.brewer.all(type='div')
#查看离散颜色色板。
display.brewer.all(type='qual')
#查看BrBG色板中的9个颜色。
display.brewer.pal(9,'BrBG')
#查看BrBG色板中的9个颜色的色值。
brewer.pal(9,'BrBG')
display.brewer.pal(9,'RdYlBu')
brewer.pal(9,'RdYlBu')
pal<-brewer.pal(9,'RdGy')
mycolor<-colorRampPalette(pal)
image(x=1:9,y=1,z=as.matrix(1:9),col=mycolor(9))
##颜色处理
library(scales)
show_col(hue_pal()(n)) #ggplot 默认颜色
show_col(hue_pal(h = c(0, 180) + 30, c = 100, l = 65)(3))
p0<-p2+scale_fill_hue(h = c(0, 180) + 30, c = 100, l = 65)#语法
##ggplot画火山图并标记差异基因(也可用ggscatter)
setwd('C:/Users/l1047/Desktop/1124')
data<-read.csv("C5_D0-vs-C5_D6.edgeR.all(1).xls",sep = '\t',header = T)
# 设置p_value和logFC的阈值
cut_off_pvalue = 0.0000001 #统计显著性
cut_off_logFC = 1 #差异倍数值
data$change <-ifelse(data$PValue< cut_off_pvalue & abs(data$log2.fc.) >= cut_off_logFC,
ifelse(data$log2.fc.> cut_off_logFC ,'Up','Down'),'Stable')
data2<-data[,c('id','PValue','log2.fc.','change')]
library(ggplot2)
p<- ggplot(data2, aes(x = log2.fc., y = -log10(PValue), colour=change)) +
geom_point(alpha=0.4, size=3.5) +
scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757"))+
# 辅助线
geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
geom_hline(yintercept = -log10(cut_off_pvalue),lty=4,col="black",lwd=0.8) +
# 坐标轴
labs(x="log2(fold change)",
y="-log10 (p-value)")+
theme_bw()+
# 图例
theme(plot.title = element_text(hjust = 0.5),
legend.position="right",
legend.title = element_blank())
library(ggrepel)###加上标签
p1<-p+ geom_label_repel(data = data2, aes(x = log2.fc.,
y = -log10(data2$PValue),
label = id),
size = 3, box.padding = unit(0.5, "lines"),
point.padding = unit(0.8, "lines"),
segment.color = "black",
show.legend = FALSE)
##基因名互相转换
library(org.Hs.eg.db)
library(hgu95av2.db)
library(dplyr)
library(clusterProfiler)
symple<-read.table('symple.txt',header = F,sep = '\t')
symple<-as.character(symple$V1)
library(org.Hs.eg.db)
eg <- bitr(symple,fromType = 'SYMBOL',
toType = 'REFSEQ',
OrgDb='org.Hs.eg.db')
#字符替换
1.gsub(原,后,文件)#可用^放在开头,表示只替换开头的而不换其他的
2.tidyverse的str_remove(文件,parttern=“”)
3.sub,同1
###循环画图
dir.create('Vlnplot for placenta') ###创建新文件夹
setwd('Vlnplot for placenta') ##改变工作路径
for(i in rownames(mydata)){ #读取样品名
result_name = paste(i,".pdf",sep="") #根据样品名生成PDF文件
pdf(result_name,10,7)
p<-VlnPlot(mydata,features = i,group.by = 'cell_type')
print(p)
dev.off()
}
###行名有重复值报错
row.names(newdata)<-make.names(newdata[,1],TRUE)#行名有重复值时的办法
##给R分配多个核同时运算
##让R多芯片并行运算
#加载parallel包
library(parallel)
#detectCores函数可以告诉你你的CPU可使用的核数
clnum<-detectCores()
#设置参与并行的CPU核数目,这里我们使用了所有的CPU核,也就是我们刚才得到的clnum,具体到这个案例,clnum=4
cl <- makeCluster(getOption("cl.cores", 32))
#关闭并行计算
#stopCluster(cl)
###数据框转换成数值型
data5 <- as.data.frame(lapply(data4[,1:64],as.numeric))
###将右表匹配进左表
library(dplyr)
data1 <- left_join(list,data,by=c('V1'='gene'))##将右表内容匹配进左表
all(data1$V1==data1$gene) #查看两列是否匹配
data2<- na.omit(data1)##删除缺失行
###长宽数据转换(reshape2包)
宽转长,用melt函数
- id.vars:选择用来做主键的列,对于表2来说就是
date
列,当然也可以有多列 - measure.vars:观测值的列,默认是非
id.vars
那些列,表2中就是除第一列的其他列 - variable.name:
measure.vars
转换成一列后的列名,这里就是“stock” - value.name:观测值转换成一列后的列名,这里是“value”
- na.rm:是否删除
NA
数据