终于整理好了
P1 介绍R语言及Rstudio编辑器
00:45 下载R语言软件 for win
00:53 下载Rstudio Rstudio
-
01:47 Rstudio界面
-
03:45 安装R包
复制代码到左上Source区块。每点击一行，Console区块出现一行命令。Or, 一次性把所有命令粘贴到Console区块，但无法查看/控制进程。键盘↑↓键可以按顺序调出已经输入的命令。
-
5:20
plot()画图在未开窗口，e.g.
png(’temp.png‘)，的情况下，plot()绘制出的图像在右下“Plots”窗口显示。 06:18
.libPaths()查看包的安装路径-
08:51
?function查询函数帮助文档e.g.
?subString 10:25 查看当前路径
getwd()，设置工作路径setwd()
P2 R语言基础变量讲解
01:55 创建向量，
c(), e.g.a=c(1,2,3)01:58 查看向量类型，
class()-
02:54 创建向量
a=1:10a=LETTERS（R语言内置常量Built-in Constants） -
03:20 取内置常量的一部分
LETTERS[1:10] 04:00 字符串需要加引号，一般情况下单双引号没有区别
-
04:56 向量→矩阵
> a=1:10 > dim(a)=c(2,5) -
06：26
pheatmap::pheatmap(a) -
06:36 赋值
> a[1,2]='5'ps.
'5'为character导致
此时则不能生成热图
-
07:09 判断工作对象 (object) 的性质
> class(a) > str(a) -
09:11 判断 obejct 是否是
matrix或其他> is.matrix(a) > is.data.frame(a)
-
11:12 将 object 转换为
data.frame格式> as.data.frame(a) > b=as.data.frame(a) > b[1,2]='5' > pheatmap::pheatmap(b)Error in cut.default(x, breaks = breaks, include.lowest = T) :
'x' must be numeric要保证执行操作前输入数据的性质符合函数要求。
-
17:04
$可查看数据框中的列> d=options() -
18:39 查看
d$repos> e=d$repos > class(e) > mode(e) > typeof(e) -
19:16 查看 d 中每一个元素的长度
> lapply(d,length)d 中一共有87个元素，这样查看很不直观/方便。
-
19:56 将 d 变为非 list 数据，查看每个元素的长度
> unlist(lapply(d,length)) -
20:19 只查看长度，不显示名称
> as.numeric(unlist(lapply(d,length))) -
24:15 学习在数据框中抓取元素（取子集）的数据准备
> a=read.table('YourPath/SraRunTable.txt',stringsAsFactors=F,header=T,sep='\t') -
24:16 查看 a 的性质
-
24:58 任意取几行数据成为一个
temp, 以 ‘Assy_Type’ 为 ‘RNA-Seq’ 为搜索项不过这份数据里都是 RNA-Seq, 不是非常能体现意义
> temp=a[c(1,2,7),] -
25:25 上面的方法需要用肉眼判断每个元素的坐标，而采取 TRUE or FALSE 的方法更高效
> grep('RNA-Seq',a$Assay_Type)取到了 'Assay_Type' 为 ‘RNA-Seq’ 的所有元素的下标
> grepl('RNA-Seq',a$Assay_Type)判断每个元素的 'Assay_Type' 是否为 ‘RNA-Seq’
-
28:00 两种方法核对
> index1=grep('RNA-Seq',a$Assay_Type) > index2=grepl('RNA-Seq',a$Assay_Type) > b=a[index1,] > b=a[index2,] -
29:11 两种取子集方法在 list 中的应用
接着之前的
d=options()> as.numeric(unlist(lapply(d,length))) > 2> table(as.numeric(unlist(lapply(d,length))) > 2)> index3=as.numeric(unlist(lapply(d,length))) > 2 > d=d[index3]list 是层层叠加的，所以在
d=d[index3]中没有代表“行、列”的逗号。 -
30:47 验证 list 的层层叠加
> d[4] > class(d[4])
P3 外部数据导入导出
-
02:12 直接导入
-
02:37
read.table()读取表格文件参数中：
sep='\t'，以制表符为分隔符header=T，设第一行为表头（第一行比其他行少一列时）> a=read.table('SraRunTable.txt',sep='\t',header=T)image
comment.char=''，' ' 后的内容不读取
> b=read.table('GSE17215_series_matrix.txt.gz',sep = '\t',comment.char='!',header = T)
-
08:53
write.csv()创建/另存为 csv 文件> write.csv(b,'GSE17215_series_matrix.csv') -
11:08
rownames()=[,1]，第一列设为行名
[,-1]，去掉行名> rownames(b)=b[,1]> b=b[,-1] -
12:05 随机取数据绘制热图
> pheatmap::pheatmap(b[1:10,])以log2值呈现
> b=log2(b) > pheatmap::pheatmap(b[1:10,]) -
12:44 储存为 Rdata 格式
> save(b,file = 'b_input.Rdata')
P4 中级变量操作
-
02:36 去行名
> write.csv(b,'GSE17215_series_matrix.csv',row.names = F) -
05:08 筛选
> sort(a$MBases) > sort(a$MBases)[1] > sort(a$MBases,decreasing = T)[1][1], 为第一个/最小值
decreasing = T，从大到小 -
06:44 最大值、最小值、五数概括法（最小值；第1四分位数(Q1)；中位数(Q2)；第3四分位数(Q3)；最大值）
> max(a$MBases) > min(a$MBases) > fivenum(a$MBases) -
08:17 筛选
> a$MBases > 5 > a$MBases < 5查看筛选结果的数量
> table(a$MBases < 5) -
09:51 根据
Assay_Type分组> boxplot(a$MBases~a$Assay_Type)这份数据里只有 'RNA-Seq' 一种 Assay_type
-
11:56 分别取
fivenum()（有不同 Assay_type 的情况下）> wes=a[a$Assay_type=='WXS'] > rna=a[a$Assay_type=='RNA-Seq'] > fivenum(wes$MBases) > fivenum(rna$MBases) -
16:19 取行平均值
> as.numeric(b[1,]) > rowMeans(b) -
16:52 引入循环
for循环，得到所有行的平均值> for (i in 1:nrow(b)) { print(mean(as.numeric(b[i,]))) }apply循环，得到所有行的平均值> x=as.numeric(b[1,]) > apply(b,1,function(x){mean(x)})for循环，得到所有行的最大值> for (i in 1:nrow(b)) { print(max(as.numeric(b[i,]))) }apply循环，得到所有行的最大值> apply(b,1,function(x){max(x)}) -
25:05 任意设置一个函数
> testf <- function(b){ for(i in 1:nrow(b)){ x=as.numeric(b[i,]) y=x[1]+x[2]-x[3]+x[4]-x[5]+x[6] print(y) } } > testf(b) -
26:07 求方差
求每一行的方差
> apply(b,1,sd)方差最大的前50个探针
> sort(apply(b,1,sd),decreasing = T)[1:50]得到这些探针的名字，绘制热图
> cg=names(sort(apply(b,1,sd),decreasing = T)[1:50])> pheatmap::pheatmap(b[cg,])
P5 热图
-
01:18 直接用矩阵画图
> a1=rnorm(100) > dim(a1)=c(5,20) > library(pheatmap) > pheatmap(a1)> a2=rnorm(100)+2 > dim(a2)=c(5,20) > pheatmap(a2) -
02:25 将 a1, a2 两个矩阵绘制在同一副热图
> pheatmap(cbind(a1,a2))> pheatmap(cbind(a1,a2),cluster_cols = F)cluster_cols = F，不排序 -
03:50 转换为
data.frame格式，添加列名> b=cbind(a1,a2) > b=as.data.frame(b) > c(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) -
12:57 加上
annotation_col参数> tmp=data.frame(group=c(rep('a1',20),rep('a2',20)))rep()，重复
给 tmp 加上列名，画图
> row.names(tmp)=colnames(b)
> pheatmap(b,annotation_col = tmp)
P6 选取差异明显的基因的表达量矩阵绘制热图
-
01:00 选1000个基因绘制热图
> cg=names(tail(sort(apply(dat, 1, sd)),1000)) > library(pheatmap) > pheatmap(dat[cg,],show_colnames = F,show_rownames = F) -
01:56 归一化
> n=t(scale(t(dat[cg,]))) > n[n>2]=2 > n[n< -2]= -2 > n[1:4,1:4] > pheatmap(n,show_colnames = F,show_rownames = F)加注释
> ac=data.frame(g=group_list) > rownames(ac)=colnames(n) > pheatmap(n,show_colnames = F,show_rownames = F,annotation_col = ac)
P7 ID转换
数据准备：Ensembl gene ID （似乎需要科学上网）
-
01:48 读入
> options(stringsAsFactors = F) > a=read.table('ensembl.txt') -
02:56 取 ‘.’ 号之前的部分（基因名）
> strsplit('ENSG00000000003.14','[.]') > strsplit('ENSG00000000003.14','[.]')[[1]] > strsplit('ENSG00000000003.14','[.]')[[1]][1] -
03:48 批量取基因名
> library(stringr) > str_split(a$V1,'[.]')> unlist(str_split(a$V1,'[.]'))这一步得到的是
character所以需要在之前的代码中增加一个参数
> unlist(str_split(a$V1,'[.]',simplify = T)) -
05:13 批量取第一列（基因名）
> a$ensembl_id=str_split(a$V1,'[.]',simplify = T)[,1] -
05:50 转换
> library(org.Hs.eg.db) > g2s=toTable(org.Hs.egSYMBOL) > g2e=toTable(org.Hs.egENSEMBL) -
07:15 关联
> b=merge(a,g2e,by='ensembl_id',all.x=T) > d=merge(b,g2s,by='gene_id',all.x=T) -
08:46 调序，调成和 ensembl.txt 中一致的顺序
搜索 d 中 ensembl_id 重复的基因
> table(table(d$ensembl_id)>1) > table(d$ensembl_id)[table(d$ensembl_id)>1]去重，调序
> d=d[order(d$V1),] > d=d[!duplicated(d$V1),] > d=d[match(a$V1,d$V1),]
P8 任意基因任意癌症表达量分组的生存分析
-
02:19 读取
> options(stringsAsFactors = F) > a=read.table('LGG_93663_50_50.csv',header = T,sep = ',',fill = T) -
03:23
ggbetweenstats作图> library(ggstatsplot) > ggbetweenstats(data = dat,x=Group,y=Expression) -
05:19 生存曲线绘制
> 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 = TRUE) -
06:09 调整图片
> ggsurvplot(sfit,palette = c("#E7B800","#2E9FDF"), + risk.table=TRUE,xlab="Time in months", + ggthemes=theme_light(), + ncensor.plot=TRUE)
P9 任意基因任意癌症表达量和临床形状关联
-
04:35 读取文件，更改列名
> a=read.table('plot-data-ARHGAP18-TCGA-0V-cbioportal.txt',header = T,sep = '\t',fill = T) > colnames(a)=c('id','stage','gene','mut') -
05:37
ggbetweenstats作图> library(ggstatsplot) > ggbetweenstats(data=dat,x=stage,y=gene)
P10 表达矩阵的样本的相关性
-
02:05 从
airway获得表达矩阵> library(airway) > data(airway) > exprSet=assay(airway)查看维度
> dim(exprSet) -
03:07 查看样本相关性，绘制热图
> cor(exprSet[,1],exprSet[,2]) > cor(exprSet) > pheatmap::pheatmap(cor(exprSet)) -
05:42 给热图添加注释
> group_list=colData(airway)[,3] > tmp=data.frame(g=group_list) > rownames(tmp)=colnames(cor(exprSet)) > pheatmap::pheatmap(cor(exprSet),annotation_col = tmp) -
07:33 数据筛选
显然第2行的数据不可用
> x=exprSet[1,] > table(x>1) > sum(x>1) > sum(x>1)>5 > x=exprSet[2,] > sum(x>1)>5通过以上代码找到判断方法
> apply(exprSet, 1,function(x) sum(x>1)>5) > exprSet=exprSet[apply(exprSet, 1,function(x) sum(x>1)>5),] -
10:24 进一步缩小表达矩阵
> exprSet=log(edgeR::cpm(exprSet)+1) > exprSet=exprSet[names(sort(apply(exprSet, 1, mad),decreasing = T)[1:500]),] -
11:07 得到新的相关性矩阵，重新绘制热图
> M=cor(log2(exprSet+1)) > pheatmap::pheatmap(M,annotation_col = tmp)
P11 芯片表达矩阵下游分析
-
00:46 获得一个表达矩阵
> suppressPackageStartupMessages(library(CLL)) > data(sCLLex) > exprSet=exprs(sCLLex) > samples=sampleNames(sCLLex) > pdata=pData(sCLLex) > group_list=as.character(pdata[,2]) > dim(exprSet) > exprSet[1:5,1:5]> boxplot(exprSet) -
03:21 单个差异分析
> boxplot((exprSet[1,]~group_list)) > t.test(exprSet[1,]~group_list) -
04:01 使用
limma进行批量差异分析设计矩阵
> suppressMessages(library(limma)) > design <- model.matrix(~0+factor(group_list)) > colnames(design)=levels(factor(group_list)) > rownames(design)=colnames(exprSet) > design构造分组
> contrast.matrix <- makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design) > contrast.matrix差异分析
> fit <- lmFit(exprSet,design) > fit2 <- contrasts.fit(fit,contrast.matrix) > fit2 <- eBayes(fit2) > tempOutput=topTable(fit2,coef = 1,n=Inf) > nrDEG=na.omit(tempOutput)
P12 RNA-Seq 表达矩阵差异分析
-
00:39 加载表达矩阵
> 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 ),] > table(group_list) > exprSet[1:4,1:4] > boxplot(log(exprSet+1)) -
03:11 使用
DESeq2进行差异分析> if(T){ + library(DESeq2) + + (colData <- data.frame(row.names = colnames(exprSet), + group_list=group_list)) + dds <- DESeqDataSetFromMatrix(countData=exprSet, + colData=colData, + design=~group_list) + tmp_f='airway_DESseq2-dds.Rdata'} > if(!file.exists(tmp_f)){ + dds <- DESeq(dds) + save(dds,file = tmp_f) + } estimating size factors estimating dispersions gene-wise dispersion estimates mean-dispersion relationship final dispersion estimates fitting model and testing > load(file = tmp_f) > res <- results(dds, + contrast=c("group_list","trt","untrt")) > resOrdered <- res[order(res$padj),] > head(resOrdered) > DEG=as.data.frame(resOrdered) > DESeq2_DEG=na.omit(DEG) -
04:32 下游分析
> nrDEG=DESeq2_DEG[,c(2,6)] > colnames(nrDEG)=c('logFC','P.Value') > attach(nrDEG) > plot(logFC,-log10(P.Value))> library(ggpubr) > df=nrDEG > df$v= -log10(P.Value) > ggscatter(df,x="logFC",y="v",size = 0.5)更改设置
> df$g=ifelse(df$P.Value>0.01,'stable', ifelse(df$logFC > 1.5,'up', ifelse(df$logFC < -1.5,'down','stable'))) > table(df$g) > df$name=rownames(df) > ggscatter(df,x="logFC",y="v",size = 0.5,color = 'g', palette = c("#00AF88","#E7B800","#FC4E07"))
最后，向大家隆重推荐生信技能树的一系列干货！
- 生信技能树全球公益巡讲：https://mp.weixin.qq.com/s/E9ykuIbc-2Ja9HOY0bn_6g
- B站公益74小时生信工程师教学视频合辑：https://mp.weixin.qq.com/s/IyFK7l_WBAiUgqQi8O7Hxw
- 招学徒：https://mp.weixin.qq.com/s/KgbilzXnFjbKKunuw7NVfw
