这是生信技能树论坛R语言的中级测试题，其中大部分是GEO数据挖掘和TCGA数据库的一些操作，虽然我目前用不到，但是却可以通过查看大佬的代码来提高自己的认识水平，也可以对其他领域有个粗浅的认识。

1.根据R包org.Hs.eg.db找到下面ensembl 基因ID 对应的基因名(symbol)
ENSG00000000003.13
ENSG00000000005.5
ENSG00000000419.11
ENSG00000000457.12
ENSG00000000460.15
ENSG00000000938.11

BiocManager::install("org.Hs.eg.db")
library(stringr)
library(org.Hs.eg.db)
options(stringsAsFactors = F)
a <- read.csv("s2-1.txt",header = F)
g2s=toTable(org.Hs.egSYMBOL)
g2e=toTable(org.Hs.egENSEMBL)
class(a)
b <- unlist(lapply(a$V1,function(x){
  str_split(x,"[.]")[[1]][1]}))
b <- as.data.frame(b)
c <- merge(b,g2e,by.x="b",by.y="ensembl_id")
d <- merge(c,g2s,by.x="gene_id",by.y="gene_id")
head(d)

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

BiocManager::install("hgu133a.db")
library(hgu133a.db)
library(stringr)
ids=toTable(hgu133aSYMBOL)
head(ids)
options(stringsAsFactors = F)
a <- read.csv("s2-2.txt",sep = "\t",header = F)
#方法1，使用merge()
b <- merge(a,ids,by.x="V1",by.y="probe_id")
head(b)
#方法二，使用match(),返回的是前面的向量在后面的坐标
c <- ids[match(a$V1,ids$probe_id),]
b==c #判断bc是否相等

3.找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量，并且绘制在 progres.-stable分组的boxplot图

options(stringsAsFactors = F)
rm(list=ls())
BiocManager::install("CLL")
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex) 
BiocManager::install("hgu95av2.db")
library(hgu95av2.db)
ids=toTable(hgu95av2SYMBOL)

pd <- pData(sCLLex) #查看分组信息
#五个以下不用写循环
boxplot(exprSet["1939_at",]~pd$Disease)
boxplot(exprSet["1974_s_at",]~pd$Disease)
boxplot(exprSet["31618_at",]~pd$Disease)

Figure3.png

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)
#fill=T，不管如何都会读取文件
a=read.table('e4-plot.txt',sep = '\t',fill = T,header = T)

colnames(a)=c('id','subtype','expression','mut')
BiocManager::install("ggstatsplot")
#画boxplot，带统计值的
library(ggstatsplot)
ggbetweenstats(data =a, x = subtype,  y = expression)
#b与a外观一致，不知a其他两列的作用是什么
b <- data.frame(subtype=a$subtype,expression=a$expression)
ggbetweenstats(data =b, x = subtype,  y = expression)
#保存
library(ggplot2)
ggsave('plot-again-BRCA1-TCGA-BRCA-cbioportal.png')

Figure4.png

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)
#将status列dead值变为1.alive值为0
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')
#查看亚型之间是否显著
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)  

#ifelse()的使用
ifelse(test, yes, no)

Figure5.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)
# 注意查看下载文件的大小，检查数据 
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'
#拆分成合适的character后进行提取
ng=strsplit(ng,' ')[[1]]
#%in%判断在不在，match也可以
table(ng %in%  rownames(dat))
#取返回值为T的
ng=ng[ng %in%  rownames(dat)]
dat=dat[ng,]
#因为热图的值比较大，所以需要log一下
dat=log2(dat)
#scale=‘row’对row进行归一化
pheatmap::pheatmap(dat,scale = 'row')

7.下载数据集GSE24673的表达矩阵计算样本的相关性并且绘制热图，需要标记上样本分组信息.【这里展示了一个操作，ctrl+F，对选定植进行一键替换】

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)
#因为这里的分组信息乱，干脆就自己创建list
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注释包，并且安装它！
注意：所有的注释R包后面都要加.db

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("请输入自己找到的R包",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 limma接受的是log后的矩阵
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)

启示
1.我们不生产代码，我们只是代码的板运工。
2.理解R包原理后会修改一些较复杂图的参数。
3.apply,%in%,match(),字符串操作(base和stringr，拆分用的较多)，sort，ifelse。
4.普通boxplot，带统计值的boxplot，生存分析。
5.与网页小工具的交替使用，可提高效率。

原文链接：http://www.bio-info-trainee.com/3750.html
B站视频解答：https://www.bilibili.com/video/av25643438?p=13

Github原始代码：https://github.com/jmzeng1314/R_bilibili