生信c语言,生信人的R使用

最新推荐文章于 2022-11-30 15:07:10 发布

weixin_39946534

最新推荐文章于 2022-11-30 15:07:10 发布

阅读量820

点赞数

文章标签：生信c语言

接下来介绍R语言：

【生信技能树】生信人应该这样学R语言

6bcaf21490ed

R语言

在你开始R之旅前，建议你看看下面这两个

1. 介绍R语言及Rstudio

了解R，Rstudio及R包;安装的包在packages中检查

.libPaths() #找安装路径

帮助文档，帮忙看路径

?substring

定位文件、设置文件位置

getwd()

setwd()

plot画板关闭dev.off()

2. R语言基础变量讲解

重点就是理解: 五种变量结构(class属性)

我也曾经写过一点这方面笔记：R语言学习笔记

grep()搜索函数

index1 = grep('RNA-Seq', a$Assay_Type)

index2 = grepl('RNA-Seq', a$Assay_Type)

b = a[index1,] # 下标

b = a[index2,] # 索引

3. 外部数据导入导出

Excel表格到R语言转换

去GEO数据库下载表达矩阵GSE17215

用excel打开后

6bcaf21490ed

GEO17215

此时用导入R会出现这种情况

> b=read.table('GSE17215_series_matrix.txt.gz',sep = '\t',comment.char = 'T',header = T)

Error in scan(file = file, what = what, sep = sep, quote = quote, dec = dec, :

line 29 did not have 2 elements

# 因为存在一些！开头的文件存在，此时经过下面这种处理：

# 带感叹号的不读取

> b=read.table('GSE17215_series_matrix.txt.gz',sep = '\t',comment.char = '!',header = T)

> View(b)

保存，进一步处理GSE17215_series_matrix.csv

# 保存为.CSV格式

write.csv(b,'GSE17215_series_matrix.csv')

# 第一列设为行名，再去掉第一行

rownames(b)=b[,1]

b=b[,-1]

经过上面处理后如下；

6bcaf21490ed

image.png

画个热图看看看

b = log2(b)

pheatmap::pheatmap(b[1:10,])

6bcaf21490ed

热图

保存b推荐采用下面这种方式进行

save(b,file = 'b_input.Rdata')

load(file = 'b_input.Rdata')

建议把excel转成csv来读取

4. 中级变量操作

所有函数有参数，很多是互通的， eg sort, max, min, fivenum(四分位值)

函数

# 最大值

> sort(b$GSM431121,decreasing = T)[1]

[1] 15.28882

> max(b$GSM431121)

[1] 15.28882

# 最小值，五分位数

> min(b$GSM431121)

[1] 3.859587

> fivenum(b$GSM431121)

[1] 3.859587 4.710004 5.952603 8.691293 15.288819

向量化操作

> table(b$GSM431121<5)

FALSE TRUE

14398 7879

> d=b[b$GSM431121<5,]

# 你可以看看你的b，是为有7879行

看b的每行的平均值

> mean(b[1,])

[1] NA

Warning message:

In mean.default(b[1, ]) : 参数不是数值也不是逻辑值：回覆NA

> mean(as.numeric(b[1,]))

[1] 10.60797

> as.numeric(b[1,])

[1] 8.911691 9.221081 11.410364 11.325483 11.418782 11.360438

> mean(as.numeric(b[1,]))

[1] 10.60797

采用rowMean函数，head前6行

> head(rowMeans(b))

1007_s_at 1053_at 117_at 121_at 1255_g_at 1294_at

10.607973 7.925899 5.193894 7.168633 4.275652 5.686036

# 采用for循环看看

for (i in 1:nrow(b)) {

print(mean(as.numeric(b[i,])))

}

for (i in 1:6) {

print(mean(as.numeric(b[i,])))

}

> for (i in 1:6) {

+ print(mean(as.numeric(b[i,])))

+ }

[1] 10.60797

[1] 7.925899

[1] 5.193894

[1] 7.168633

[1] 4.275652

[1] 5.686036

使用apply循环也可以实现

x=mean(as.numeric(b[1,]))

apply(b,1,function(x){

mean(x)

})

查找最大值

# 最大值的查找

for (i in 1:nrow(b)) {

print(max(as.numeric(b[i,])))

}

apply(b,1,function(y){

max(y)

})

apply(b,1,max)

# 自写函数rowMax查找

rowMax=function(z){

apply(z, 1, max)

}

rowMax(b)

计算每行的方差，取最大的50个，画热图

apply(b, 1, sd)

sort(apply(b, 1, sd),decreasing=T)[1:50]

cg=names(sort(apply(b, 1, sd),decreasing=T)[1:50])

pheatmap::pheatmap(b[cg,])

6bcaf21490ed

image.png

随机选取50ge，sample函数

sample(1:nrow(b),50)

pheatmap::pheatmap(b[sample(1:nrow(b),50),])

6bcaf21490ed

随机50个

+ abs, sqrt :戒对治，平方根

+ Log, log10, log2, exp: 对数与指数函数

+ Sin, cos, tan, acos, atan, atan2: 三角函数

+ sinh, cosh, tanh, asinha, acosh, aranh：双曲函数

+ 集合运算，reshape, merge总结

思考一下excel表格里面有变量类型吗

a = read.table('XXtable.txt', head = T

sep = '\t')

b = read.table('GSE17215_series_matrix.txt.gz',

comment,char = '!', head = T,

sep = '\t')

write.csv(b, 'GSE17215_series_matrix.csv')

write.table(b,'tmp.csv', sep = ',')

##把行名去掉

d = read.csv('GSE17215_series_matrix.csv')

# readline 读入之后拆分

5. 热图

随机产生a1，并产生热图

a1=rnorm(100) #随机产生100个服从正态分布的数

?rnorm

dim(a1)=c(5,20) # 添加维度属性，矩阵matrix

pheatmap(a1)

6bcaf21490ed

产生a2+热图

a2=rnorm(100)+2

dim(a2)=c(5,20)

pheatmap(a2)

6bcaf21490ed

将a1,a2横向拼接并画图

pheatmap(cbind(a1,a2),cluster_cols = F)

a3=cbind(a1,a2) #横向拼接

a4=rbind(a1,a2) #纵向拼接

6bcaf21490ed

a1,a2

6. 选取差异明显的基因的表达量矩阵绘制热图

rm(list = ls()) #魔幻操作，一键清空

library(pheatmap)

a1 = rnorm(100)

dim(a1) = c(5,20)

pheatmap(a1)

a2 = rnow(100)+2

dim(a2) = c(5,20)

library(pheatmap)

pheatmap(a1, cluster_rows = F, cluster_cols = F)

pheatmap(cbind(a1,a2))

pheatmap(cbind(a1,a2), show_rownames = F, show_colnames = F)

拉平极差值

7. ID转换

b=as.data.frame(a3)

paste('a1',1:20,sep = '_')

paste('a2',1:20,sep = '_')

names(b)=c(paste('a1',1:20,sep = '_'),paste('a2',1:20,sep = '_'))

pheatmap(b,cluster_cols = F)

ENSG基因ID处理

#strsplit('','[.]') #根据点号分割

> strsplit('ENSG0000000003.13','[.]')

[[1]]

[1] "ENSG0000000003" "13"

> strsplit('ENSG0000000003.13','[.]')[[1]]

[1] "ENSG0000000003" "13"

> strsplit('ENSG0000000003.13','[.]')[[1]][1]

[1] "ENSG0000000003"

#ID转换包

library(stringr)

a$ensemble_id=str_split(a$v1,'[1]',simplify = T)[,1]

# duplicated() #去重

一些包

org.Hs.eg.db #在包里有基因注释关系

8. 任意基因任意癌症表达量分组的生存分析

也可以下载原始数据自己分析

#读取csv文件

rm(list = ls()) #清除所有变量

options(stringsAsFactors = F)

a=read.table('LGG_93663_50_50.csv',header = T,sep =',',fill = T)

# 后续画图

colnames(a)

dat=a

# 图ggbetweenstats

library(ggstatsplot)

ggbetweenstats(data = dat,x=Group,y=Expression)

# 图ggsurvplot

library(ggplot2)

library(survival)

library(survminer)

table(dat$Status)

dat$Status

sfit

summary(sfit)

ggsurvplot(sfit, conf.int = F, pval = T)

#ggsave('survival_ARHGAP18_in_LGG.png')

ggsurvplot(sfit,palette = c("#E7B800","#2E9FDF"),

risk.table = TRUE,pval = TRUE,

conf.int = TRUE,xlab="Time in months",

ggtheme = theme_light(),

ncensor.plot=TRUE)

#ggsave('survival_ARHGAP18_in_LGG.png')

6bcaf21490ed

ggbetweenstats

6bcaf21490ed

ggsurvplot

9. 任意基因任意癌症表达量和临床性状关联

某个基因在某个癌症的表达量，关联临床信息

rm(list = ls()) #清除所有变量

options(stringsAsFactors = F)

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)

6bcaf21490ed

image.png

10. 表达矩阵样本的相关性

看两个变量的相关性

> cor(1:10,1:10)

[1] 1

> a = rnorm(10)

> b = rnorm(10)

> cor(a,b)

[1] -0.1555608

> a = rnorm(10)

> b = 10*a+rnorm(10)

> cor(a,b)

[1] 0.9971822

学习airway这个数据包

rm(list = ls()) #清除所有变量

options(stringsAsFactors = F)

library(airway)

data("airway") #加载数据

exprSet=assay(airway)

colnames(exprSet)

group_list=colData(airway)[,3]

exprSet=exprSet[apply(exprSet,1,function(x) sum(x>1)>5),]

exprSet=log(edgeR::cpm(exprSet)+1)

exprSet=exprSet[names(sort(apply(exprSet, 1,mad),decreasing = T)[1:500])]

学习上面得到的exprSet

6bcaf21490ed

image.png

> dim(exprSet)

[1] 64102 8

> cor(exprSet[,1],exprSet[,2]) #查看相关性

[1] 0.9632268

可知第一列和第二列的相关性较好：

生信c语言,生信人的R使用

“相关推荐”对你有帮助么？