接下来介绍R语言：
【生信技能树】生信人应该这样学R语言
在你开始R之旅前，建议你看看下面这两个
1. 介绍R语言及Rstudio
- 了解R，Rstudio及R包;安装的包在packages中检查
.libPaths() #找安装路径
- 帮助文档，帮忙看路径
?substring
- 定位文件、设置文件位置
2. R语言基础变量讲解
- 重点就是理解: 五种变量结构（class属性）
- 我也曾经写过一点这方面笔记：R语言学习笔记
index1 = grep('RNA-Seq', a$Assay_Type)
index2 = grepl('RNA-Seq', a$Assay_Type)
b = a[index1,] # 下标
b = a[index2,] # 索引
3. 外部数据导入导出
- 去GEO数据库下载表达矩阵GSE17215
-
用excel打开后
> 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]
b = log2(b)
pheatmap::pheatmap(b[1:10,])
save(b,file = 'b_input.Rdata')
load(file = 'b_input.Rdata')
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
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,])
sample(1:nrow(b),50)
pheatmap::pheatmap(b[sample(1:nrow(b),50),])
+ abs, sqrt :戒对治，平方根
+ Log, log10, log2, exp: 对数与指数函数
+ Sin, cos, tan, acos, atan, atan2: 三角函数
+ sinh, cosh, tanh, asinha, acosh, aranh：双曲函数
+ 集合运算，reshape, merge总结
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=rnorm(100) #随机产生100个服从正态分布的数
?rnorm
dim(a1)=c(5,20) # 添加维度属性，矩阵matrix
pheatmap(a1)
a2=rnorm(100)+2
dim(a2)=c(5,20)
pheatmap(a2)
pheatmap(cbind(a1,a2),cluster_cols = F)
a3=cbind(a1,a2) #横向拼接
a4=rbind(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)
#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 <- 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 = TRUE,pval = TRUE,
conf.int = TRUE,xlab="Time in months",
ggtheme = theme_light(),
ncensor.plot=TRUE)
#ggsave('survival_ARHGAP18_in_LGG.png')
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)
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
> dim(exprSet)
[1] 64102 8
> cor(exprSet[,1],exprSet[,2]) #查看相关性
[1] 0.9632268
- 相关性真的较好
- 存在较多的0，那么大部分相关性由这些0决定
直接查看列样本之间基因表达的相关性
> cor(exprSet)
SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517
SRR1039508 1.0000000 0.9632268 0.9431333 0.8978772 0.9476515 0.9014495
SRR1039509 0.9632268 1.0000000 0.9323701 0.9428196 0.9202060 0.9290967
SRR1039512 0.9431333 0.9323701 1.0000000 0.9592993 0.9291853 0.9203183
SRR1039513 0.8978772 0.9428196 0.9592993 1.0000000 0.8834640 0.9498761
SRR1039516 0.9476515 0.9202060 0.9291853 0.8834640 1.0000000 0.9313570
SRR1039517 0.9014495 0.9290967 0.9203183 0.9498761 0.9313570 1.0000000
SRR1039520 0.9603249 0.9398292 0.9827034 0.9413725 0.9298860 0.9030552
SRR1039521 0.9023008 0.9463991 0.9525262 0.9853059 0.8724714 0.9291902
SRR1039520 SRR1039521
SRR1039508 0.9603249 0.9023008
SRR1039509 0.9398292 0.9463991
SRR1039512 0.9827034 0.9525262
SRR1039513 0.9413725 0.9853059
SRR1039516 0.9298860 0.8724714
SRR1039517 0.9030552 0.9291902
SRR1039520 1.0000000 0.9559571
SRR1039521 0.9559571 1.0000000
- 这样直接查看是没有意义的，通常实验设计是分有control组，实验组的
-
通过热图看看
- 添加样本分组信息
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)
筛选没有意义的基因：表达量大于1，样本数小于5的
> x=exprSet[1,]
> x
SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520
679 448 873 408 1138 1047 770
SRR1039521
572
> x.1
Error: object 'x.1' not found
> x>1
SRR1039508 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=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排除0
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)
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 SRR1039513
ENSG00000000003 679 448 873 408
ENSG00000000419 467 515 621 365
ENSG00000000457 260 211 263 164
ENSG00000000460 60 55 40 35
dev.off()
boxplot(exprSet)
13. R语言小作业
准备工作--安装所需的包
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./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 CRAN
if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna./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 install
for (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 symble
g2e=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')
2.根据探针名找对应symbol ID
1053_at
117_at
121_at
1255_g_at
1316_at
1320_at
1405_i_at
1431_at
1438_at
1487_at
1494_f_at
1598_g_at
160020_at
1729_at
177_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
rm(list=ls())
options(stringsAsFactors = F)
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=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，2
boxplot(exprSet[probe_tp53[i],] ~ pdata$Disease)
4.BRCA1基因表达量
- 找到BRCA1基因在TCGA数据库的乳腺癌数据集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表达情况
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=a
library(ggstatsplot)
ggbetweenstats(data = dat,x=subtype,y=expression)
5 .TP53 生存分析
rm(list=ls())
options(stringsAsFactors = F)
a=read.table('BRCA_7157_50_50.csv',sep = ',',fill=T,header = T)
dat=a
library(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)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int = F, pval = T)
ggsave("survival_TP53_in_BRCA_taga.png")
- 和官网给的分析相差不大0.139的0.14
- 也可以进一步花复杂的图
# complex survplot
ggsurvplot(sfit,palette = c("orange", "navyblue"),
risk.table = T, pval = T,
conf.int = T, xlab = "Time of datys",
ggtheme = theme_light(),
ncensor.plot = T)
ggbetweenstats(data = dat,
x = Group,
y = Expression)
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)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int = F, pval = T)
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平台，所以下载到的是一个含有一个元素的list
a=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$symbol
dat[1:4,1:4]
dim(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')
7.相关性计算和热图
# 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平台，所以下载到的是一个含有一个元素的list
a=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.找到芯片平台对应的注释包
# 8 题
rm(list=ls())
options(stringsAsFactors = F)
options()$repos
options()$BioC_mirror
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna./CRAN/"))
BiocManager::install("hugene10sttranscriptcluster.db",ask = F,update = F) #注释包都加上.db
options()$repos
options()$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平台，所以下载到的是一个含有一个元素的list
a=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平台，所以下载到的是一个含有一个元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
# 提取分组信息
group_list=unlist(lapply(pd$title, function(x){
strsplit(x,' ')[[1]][4]
}))
exprSet=dat
exprSet[1:4,1:4] # limma接受log后的表达矩阵，这一步需要检查一下
#DEO by 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<-makeContrasts("progres.-stable",levels = design)
contrast.matrix
## 这个矩阵声明，我们要把progres.组和stable进行差异分析
## step1
fit2<-contrasts.fit(fit,contrast.matrix) #这一步很重要
fit2<- eBayes(fit2) ##default no trend
## eBayes() with trend=TRUE
## step3
tempOutput=topTable(fit2.coef=1.n=Inf)
nrDEG=na.omit(tempOutput)
#write.csv(nrDEG2,"limma_noted,results.csv",quote=F)
head(nrDEG)
|
