【R语言学习】生信人应该这样学R语言--学习笔记（2）

保存及加载数据:

save(a,file="GSE17215.Rdata")
load("GSE17215.Rdata")

清除内存变量、释放内存：

rm(list = ls())  #清除内存变量
gc()  #释放内存

画一个热图：

#参考https://www.jianshu.com/p/4e3a94bd155a 第五部分
a<-read.table("GSE17215_series_matrix.txt",comment.char = "!",header = T,sep = "\t")
write.csv(a,"a1.csv")

rownames(a)<-a[,1]
a<-a[-1]

a<-log2(a)
library(pheatmap)
pheatmap(a)

OncoLnc网站 获取某基因在TCGA数据库中的表达情况

绘制生存曲线：

#从网页(例如OncoLnc)获取数据，绘制生存曲线
a<-read.csv("LIHC_93663_50_50.csv",header = T,sep = ",",fill = T)
colnames(a)
dat<-a
library(ggstatsplot)
ggbetweenstats(data=dat,x=Group,y=Expression)
ggbetweenstats(data=dat,x=Status,y=Expression)
library(ggplot2)
library(survival)
library(survminer)
table(dat$Status)
dat$Status<-ifelse(dat$Status=="Dead",1,0)
sfit<-survfit(Surv(Days,Status)~Group,data=dat)
sfit
summary(sfit)
ggsurvplot(sfit,conf.int=F,pval=T)
ggsave("survival_ARHGAP18_in_LGG.png")

ggsurvplot(sfit,palette = c("#E7B800", "#2E9FDF"),
           risk.table=T,pval=T,conf.int=T,
           xlab="Time in months",ggtheme = theme_light(),ncensor.plot=T)
ggsave("survival_ARHGAP18_in_LGG_2.png")

cBioPortal 网站获取临床信息

表达量与临床分期画图：

a<-read.table("plot.txt",header = T,sep="\t",fill = T)
colnames(a)<-c("id","stage","gene","mut")
dat<-a
library(ggstatsplot)
ggbetweenstats(data = dat,x = stage,y = gene)
ggsave("ov_ARHGAP18.png")

#改变统计方法
res.aov <- aov(gene~stage,data=dat)
summary(res.aov)
TukeyHSD(res.aov)

相关性cor(a,b) --->a,b之间的相关性,cor(exprSet)可以查看exprSet矩阵所有列之间的相关性；

dim()`可以看某数据的维度。

表达矩阵样本的相关性：

rm(list = ls())
options(stringsAsFactors = F)
library(airway)
data("airway")
exprSet<-assay(airway)  #h获得表达矩阵

colnames(exprSet)
cor(exprSet)
pheatmap::pheatmap(col(exprSet))

group_list<-colData(airway)[,3]
tmp<-data.frame(g<-group_list)
rownames(tmp)<-colnames(exprSet)
pheatmap::pheatmap(col(exprSet),annotation_col = tmp)
dim(exprSet)
exprSet<-exprSet[apply(exprSet,1,function(x) sum(x>1) >5),]
dim(exprSet)

exprSet<-log(edgeR::cpm(exprSet)+1)   #去除文库大小差异，并得到新矩阵
exprSet<-exprSet[names(sort(apply(exprSet,1, mad),decreasing = T)[1:500]),]
dim(exprSet)
M<-cor(log2(exprSet+1))
tmp<-data.frame(g<-group_list)
rownames(tmp)<-colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp,filename = "cor.png")