Day15-R学习-生信人应该这样学R语言-笔记(转载)

01. 介绍r语言及rstudio编辑器

安装r及rstudio
打开rstudio编辑器:4块
1.新建脚本,markdown, 写source代码,点任何地方会进行下一步。
2.命令在console运行,初学一步一步运行代码,上下左右键看历史命令记录。
3.history为历史记录,可以to source 或 to console进行重新运行,environment环境变量,赋值的变量会出现在这个位置。
4.files查看有哪些文件或浏览文件
plot显示画出的图,dev.off()关闭画板
packages显示安装的所有包,.libPaths()显示包安装路径
Help 帮助文档查看运行函数example,体会如何使用该函数

定位当前文件位置:getwd()

02. R语言基础变量讲解

变量类型
**1.向量vector **
a=c(1,2,3)
b=c(1,'2',3)
class() 查看变量类型(a为numeric,b为character)

可以用函数创建向量,或直接使用内置变量,左边变量,右边值。
数字不加引号,字符串加引号(*单引号或双引号一般可通用,当值内本身含有单引号时,需用双引号。例如:a='Hello World!'a="Hello World!"输出内容相同.b="This's a nice example."b='This's a nice example.'不同。)

2.矩阵matix
向量加维度为矩阵 dim()加维度,矩阵中任何一个元素类型发生变化,其他元素的类型也会变化。
取元素的方式:1.通过下标来取,#逗号左边是行右边是列 a[,]; 2.逻辑符取
判断元素类型: class(), str()

3.数组array

4.data.frame()数据框
转变其中一个元素的类型,其他元素的类型不变

is.系列函数,as.系列函数
是什么函数和转变为什么函数(用Table键补全时可看到很多函数)

提问方式
新建文件夹,放入需要提问/报错的“.Rproj”以及代码".R"标注报错的位置,发压缩包。
如果其中有变量,报错之前一步写
save(filename,file='input.Rdata')
load(file='input.Rdata')
然后注释掉load之前的代码(加上#)

5.列表list
用'$'从列表中取出的是元素,数据框中取出的是一列
length()有多少元素,lapply()对每个元素操作,返回list,用unlist()让返回值为character,用as.numeric()将其转变为数值型。

取下标

  1. 用下标取索引
    b[,]逗号左行右列。若没有该元素,返回NA;若有,返回对应元素
  2. 用判断方式取下标
    b[,c(F,F,T,F,F,)]取b第三列的元素

read.table()函数读取table
grep()#搜索函数,例grep('RNA-Seq',a$Assay_Type)搜索在a的Assay_Type一列中,搜索含有RNA-Seq的下标。
grepl()#取匹配到的所有的行,返回TRUE/FALSE,例grepl('RNA-Seq',a$Assay_Type)

取list
list[]取元素,不能取元素里的内容,拿出的元素可能还是一个list,如果用[[*]][*]取出的是元素里面的内容

03. 外部数据导入导出

读入
1.直接用import dataset
2.代码读取
read.table('filename',sep='\t',header = T)#参数与参数之间用逗号分开,sep指定制表符,header加上表头。

write.csv()#读出为csv文件

取b第一列作为b的行名,然后去掉第一列
row.names(b)=b[,1]
b=[,-1]

保存为R读取的格式,避免读入读出时格式的更改
save(b,file = 'b_input.Rdata')
load(file = 'b_input.Rdata')

04. 中级变量操作

读入数据
read.table(filename,header = T,sep = '\t')
read.table(tablename,comment.char = "!",header = T,sep = '\t')#手动指定注释符,以!开头的不读入comment.char = "!"

读出csv文件,并且去掉行名两种方式
write.csv(b,'tmp.csv',row.names=F)
write.table(b,'tmp.csv',sep = ',')

英文单词为函数,函数有参数,互通
sort()#分类
max()#最大值
min()#最小值
fivenum()#包括最小值, 25%分位数, 中位数, 75%分位数, 最大值

table()#查看多少个元素
boxplot(y~x)#分组取域值,可用boxplot查看峰位数
先分组,分别赋值,"="是赋值,"=="是判断,
rna=a[a$Assay_type=='Rna-Seq',]#a中Assay_type这一列中Rna-Seq对应的行
wxs=a[a$Assay_type=='WXS',]#a中Assay_type这一列中WXS对应的行

表达矩阵,不同的表达量
str(mean(b[1,]))#取b的第一行平均值报错,查看其类型
mean(as.numeric(b[1,]))#转换为数值型取均值
head(rowMeans(b))#对所有的行取均值,查看前10个

取b中行的均值
rowMeans(b)#直接用函数
for (i in 1:nrow(b)){ print(mean(as.numeric(b[i,])))}#用循环
apply(b,1,function(x){mean(x)})#apply函数
apply(b,1,mean)#与上一个相同

有函数可直接用函数,没有函数可以自定义后进行操作
rowMax=function(x){apply(x,1,max)}#定义一个函数
rowMax(b)#用定义的函数进行操作

向量取元素直接写下标,data.frame()list()取元素[,]左边取行右边取列

top50的方差对应的基因名做热图
cg=names(sort(apply(b,1,sd),decreasing = T)[1:50])
library(pheatmap)
pheatmap(b[1:50,])
pheatmap(sample(1:nrow(b),50))#随机选取50个方差

05. 热图

不同大小的数据达标不同颜色的深浅
rnorm(n, mean = 0, sd = 1)#随机正态分布

画个热图

a1=rnorm(100)
#dim加维度
dim(a1)=c(5,20)
a2=rnorm(100)+2
dim(a2)=c(5,20)
#as.data.frame将向量变为矩阵
b=as.data.frame(cbind(a1,a2))
#paste函数给a1和a2取名并以此作为b的名字看热图的变化,
names(b)=c(paste('a1',1:20,sep = '_'),paste('a2',1:20,sep = '_'))
#加载pheatmap包
library(pheatmap)
#加注释信息
tmp=data.frame(group=c(rep('a1',20),rep('a2',20)))
rownames(tmp)=colnames(b)
pheatmap(b,cluster_cols = F) #cluster_cols = F 不进行排序

pheatmap学习,查看example
?pheatmap

06. 选取差异名明显的基因的表达矩阵绘热图

#均一化(normalization)调整差异过大的数,转置-scale-转置。
n=t(scale(t(dat[cg,])))
n[n>2]=2#大于二的值均等于二
n[n<-2]=-2#小于负二的值均等于负二

?scale()#查看scale函数的用法

若代码过长,可以定义一个函数将代码包裹起来,后面如果有需要可以直接用不重复写代码
d_h <- function(dat,group_list){*}
d_h(dat,group_list)

07. id转换

#矩阵a的V1列元素含有以'.'分割的内容,想要去除‘.’后的内容
library(stringr)#载入stringr
str_split(a$V1,'[.]',simplify = T)[,1]#用str_split()函数,simplify=T返回值为字符型矩阵,F返回为字符型向量,[,1]取第一列

library(org.Hs.eg.db)#加载org.Hs.eg.db包
toTable()#org.Hs.eg.db包规定读入数据用toTable

#两个数据框做关联,两个数据框共有相同的列/行名,保留没有关联上的元素 
merge(a,b,by='*',all.x=T)
#出现的频率,对大于1的做统计
table(table(*)>1)
table(*)[table(*)>1]
#dataframe(d)的行去重
d=d[!duplicated(),]
#a所在的顺序放在b这边来
match(a,b)

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

一个基因在TCGA数据库中的各个癌症的生存分析
网页中输入感兴趣的基因,分组(高低表达量),下载数据,
Rstudio中读入数据,ggstatsplot包中的代码进行分析

a=read.table('*.csv',header = T,sep = ',',fill = T)
ggbetweenstats(data=a,x='Group',y='Expression')

添加Group和Expression最好复制加'',避免打错和不识别。
学会基础变量及函数操作,然后对包的说明书进行学习,数据应用

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

网页中输入感兴趣的基因,下载数据,导入Rstudio
a=read.table('*.txt',header = T,sep = '\t',fill = T)
列名赋予简单的值,dat与之前写过的代码中一致,不需改动代码可直接运行(个人习惯)

colnames(a)=c('id','stage','gene','mut')
dat=a

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

cor() 关联函数

#先安装Bioconductor的代码,需要R版本3.5.0
if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
BiocManager::install()

Bioconductor三个包:airway(数据包),Annotatiobdbi(注释包)和GenomicFeatures(功能函数包)

#加载数据,获取表达矩阵,查看表达,用dim()维度函数获得矩阵行列情况
data(airway)
exprSet=assay(airway)
colnames(exprSet)
dim(exprSet)
#样本与样本之间的相关性,
cor(exprSet[,1],exprSet[,2])

相关性很高警惕两种情况
1.同一个样本技术重复
2.大多数值为低表达值或零,掩盖真实相关性
cor(exprSet)#整个样本做相关性

#筛选数据,例如一个基因表达量大于一样本量少于5,用dim()维度函数获得矩阵行列情况,可以与之前的dim()做比较
x= exprSet[1,]
table(x>1)
true=1,false=0
sum(x>1)>5
exprSet= exprSet[apply(exprSet,1,sum(x>1)>5),1]
#继续筛选,edgeR::cpm()函数去除文库大小差异,mad的前50个变化最大的基因
#最后cor做分析,画热图
M=cor(log2(exprSet+1))

#画heatmap 加注释信息,先分组(group_list)并转为一个数据框,加注释信息(多少列不重要,行名等于矩阵的列名很重要)tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
library(pheatmap)
pheatmap(cor(exprSet),annotation_col = tmp)

11.芯片表达矩阵下游分析

构建表达矩阵,分组信息

#单个做差异分析
t.test(exprSet[1,]~ group_list)
boxplot(exprSet[1,]~ group_list)
#或用limma()包,构造一个design矩阵,矩阵中显示分组信息属于(值为1)不属于(值为0)

12.RNA-seq表达矩阵差异分析

airway()
assay()
表达矩阵,分组准备,用DEseq2()
统计学相关的知识:statquest

原文

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 204,684评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,143评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,214评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,788评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,796评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,665评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,027评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,679评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 41,346评论 1 299
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,664评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,766评论 1 331
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,412评论 4 321
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,015评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,974评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,203评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,073评论 2 350
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,501评论 2 343