1 介绍R及R Studio
右上角history，To Console选择到左下角，To Source选择到左上编辑器内。
画图没出图：dev.off
包的安装路径：.libPaths()
善用帮助文档：?函数
2 R语言基础变量
变量结构：向量、矩阵、数组、数据框、列表
当创建向量元素出现字符，那么全部元素都是字符型。
重要函数：class( ) str( ) 报错的时候查看元素性质
is. 函数：判断
as. 函数：转换成
索引：取数据 b[,3] 等价于 b[,F,F,T,F]
直接搜索：grep('RNA-Seq', a$Assay_type)、grepl() 结果判断为TRUE FALSE、
table() 计数
3 外部数据导入导出
把txt文件导入R
read.table('SraRunTable.txt')
R无法识别正确的列数
增加参数：
a <- read.table('SraRunTable.txt', sep = '\t')
绘制表头：
a <- read.table('SraRunTable.txt', header = TRUE,
sep = '\t')
image.png
更复杂的文件导入方式：!：有感叹号开头的内容不要
b <- read.table('R_bilibili-master/R语言小作业/GSE17215_series_matrix.txt.gz',
comment.char = "!",
header = TRUE,
sep = '\t')
把刚才文件导出成csv文件
write.csv(b,'GSE17215_series_matrix.csv')
把第一列名去掉
参数：row.name=FALSE
读取文件：
write.table(b,'tmp.csv',sep=',')
把第一列添加到列名，且删除原本第一列
rownames(b) <- b[,1]
b <- b[,-1]
在配置好R包后，绘制图像：
pheatmap::pheatmap(b[1:10,])
image.png
b数值太大了，取对数
b <- log2(b)
pheatmap::pheatmap(b[1:10,])
image.png
导出数据最好保存成Rdata文件
save(b, file = 'b_input.Rdata')
读取
load(file='b_input.Rdata')
4 中级变量操作
a <- read.table('SraRunTable.txt', header = TRUE,
sep = '\t')
以a为例子。
选取最大值和最小值
> sort(a$MBases)[1]
[1] 2128
> sort(a$MBases,decreasing = TRUE)[1]
[1] 19506
也可以直接用函数
max(a$MBases)
min(a$MBases)
也可以用统计学五分位数：
> fivenum(a$MBases)
[1] 2128.0 5513.0 6873.5 8415.0 19506.0
统计该变量小于5000的个数
> table(a$MBases < 5000)
FALSE TRUE
151 31
把小于5000的数据新建一个文件(列数据)
d <- a[a$MBases < 5000,]
发现吧绝大部分RNA-Seq的数据取出来了，我们看一下外显子和RNA-Seq的数据分布
boxplot(a$MBases~a$Assay_Type)
image.png
差异很大，所以要分组
wes <- a[a$Assay_Type=='Wxs',]
rna <- a[a$Assay_Type=='RNA-Seq',]
然后再根据每个组的数据分布来筛选
b <- read.table('R_bilibili-master/R语言小作业/GSE17215_series_matrix.txt.gz',
comment.char = "!",
header = TRUE,
sep = '\t')
rownames(b) <- b[,1]
b <- b[,-1]
b <- log2(b)
以b为例子。
计算b第一行的数据，但是查看格式
> str(b[1,])
'data.frame': 1 obs. of 6 variables:
$ GSM431121: num 8.91
$ GSM431122: num 9.22
$ GSM431123: num 11.4
$ GSM431124: num 11.3
$ GSM431125: num 11.4
$ GSM431126: num 11.4
是数据框格式，要转变为数字
> 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
想要批量计算均值？采用函数（相对最简单）
> 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
也可以建立循环，类似C语言的for循环
for (i in 1:nrow(b)) {
print(mean(as.numeric(b[1,])))
}
或者
apply(b, 1, function(x){
x <- as.numeric(b[1,])
mean(x)
})
自定义函数：选取每一列最大值
rowMax=function(x){
apply(x,1,max)
}
计算每一行方差值
apply(b,1,sd)
按照递减顺序，取前50，并拿到基因名字
cg=names(sort(apply(b,1,sd),decreasing=TRUE)[1:50])
可以随机取50个画热图
pheatmap::pheatmap(b[sample(1:nrow(b),50),])
用cg来画热图
pheatmap::pheatmap(b[cg,])
image.png
5 热图
生成向量，定义成矩阵进行画图
a1 <- rnorm(100)
dim(a1) <- c(5,20)
a2 <- rnorm(100)+2
dim(a2) <- c(5,20)
library(pheatmap)
pheatmap(a1,cluster_rows = FALSE,show_colnames = FALSE)
pheatmap(cbind(a1,a2))
pheatmap(cbind(a1,a2),show_rownames = FALSE, show_colnames = FALSE)
不排序
pheatmap(cbind(a1,a2),cluster_cols = F)
可以看到a2整体比a1大（左侧蓝a1偏小，右侧红a2偏大）
image.png
学习函数paste
> paste('a1',1:20,sep = '_')
[1] "a1_1" "a1_2" "a1_3" "a1_4" "a1_5" "a1_6" "a1_7" "a1_8" "a1_9"
[10] "a1_10" "a1_11" "a1_12" "a1_13" "a1_14" "a1_15" "a1_16" "a1_17" "a1_18"
[19] "a1_19" "a1_20"
b <- cbind(a1,a2)
b <- as.data.frame(b)
names(b) <- c(paste('a1',1:20,sep = '_'),paste('a2',1:20,sep = '_'))
pheatmap(b,cluster_cols = F)
image.png
图就有横坐标了。
增加group名
b <- cbind(a1,a2)
b <- as.data.frame(b)
pheatmap(cbind(a1,a2),cluster_cols = F)
names(b) <- c(paste('a1',1:20,sep = '_'),paste('a2',1:20,sep = '_'))
tmp <- data.frame(group=c(rep('a1',20), rep('a2',20)))
rownames(tmp) <- colnames(b)
pheatmap(b,annotation_col = tmp)
image.png
根据help文件里面的example，可以不断丰富热图。
6 选取差异明显的基因的表达量矩阵绘制热图
1转置
sort选取极值，用scale标准化。拉平大值和小值
n[n>2] <- 2
n[n<-2] <- -2
2再绘制
names(tail((sort(apply(dat, 1, sd))),1000)),t(转置)
7 id转换
首先导入ENSG文本
a <- read.table('e1.txt')
ENSG00000000003.13
分割点号：
> strsplit('ENSG00000000003.13', '[.]')
[[1]]
[1] "ENSG00000000003" "13"
取这个向量第一个元素，再取第一个元素
> strsplit('ENSG00000000003.13', '[.]')[[1]][1]
[1] "ENSG00000000003"
合并成功
引入循环，用stringr包
library(stringr)
a$ensembl_id=str_split(a$V1, '[.]', simplify = T)[,1]
我们就在a中加入了名为ensembl_id，合并数据的一列
image.png
载入包
library(org.Hs.eg.db)
ls("package:org.Hs.eg.db")
g2s=toTable(org.Hs.egSYMBOL);head(g2s)
g2e=toTable(org.Hs.egENSEMBL);head(g2e)
把a数据框与g2e关联
b <- merge(a,g2e,by='ensembl_id', all.x=T)
d <- merge(b,g2s,by='gene_id', all.x=T)
最后在d中恢复与a相同的排序
d <- d[order(d$V1),]
image.png
去除d中的重复序号
d <- d[!duplicated(d$V1)，]
合并d与a的顺序
d <- d[match(a$V1,d$V1),]
最后导出成csv文件
write.csv(d,'geneID2symbol.csv')
image.png
8 任意基因任意癌症表达量分组的生存分析
从生存分析网页工具下载数据
http://www.oncolnc.org/kaplan/?lower=50&upper=50&cancer=LGG&gene_id=93663&raw=ARHGAP18&species=mRNA
导入R
options(stringsAsFactors = F)
a <- read.table('R_bilibili-master/R语言小作业/LGG_93663_50_50.csv',
header = T,
sep = ',',
fill = T)
colnames(a)
dat <- a
绘制第一个图
library(ggstatsplot)
ggbetweenstats(data = dat, x= Group, y= Expression)
image.png
绘制生存曲线
table(dat$Status)
dat$Status <- ifelse(dat$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=dat)
summary(sfit)
ggsurvplot(sfit, conf.int = F, pval = T)
image.png
更进一步：
ggsurvplot(sfit, palette = c("#E7B800", "#2E9FDF"),
risk.table = T, pval = T,
conf.int = T,
xlab="Time in months",
ggtheme = theme_light(),
ncensor.plot=T)
保存图片
ggsave('survial_ARHGAP18_in_LGG.png')
image.png
9 任意基因任意癌症表达量和临床性关联
下载plot文件。
options(stringsAsFactors = F)
a=read.table('R_bilibili-master/R语言小作业/plot.txt',
sep = '\t',fill = T,header = T)
a <- a[,-5]
colnames(a)=c('id','subtype','expression','mut')
dat=a
library(ggstatsplot)
ggbetweenstats(data = dat, x = subtype, y = expression)
image.png
中间出现了一段报错，原因是Package 'PMCMRplus'没有安装，报错显示安装即可画出。
10 表达矩阵的样本相关性
获取airway数据：
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("airway")
library(airway)
options(stringsAsFactors = F)
library(airway)
data("airway")
exprSet <- assay(airway)
colnames(exprSet)
查看样本：
> dim(exprSet)
[1] 64102 8
8列，6万多个值。
第一列和第二列相关性
> cor(exprSet[,1],exprSet[,2])
[1] 0.9632268
group_list <- colData(airway)[,3]
tmp <- data.frame(g=group_list)
rownames(tmp) <- colnames(exprSet)
pheatmap::pheatmap(cor(exprSet),annotation_col = tmp)
image.png
删除数据里的空白数值，判断每一行是否都有数据。
exprSet <- exprSet[apply(exprSet,1, function(x) sum(x>1)>5),]
筛选后从原先
> dim(exprSet)
[1] 64102 8
6万多个基因变成
> dim(exprSet)
[1] 19481 8
接近2万个。
进一步整理数据，需要安装edgeR包：
> BiocManager::install("edgeR")
取mad数值最大前500个基因画图
exprSet <- log(edgeR::cpm(exprSet)+1)
exprSet <- exprSet[names(sort(apply(exprSet,1,mad),decreasing = T)[1:500]),]
M <- cor(log2(exprSet+1))
tmp <- data.frame(g=group_list)
rownames(tmp) <- colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)
调整后的图形：
image.png
11 芯片表达矩阵下游分析
获取CLL包的数据：
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("CLL")
正式开始：
library(CLL)
data(sCLLex)
sCLLex
exprSet <- exprs(sCLLex)
samples <- sampleNames(sCLLex)
pdata <- pData(sCLLex)
group_list <- as.character(pdata[,2])
dim(exprSet)
[1] 12625 22
exprSet[1:5,1:5]
CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL
1000_at 5.743132 6.219412 5.523328 5.340477 5.229904
1001_at 2.285143 2.291229 2.287986 2.295313 2.662170
1002_f_at 3.309294 3.318466 3.354423 3.327130 3.365113
1003_s_at 1.085264 1.117288 1.084010 1.103217 1.074243
1004_at 7.544884 7.671801 7.474025 7.152482 6.902932
获得差异矩阵
boxplot(exprSet)
构造比较矩阵
#DEG by limma
suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(exprSet)
design
progres. stable
CLL11.CEL 1 0
CLL12.CEL 0 1
CLL13.CEL 1 0
CLL14.CEL 1 0
CLL15.CEL 1 0
略
把progress组与stable组进行差异分析比较
contrast.matrix <- makeContrasts(paste0(unique(group_list),collapse="-"),levels=design)
contrast.matrix
Contrasts
Levels progres.-stable
progres. 1
stable -1
差异分析第一步：
fit <- lmFit(exprSet,design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
tempOutput <- topTable(fit2, coef=1, n=Inf)
nrDEG <- na.omit(tempOutput)
最后得到nrDEG，画火山图，富集分析等。
image.png
12 RNA-Seq表达矩阵差异分析
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)
方法：#DESeq2
if(T){
library(DESeq2)
(colData <- data.frame(row.names = colnames(exprSet),
group_list = group_list))
dss <- DESeqDataSetFromMatrix(countData = exprSet,
colData = colData,
design = ~group_list)
tmp_f <- 'airway_DESeq2-dss.Rdata'
if(!file.exists(tmp_f)){
dss <- DESeq(dss)
save(dss, file = tmp_f)
}
load(file=tmp_f)
res <- results(dds,
contrasts=c("group_list","trt","untrt"))
resOrdered <- res[order(res$padj),]
head(resOrdered)
DEG <- as.data.frame(resOrdered)
DESeq2_DEG <- na.omit(DEG)
nrDEG <- DESeq2_DEG[,c(2,6)]
colnames(nrDEG) <- c('log2FoldChange','pvalue')
}
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')
ggscatter(df, x="logFC", y='v',color='g', size=0.5,
palette=c("#00AFBB","#E7B800","#FC4E07"))
13 R语言习题
作业 1
根据R包org.Hs.eg.db找到下面ensembl 基因ID 对应的基因名(symbol)
ENSG00000000003.13
ENSG00000000005.5
ENSG00000000419.11
ENSG00000000457.12
ENSG00000000460.15
ENSG00000000938.11
提示：
library(org.Hs.eg.db)
g2s=toTable(org.Hs.egSYMBOL)
g2e=toTable(org.Hs.egENSEMBL)
rm(list = ls())
options(stringsAsFactors = F)
a <- read.table('R_bilibili-master/R语言小作业/e1.txt')
library(org.Hs.eg.db)
g2s=toTable(org.Hs.egSYMBOL)
g2e=toTable(org.Hs.egENSEMBL)
head(g2e)
library(stringr)
a$ensembl_id <- unlist(lapply(a$V1,function(x){
strsplit(x,'[.]')[[1]][1]
})
)
tmp = merge(a,g2e,by='ensembl_id')
tmp = merge(tmp,g2s,by='gene_id')
作业 2
根据R包hgu133a.db找到下面探针对应的基因名(symbol)。
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
提示：
library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
两种方法，tmp1修改a列名，tmp2不修改分别设置
rm(list = ls())
options(stringsAsFactors = F)
a=read.table('R_bilibili-master/R语言小作业/e2.txt')
colnames(a)='probe_id'
library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
tmp1=merge(ids,a,by='probe_id')
tmp2=ids[match(a$probe_id,ids$probe_id),]
比较tmp1tmp2：tmp1==tmp2
作业 3
找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量，并且绘制在 progres.-stable分组的boxplot图。
提示：
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex)
library(hgu95av2.db)
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex)
library(hgu95av2.db)
pd=pData(sCLLex)
ids=toTable(hgu95av2SYMBOL)
head(ids)
boxplot(exprSet['1939_at',] ~ pd$Disease)
boxplot(exprSet['1974_s_at',] ~ pd$Disease)
boxplot(exprSet['31618_at',] ~ pd$Disease)
作业 4
找到BRCA1基因在TCGA数据库的乳腺癌数据集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表达情况。
提示：使用http://www.cbioportal.org/index.do 定位数据集：http://www.cbioportal.org/datasets
rm(list = ls())
options(stringsAsFactors = F)
a=read.table('R_bilibili-master/R语言小作业/e4-plot.txt',
sep = '\t',
fill = T,
header = T)
colnames(a)=c('id','subtype','expression','mut')
dat=a
dat=dat[,-4]
library(ggstatsplot)
ggbetweenstats(data = dat, x = subtype, y = expression)
lastlibrary(ggplot2)
ggsave('plot-again-BRCA1-TCGA-BRCA-cbioportal.png')
tips:如果出现报错，Names must be unique.，那么将expression列名改为exp即可成功绘图。（可能是expression与某一函数同名会起冲突）
作业 5
找到TP53基因在TCGA数据库的乳腺癌数据集的表达量分组看其是否影响生存
提示使用：http://www.oncolnc.org/
rm(list = ls())
options(stringsAsFactors = F)
a=read.table('BRCA_7157_50_50.csv',sep = ',',fill = T,header = T)
dat=a
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)
ggsave('survival_TP53_in_BRCA_TCGA.png')
### 分割线
head(a)
b=read.table('e4-plot.txt',sep = '\t',fill = T,header = T)
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)
type='BRCA_LumB'
x=tmp[tmp$subtype==type,]
library(ggplot2)
library(survival)
library(survminer)
#table(x$Status)
x$Status=ifelse(x$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=x)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
table(tmp$subtype)
type='BRCA_Normal'
x=tmp[tmp$subtype==type,]
library(ggplot2)
library(survival)
library(survminer)
#table(x$Status)
x$Status=ifelse(x$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=x)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
table(tmp$subtype)
type='BRCA_Basal'
x=tmp[tmp$subtype==type,]
library(ggplot2)
library(survival)
library(survminer)
#table(x$Status)
x$Status=ifelse(x$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=x)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
table(tmp$subtype)
type='BRCA_Her2'
x=tmp[tmp$subtype==type,]
library(ggplot2)
library(survival)
library(survminer)
#table(x$Status)
x$Status=ifelse(x$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=x)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
table(tmp$subtype)
type='BRCA_LumA'
x=tmp[tmp$subtype==type,]
library(ggplot2)
library(survival)
library(survminer)
#table(x$Status)
x$Status=ifelse(x$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=x)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
作业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)
# 注意查看下载文件的大小，检查数据
f='GSE17215_eSet.Rdata'
library(GEOquery)
# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
gset <- getGEO('GSE17215', destdir=".",
AnnotGPL = F, ## 注释文件
getGPL = F) ## 平台文件
save(gset,file=f) ## 保存到本地
}
load('GSE17215_eSet.Rdata') ## 载入数据
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[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
下载数据集GSE24673的表达矩阵计算样本的相关性并且绘制热图，需要标记上样本分组信息
rm(list = ls())
options(stringsAsFactors = F)
# 注意查看下载文件的大小，检查数据
f='GSE24673_eSet.Rdata'
library(GEOquery)
# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
gset <- getGEO('GSE24673', destdir=".",
AnnotGPL = F, ## 注释文件
getGPL = 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)
pd=pData(a)
group_list=c('rbc','rbc','rbc',
'rbn','rbn','rbn',
'rbc','rbc','rbc',
'normal','normal')
dat[1:4,1:4]
M=cor(dat)
pheatmap::pheatmap(M)
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)
作业8
找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 对应的R的bioconductor注释包，并且安装它
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)
options()$repos
options()$BioC_mirror
作业9
下载数据集GSE42872的表达矩阵，并且分别挑选出 所有样本的(平均表达量/sd/mad/)最大的探针，并且找到它们对应的基因。
rm(list = ls())
options(stringsAsFactors = F)
# 注意查看下载文件的大小，检查数据
f='GSE42872_eSet.Rdata'
library(GEOquery)
# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
gset <- getGEO('GSE42872', destdir=".",
AnnotGPL = F, ## 注释文件
getGPL = 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/)最大的探针
boxplot(dat)
sort(apply(dat,1,mean),decreasing = T)[1]
sort(apply(dat,1,sd),decreasing = T)[1]
sort(apply(dat,1,mad),decreasing = T)[1]
作业10
下载数据集GSE42872的表达矩阵，并且根据分组使用limma做差异分析，得到差异结果矩阵
rm(list = ls())
options(stringsAsFactors = F)
# 注意查看下载文件的大小，检查数据
f='GSE42872_eSet.Rdata'
library(GEOquery)
# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
gset <- getGEO('GSE42872', destdir=".",
AnnotGPL = F, ## 注释文件
getGPL = 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/)最大的探针
boxplot(dat)
group_list=unlist(lapply(pd$title,function(x){
strsplit(x,' ')[[1]][4]
}))
exprSet=dat
exprSet[1:4,1:4]
# DEG 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
fit <- lmFit(exprSet,design)
##step2
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_notrend.results.csv",quote = F)
head(nrDEG)
