接下来介绍R语言：

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

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打开后

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] 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,])

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个服从正态分布的数

?rnorm

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

pheatmap(a1)

a1

产生a2+热图

a2=rnorm(100)+2

dim(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)+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

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')

ggbetweenstats

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)

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 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组，实验组的

通过热图看看

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,]

> 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处理

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)

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 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)

image.png

13. R语言小作业

image.png

准备工作--安装所需的包

cran_packages

'ggpubr',

'ggstatsplot'

Biocductor_packages

'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 CRAN

if(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 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')

答案1

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

symbol

rm(list=ls())

options(stringsAsFactors = F)

suppressPackageStartupMessages(library(CLL))

data(sCLLex)

sCLLex

exprSet=exprs(sCLLex)

library(hgu95av2.db) #这个包找到目的基因

ids=toTable(hgu95av2SYMBOL)

exprSet

pdata

p2s

p3

# boxplot [find TP53 has 3 probe IDs]

probe_tp53

i = 3 #可换1，2

boxplot(exprSet[probe_tp53[i],] ~ pdata$Disease)

image.png

4.BRCA1基因表达量

找到BRCA1基因在TCGA数据库的乳腺癌数据集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表达情况

该基因有四个亚型，用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=a

library(ggstatsplot)

ggbetweenstats(data = dat,x=subtype,y=expression)

image.png

5 .TP53 生存分析

找到TP53基因在TCGA数据库的乳腺癌数据集的表达量分组看其是否影响生存

思路：生存分析，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

sfit

sfit

summary(sfit)

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

ggsave("survival_TP53_in_BRCA_taga.png")

TP53

和官网给的分析相差不大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)

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

sfit

sfit

summary(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

f=GSE17125

#这个包需要注意两个配置，一般来说自动化的配置是足够的。

#Setting options('download.file.method.GEOquery'='auto')

#Setting options('GEOquery.inmemory.gpl'=FLASE)

if(!file.exists(f)){

gset

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)

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

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第八列做一个分组信息矩阵

grp

'rbn','rbn','rbn',

'rbc','rbc','rbc',

'normal','normal')

grp_df

rownames(grp_df)

new_cor

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.tsinghua.edu.cn/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

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

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

colnames(design)=levels(factor(group_list))

rownames(design)=colnames(exprSet)

design

contrast.matrix

contrast.matrix

contrast.matrix

## 这个矩阵声明，我们要把progres.组和stable进行差异分析

## step1

fit2

fit2

## 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)