数据科学与R语言R语言R语言做生信

生信人的R语言

2018-12-29  本文已影响16人  泥人吴

生信技能树联盟入门视频:生信技能树视频课程学习路径,这么好的视频还免费!

专题历史目录:3个学生的linux视频学习笔记

接下来介绍R语言:

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

R语言
在你开始R之旅前,建议你看看下面这两个

1. 介绍R语言及Rstudio

    getwd()
    setwd()  

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


    GEO17215
> 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. 中级变量操作

函数
# 最大值
> 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,])
image.png
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总结
    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)
a1
a2=rnorm(100)+2
dim(a2)=c(5,20)
pheatmap(a2)
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)
    #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()  #去重

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')
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
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=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表达矩阵差异分析

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 <- 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 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
#将以上基因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
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 <- 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)
image.png

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)
image.png

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")
TP53
# 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 <- ifelse(dat$Status == "Dead", 1, 0)

sfit <- survfit(Surv(Days,Status)~Group, data=dat)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int = F, pval = T)
image.png

6.从表达矩阵中提取指定基因画热图

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
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)
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.相关性计算和热图

# 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.tsinghua.edu.cn/CRAN/"))
BiocManager::install("hugene10sttranscriptcluster.db",ask = F,update = F) #注释包都加上.db
options()$repos
options()$BioC_mirror

9.找到指定探针和对应的基因

# 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 差异分析

# 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)
上一篇 下一篇

猜你喜欢

热点阅读