生信技能树联盟入门视频：
生信技能树视频课程学习路径，这么好的视频还免费！专题历史目录：
3个学生的linux视频学习笔记接下来介绍R语言：
【生信技能树】生信人应该这样学R语言
R语言
在你开始R之旅前，建议你看看下面这两个
1.
课前准备：周末生物信息学培训班准备工作2.
致学员：无敌无脑的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打开后
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]
经过上面处理后如下；
image.png
画个热图看看看
b = log2(b) pheatmap::pheatmap(b[1:10,])
热图
保存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] NAWarning 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,])
image.png
随机选取50ge，sample函数
sample(1:nrow(b),50)pheatmap::pheatmap(b[sample(1:nrow(b),50),])
随机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个服从正态分布的数?rnormdim(a1)=c(5,20) # 添加维度属性，矩阵matrixpheatmap(a1)
a1
产生a2+热图
a2=rnorm(100)+2dim(a2)=c(5,20)pheatmap(a2)
a2
将a1,a2横向拼接并画图
pheatmap(cbind(a1,a2),cluster_cols = F)a3=cbind(a1,a2) #横向拼接a4=rbind(a1,a2) #纵向拼接
a1,a2
6. 选取差异明显的基因的表达量矩阵绘制热图
rm(list = ls()) #魔幻操作，一键清空library(pheatmap)a1 = rnorm(100)dim(a1) = c(5,20)pheatmap(a1)a2 = rnow(100)+2dim(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. 任意基因任意癌症表达量分组的生存分析
在线制作生存曲线的工具，oncolnc【
http://www.oncolnc.org/】
也可以下载原始数据自己分析
#读取csv文件rm(list = ls()) #清除所有变量options(stringsAsFactors = F)a=read.table('LGG_93663_50_50.csv',header = T,sep =',',fill = T)# 后续画图colnames(a)dat=a# 图ggbetweenstatslibrary(ggstatsplot)ggbetweenstats(data = dat,x=Group,y=Expression)# 图ggsurvplotlibrary(ggplot2)library(survival)library(survminer)table(dat$Status)dat$Status <- ifelse(dat$Status == "Dead", 1, 0)sfit <- survfit(Surv(Days,Status)~Group, data=dat)sfitsummary(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')
ggbetweenstats
ggsurvplot
9. 任意基因任意癌症表达量和临床性状关联
某个基因在某个癌症的表达量，关联临床信息
Cbioportal【
http://www.cbioportal.org/】，同样下载ARHGAP18
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=alibrary(ggstatsplot)ggbetweenstats(data = dat,x=stage,y=gene)
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
image.png
> dim(exprSet)[1] 64102 8> cor(exprSet[,1],exprSet[,2]) #查看相关性[1] 0.9632268
可知第一列和第二列的相关性较好：
相关性真的较好
存在较多的0，那么大部分相关性由这些0决定
直接查看列样本之间基因表达的相关性
> cor(exprSet) SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517SRR1039508 1.0000000 0.9632268 0.9431333 0.8978772 0.9476515 0.9014495SRR1039509 0.9632268 1.0000000 0.9323701 0.9428196 0.9202060 0.9290967SRR1039512 0.9431333 0.9323701 1.0000000 0.9592993 0.9291853 0.9203183SRR1039513 0.8978772 0.9428196 0.9592993 1.0000000 0.8834640 0.9498761SRR1039516 0.9476515 0.9202060 0.9291853 0.8834640 1.0000000 0.9313570SRR1039517 0.9014495 0.9290967 0.9203183 0.9498761 0.9313570 1.0000000SRR1039520 0.9603249 0.9398292 0.9827034 0.9413725 0.9298860 0.9030552SRR1039521 0.9023008 0.9463991 0.9525262 0.9853059 0.8724714 0.9291902 SRR1039520 SRR1039521SRR1039508 0.9603249 0.9023008SRR1039509 0.9398292 0.9463991SRR1039512 0.9827034 0.9525262SRR1039513 0.9413725 0.9853059SRR1039516 0.9298860 0.8724714SRR1039517 0.9030552 0.9291902SRR1039520 1.0000000 0.9559571SRR1039521 0.9559571 1.0000000
这样直接查看是没有意义的，通常实验设计是分有control组，实验组的
通过热图看看
image.png
添加样本分组信息
exprSet=assay(airway)colnames(exprSet)group_list=colData(airway)[,3]tmp=data.frame(g=group_list)rownames(tmp)=colnames(exprSet)pheatmap::pheatmap(cor(exprSet),annotation_col = tmp)
image.png
筛选没有意义的基因：表达量大于1，样本数小于5的
先看看第一行的情况
> x=exprSet[1,]> xSRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 679 448 873 408 1138 1047 770 SRR1039521 572 > x.1Error: object 'x.1' not found> x>1SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 TRUE TRUE TRUE TRUE TRUE TRUE TRUE SRR1039521 TRUE > sum(x>1)[1] 8# TRUE的逻辑值为1，FLASE的逻辑值为0
按照上面的规则进行筛选
> dim(exprSet)[1] 64102 8> exprSet=exprSet[apply(exprSet, 1,function(x)sum(x>1)>5),]> dim(exprSet)[1] 19481 8
进一步对exprSet处理
exprSet=log(edgeR::cpm(exprSet)+1) #edgeR::cpm(exprSet)去除文库大小差异dim(exprSet)exprSet=exprSet[names(sort(apply(exprSet, 1,mad),decreasing = T)[1:500]),]dim(exprSet)M=cor(log2(exprSet+1)) # +1排除0tmp=data.frame(g=group_list)rownames(tmp)=colnames(M)pheatmap::pheatmap(M,annotation_col = tmp)
image.png
view() dim() #看维度
11. 芯片表达矩阵下游分析
12. RNA-seq表达矩阵差异分析
芯片的表达矩阵 exprSet = exprs(sCLLex)
RNA-seq表达矩阵 exprSet = assay(ariway) # assay获取表达矩阵
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),]table(group_list)> table(group_list)group_list trt untrt 4 4 exprSet[1:4,1:4]> exprSet[1:4,1:4] SRR1039508 SRR1039509 SRR1039512 SRR1039513ENSG00000000003 679 448 873 408ENSG00000000419 467 515 621 365ENSG00000000457 260 211 263 164ENSG00000000460 60 55 40 35dev.off()boxplot(exprSet)
image.png
多个探针对应同一个基因取最大表达量探针极简代码 13. R语言小作业
image.png
准备工作--安装所需的包
cran_packages <- c('tidyverse', 'ggpubr', 'ggstatsplot'Biocductor_packages <- c('org.Hs.eg.db', 'hgu133a.db', 'CLL', 'hgu95av2.db', 'survminer', 'survival', 'hugene10sttranscriptcluster', 'limma')if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")for (pkg in cran_packages){ if (! require(pkg,character.only=T) ) { install.packages(pkg,ask = F,update = F) require(pkg,character.only=T) }}# first prepare BioManager on CRANif(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)if(length(getOption("BioC_mirror"))==0) options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")# use BiocManager to installfor (pkg in Biocductor_packages){ if (! require(pkg,character.only=T) ) { BiocManager::install(pkg,ask = F,update = F) require(pkg,character.only=T) }}
1.找到指定ensembl ID与symbol的对应关系
> ENSG00000000003.13> ENSG00000000005.5> ENSG00000000419.11> ENSG00000000457.12> ENSG00000000460.15> ENSG00000000938.11
思路：在注释包中有gene_id与symbol、gene_id与ensembl_id的对应关系。
#将以上基因id保存在a.txt，存放于工作目录下。rm(list=ls())options(stringsAsFactors = F)read.table('R 10/R1')a=read.table('R 10/R1')library(org.Hs.eg.db)g2s=toTable(org.Hs.egSYMBOL) #gene ID symbleg2e=toTable(org.Hs.egENSEMBL)library(stringr)unlist(lapply(a$V1,function(x){ strsplit(as.character(x),'[.]')[[1]][1]}))a$ensembl_id=unlist(lapply(a$V1,function(x){ strsplit(as.character(x),'[.]')[[1]][1]}))tmp=merge(a,g2e,by='ensembl_id')tmp=merge(tmp,g2s,by='gene_id')
答案1
2.根据探针名找对应symbol ID
1053_at117_at121_at1255_g_at1316_at1320_at1405_i_at1431_at1438_at1487_at1494_f_at1598_g_at160020_at1729_at177_at
思路：找到注释包中探针与symbol的对应关系然后取子集
rm(list=ls())options(stringsAsFactors = F)a=read.table('e2')library(hgu133a.db)ids=toTable(hgu133aSYMBOL)head(ids)colnames(a)='probe_id'tmp=merge(ids,a,by='probe_id')# 第二种方法match(a$probe_id,ids$probe_id)tmp2=ids[match(a$probe_id,ids$probe_id),]
附：判断得到的两组结果是否一致
# 法一：identical(tmp1,tmp2) #返回逻辑值# 法二：dplyr::setdiff(tmp1,tmp2) #返回两组的差别【没差就返回空】
3.根据symbol找基因表达量并作图
找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量，并且绘制在 progres.-stable分组的boxplot图,并通过 ggpubr 进行美化。
+探针的三大内容：表达矩阵assay/exprs、探针信息featureData、样本信息phenoData
symbol
rm(list=ls())options(stringsAsFactors = F)suppressPackageStartupMessages(library(CLL))data(sCLLex)sCLLexexprSet=exprs(sCLLex) library(hgu95av2.db) #这个包找到目的基因ids=toTable(hgu95av2SYMBOL)exprSet <- exprs(sCLLex) #探针的表达量pdata <- pData(sCLLex) #sampleID与disease的对应关系p2s <- toTable(hgu95av2SYMBOL) #探针与symbol的对应关系p3 <- filter(p2s,symbol=='TP53')# boxplot [find TP53 has 3 probe IDs]probe_tp53 <- c("1939_at","1974_s_at","31618_at")i = 3 #可换1，2boxplot(exprSet[probe_tp53[i],] ~ pdata$Disease)
image.png
4.BRCA1基因表达量
找到BRCA1基因在TCGA数据库的乳腺癌数据集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表达情况
提示：使用
http://www.cbioportal.org/index.do 定位数据集：
http://www.cbioportal.org/datasets该基因有四个亚型，用ggbetweenstats作图比较一下。
image.png
rm(list=ls())options(stringsAsFactors = F)a=read.table('plot.txt',sep = '\t',fill=T,header = T)colnames(a)=c('id','subtype','expression','mut')dat=alibrary(ggstatsplot)ggbetweenstats(data = dat,x=subtype,y=expression)
image.png
5 .TP53 生存分析
找到TP53基因在TCGA数据库的乳腺癌数据集的表达量分组看其是否影响生存
提示使用：
http://www.oncolnc.org/思路：生存分析，TP53表达量分为高低两组做图比较
官网生存分析
rm(list=ls())options(stringsAsFactors = F)a=read.table('BRCA_7157_50_50.csv',sep = ',',fill=T,header = T)dat=alibrary(ggstatsplot)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)sfitsummary(sfit)ggsurvplot(sfit, conf.int = F, pval = T)ggsave("survival_TP53_in_BRCA_taga.png")
TP53
和官网给的分析相差不大0.139的0.14
也可以进一步花复杂的图
# complex survplotggsurvplot(sfit,palette = c("orange", "navyblue"), risk.table = T, pval = T, conf.int = T, xlab = "Time of datys", ggtheme = theme_light(), ncensor.plot = T)
image.png
ggbetweenstats(data = dat, x = Group, y = Expression)
image.png
b=read.table('BRCA.txt',sep = '\t',fill=T,header = T) #txt，用'\t'；excel用','。colnames(b)=c('Patient','subtype','expression','mut')head(b)b$Patient=substring(b$Patient,1,12)tmp=merge(a,b,by='Patient')table(tmp$subtype)dat=tmp[tmp$subtype=='BRCA_Basal',]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)sfitsummary(sfit)ggsurvplot(sfit, conf.int = F, pval = T)
image.png
6.从表达矩阵中提取指定基因画热图
下载数据集GSE17215的表达矩阵并且提取下面的基因画热图
ACTR3B ANLN BAG1 BCL2 BIRC5 BLVRA CCNB1 CCNE1 CDC20 CDC6 CDCA1 CDH3 CENPF CEP55 CXXC5 EGFR ERBB2 ESR1 EXO1 FGFR4 FOXA1 FOXC1 GPR160 GRB7 KIF2C KNTC2 KRT14 KRT17 KRT5 MAPT MDM2 MELK MIA MKI67 MLPH MMP11 MYBL2 MYC NAT1 ORC6L PGR PHGDH PTTG1 RRM2 SFRP1 SLC39A6 TMEM45B TYMS UBE2C UBE2T
提示：根据基因名拿到探针ID，缩小表达矩阵绘制热图，没有检查到的基因直接忽略即可。
rm(list=ls())options(stringsAsFactors = F)#下载和表达矩阵library(GEOquery)GSE17125 <- "GSE17215"f=GSE17125#这个包需要注意两个配置，一般来说自动化的配置是足够的。#Setting options('download.file.method.GEOquery'='auto')#Setting options('GEOquery.inmemory.gpl'=FLASE)if(!file.exists(f)){ gset <- getGEO('GSE17215', destdir = ".", getGPL = F, AnnotGPL = F) save(gset, file = f)}load('GSE17215') #载入数据class(gset)length(gset)class(gset[[1]])# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的lista=gset[[1]]dat=exprs(a)dim(dat)#改基因library(hgu133a.db)ids=toTable(hgu133aSYMBOL)head(ids)dat=dat[ids$probe_id,]#dat=dat[1:4,1:4]ids$median=apply(dat,1,median)ids=ids[order(ids$symbol,ids$median,decreasing = T),]ids=ids[!duplicated(ids$symbol),]dat=dat[ids$probe_id,]rownames(dat)=ids$symboldat[1:4,1:4]dim(dat)
dat
画出上面罗列基因的热图
ng='ACTR3B ANLN BAG1 BCL2 BIRC5 BLVRA CCNB1 CCNE1 CDC20 CDC6 CDCA1 CDH3 CENPF CEP55 CXXC5 EGFR ERBB2 ESR1 EXO1 FGFR4 FOXA1 FOXC1 GPR160 GRB7 KIF2C KNTC2 KRT14 KRT17 KRT5 MAPT MDM2 MELK MIA MKI67 MLPH MMP11 MYBL2 MYC NAT1 ORC6L PGR PHGDH PTTG1 RRM2 SFRP1 SLC39A6 TMEM45B TYMS UBE2C UBE2T'ng=strsplit(ng,' ')[[1]]# 过滤不再基因名的table(ng %in% rownames(dat))ng=ng[ng %in% rownames(dat)]dat=dat[ng,]dat=log2(dat)pheatmap::pheatmap(dat,scale = 'row')
image.png
7.相关性计算和热图
下载数据集GSE24673的表达矩阵计算样本的相关性并且绘制热图，需要标记上样本分组信息
相关性分析：
# 7题rm(list=ls())options(stringsAsFactors = F)f='GSE24673_eSet.Rdata'#下载和表达矩阵library(GEOquery)#这个包需要注意两个配置，一般来说自动化的配置是足够的。#Setting options('download.file.method.GEOquery'='auto')#Setting options('GEOquery.inmemory.gpl'=FLASE)if(!file.exists(f)){ gset <- getGEO('GSE24673', destdir = ".", getGPL = F, AnnotGPL = F) save(gset, file = f)}load('GSE24673_eSet.Rdata') #载入数据class(gset)length(gset)class(gset[[1]])# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的lista=gset[[1]]dat=exprs(a)dim(dat)expr[1:4,1:4]pdata <- pData(geo[[1]]) # 自己根据pdata第八列做一个分组信息矩阵grp <- c('rbc','rbc','rbc', 'rbn','rbn','rbn', 'rbc','rbc','rbc', 'normal','normal')grp_df <- data.frame(group=grp)rownames(grp_df) <- colnames(expr)new_cor <- cor(expr)pheatmap::pheatmap(new_cor, annotation_col = grp_df)
8.找到芯片平台对应的注释包
找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 对应的R的bioconductor注释包，并且安装它！
https://mp.weixin.qq.com/s/sVSsV2fWZOQmNd72Vb3SmQhttps://www.jianshu.com/p/40b560755cdf# 8 题rm(list=ls())options(stringsAsFactors = F)options()$reposoptions()$BioC_mirror options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))BiocManager::install("hugene10sttranscriptcluster.db",ask = F,update = F) #注释包都加上.dboptions()$reposoptions()$BioC_mirror
9.找到指定探针和对应的基因
下载数据集GSE42872的表达矩阵，并且分别挑选出所有样本的(平均表达量/sd/mad/)最大的探针，并且找到它们对应的基因。
# 9题rm(list=ls())options(stringsAsFactors = F)f='GSE42872_eSet.Rdata'#下载和表达矩阵library(GEOquery)#这个包需要注意两个配置，一般来说自动化的配置是足够的。#Setting options('download.file.method.GEOquery'='auto')#Setting options('GEOquery.inmemory.gpl'=FLASE)if(!file.exists(f)){ gset <- getGEO('GSE42872', destdir = ".", getGPL = F, AnnotGPL = F) save(gset, file = f)}load('GSE42872_eSet.Rdata') #载入数据class(gset)length(gset)class(gset[[1]])# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的lista=gset[[1]]dat=exprs(a)dim(dat)pd=pData(a)# 平均表达量/sd/mad 最大探针sort(apply(dat, 1, mean),decreasing = T)[1]sort(apply(dat, 1, sd),decreasing = T)[1]sort(apply(dat, 1, mad),decreasing = T)[1]> # 平均表达量/sd/mad 最大探针> sort(apply(dat, 1, mean),decreasing = T)[1] 7978905 14.53288 > sort(apply(dat, 1, sd),decreasing = T)[1] 8133876 3.166429 > sort(apply(dat, 1, mad),decreasing = T)[1] 8133876 4.268561
10.limma 差异分析
下载数据集GSE42872的表达矩阵，并且根据分组使用limma做差异分析，得到差异结果矩阵
接第九题，得到表型信息，然后用limma做差异分析
# 10题rm(list=ls())options(stringsAsFactors = F)f='GSE42872_eSet.Rdata'#下载和表达矩阵library(GEOquery)#这个包需要注意两个配置，一般来说自动化的配置是足够的。#Setting options('download.file.method.GEOquery'='auto')#Setting options('GEOquery.inmemory.gpl'=FLASE)if(!file.exists(f)){ gset <- getGEO('GSE42872', destdir = ".", getGPL = F, AnnotGPL = F) save(gset, file = f)}load('GSE42872_eSet.Rdata') #载入数据class(gset)length(gset)class(gset[[1]])# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的lista=gset[[1]]dat=exprs(a)dim(dat)pd=pData(a)# 提取分组信息group_list=unlist(lapply(pd$title, function(x){ strsplit(x,' ')[[1]][4]}))exprSet=datexprSet[1:4,1:4] # limma接受log后的表达矩阵，这一步需要检查一下#DEO by limmasuppressMessages(library(limma))design<- model.matrix(~0+factor(group_list))colnames(design)=levels(factor(group_list))rownames(design)=colnames(exprSet)designcontrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)contrast.matrix<-makeContrasts("progres.-stable",levels = design)contrast.matrix## 这个矩阵声明，我们要把progres.组和stable进行差异分析## step1fit2<-contrasts.fit(fit,contrast.matrix) #这一步很重要fit2<- eBayes(fit2) ##default no trend## eBayes() with trend=TRUE## step3tempOutput=topTable(fit2.coef=1.n=Inf)nrDEG=na.omit(tempOutput)#write.csv(nrDEG2,"limma_noted,results.csv",quote=F)head(nrDEG)
