這是生信技能樹論壇R語言的中級(jí)測(cè)試題芳绩，其中大部分是GEO數(shù)據(jù)挖掘和TCGA數(shù)據(jù)庫(kù)的一些操作路召，雖然我目前用不到偷线，但是卻可以通過查看大佬的代碼來提高自己的認(rèn)識(shí)水平柱嫌，也可以對(duì)其他領(lǐng)域有個(gè)粗淺的認(rèn)識(shí)历帚。

1.根據(jù)R包org.Hs.eg.db找到下面ensembl 基因ID 對(duì)應(yīng)的基因名(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.根據(jù)R包hgu133a.db找到下面探針對(duì)應(yīng)的基因名(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(),返回的是前面的向量在后面的坐標(biāo)
c <- ids[match(a$V1,ids$probe_id),]
b==c #判斷bc是否相等

3.找到R包CLL內(nèi)置的數(shù)據(jù)集的表達(dá)矩陣?yán)锩娴腡P53基因的表達(dá)量挽牢，并且繪制在 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) #查看分組信息
#五個(gè)以下不用寫循環(huán)
boxplot(exprSet["1939_at",]~pd$Disease)
boxplot(exprSet["1974_s_at",]~pd$Disease)
boxplot(exprSet["31618_at",]~pd$Disease)

Figure3.png

4.找到BRCA1基因在TCGA數(shù)據(jù)庫(kù)的乳腺癌數(shù)據(jù)集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表達(dá)情況

提示：使用http://www.cbioportal.org/index.do 定位數(shù)據(jù)集：http://www.cbioportal.org/datasets

rm(list = ls())
options(stringsAsFactors = F)
#fill=T，不管如何都會(huì)讀取文件
a=read.table('e4-plot.txt',sep = '\t',fill = T,header = T)

colnames(a)=c('id','subtype','expression','mut')
BiocManager::install("ggstatsplot")
#畫boxplot摊求，帶統(tǒng)計(jì)值的
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數(shù)據(jù)庫(kù)的乳腺癌數(shù)據(jù)集的表達(dá)量分組看其是否影響生存

提示使用：http://www.oncolnc.org/.整體不顯著的時(shí)候要看其亞型是否顯著。

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值變?yōu)?.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.下載數(shù)據(jù)集GSE17215的表達(dá)矩陣并且提取下面的基因畫熱圖
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
提示：根據(jù)基因名拿到探針I(yè)D室叉，縮小表達(dá)矩陣?yán)L制熱圖睹栖，沒有檢查到的基因直接忽略即可。

rm(list = ls())  ## 魔幻操作茧痕，一鍵清空~
options(stringsAsFactors = F)
# 注意查看下載文件的大小野来，檢查數(shù)據(jù) 
f='GSE17215_eSet.Rdata'

library(GEOquery)
# 這個(gè)包需要注意兩個(gè)配置，一般來說自動(dòng)化的配置是足夠的凿渊。
#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)       ## 平臺(tái)文件
  save(gset,file=f)   ## 保存到本地
}
load('GSE17215_eSet.Rdata')  ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺(tái)梁只，所以下載到的是一個(gè)含有一個(gè)元素的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后進(jìn)行提取
ng=strsplit(ng,' ')[[1]]
#%in%判斷在不在缚柳，match也可以
table(ng %in%  rownames(dat))
#取返回值為T的
ng=ng[ng %in%  rownames(dat)]
dat=dat[ng,]
#因?yàn)闊釄D的值比較大埃脏，所以需要log一下
dat=log2(dat)
#scale=‘row’對(duì)row進(jìn)行歸一化
pheatmap::pheatmap(dat,scale = 'row')

7.下載數(shù)據(jù)集GSE24673的表達(dá)矩陣計(jì)算樣本的相關(guān)性并且繪制熱圖，需要標(biāo)記上樣本分組信息.【這里展示了一個(gè)操作秋忙，ctrl+F彩掐，對(duì)選定植進(jìn)行一鍵替換】

rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
# 注意查看下載文件的大小灰追，檢查數(shù)據(jù) 
f='GSE24673_eSet.Rdata'

library(GEOquery)
# 這個(gè)包需要注意兩個(gè)配置堵幽，一般來說自動(dòng)化的配置是足夠的。
#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)       ## 平臺(tái)文件
  save(gset,file=f)   ## 保存到本地
}
load('GSE24673_eSet.Rdata')  ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺(tái)弹澎，所以下載到的是一個(gè)含有一個(gè)元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)
#查看分組信息
pd=pData(a)
#因?yàn)檫@里的分組信息亂朴下，干脆就自己創(chuàng)建list
group_list=c('rbc','rbc','rbc',
             'rbn','rbn','rbn',
             'rbc','rbc','rbc',
             'normal','normal')
#做相關(guān)性之前要查看部分?jǐn)?shù)據(jù)的相關(guān)性
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 對(duì)應(yīng)的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("請(qǐng)輸入自己找到的R包",ask = F,update = F)
options()$repos
options()$BioC_mirror

9.下載數(shù)據(jù)集GSE42872的表達(dá)矩陣殴胧，并且分別挑選出 所有樣本的(平均表達(dá)量/sd/mad/)最大的探針，并且找到它們對(duì)應(yīng)的基因佩迟。

rm(list = ls())  ## 魔幻操作团滥，一鍵清空~
options(stringsAsFactors = F)
# 注意查看下載文件的大小，檢查數(shù)據(jù) 
f='GSE42872_eSet.Rdata'

library(GEOquery)
# 這個(gè)包需要注意兩個(gè)配置报强，一般來說自動(dòng)化的配置是足夠的灸姊。
#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)       ## 平臺(tái)文件
  save(gset,file=f)   ## 保存到本地
}
load('GSE42872_eSet.Rdata')  ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺(tái)，所以下載到的是一個(gè)含有一個(gè)元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
# (平均表達(dá)量/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.下載數(shù)據(jù)集GSE42872的表達(dá)矩陣秉溉，并且根據(jù)分組使用limma做差異分析力惯，得到差異結(jié)果矩陣

rm(list = ls())  ## 魔幻操作碗誉，一鍵清空~
options(stringsAsFactors = F)
# 注意查看下載文件的大小，檢查數(shù)據(jù) 
f='GSE42872_eSet.Rdata'

library(GEOquery)
# 這個(gè)包需要注意兩個(gè)配置父晶，一般來說自動(dòng)化的配置是足夠的诗充。
#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)       ## 平臺(tái)文件
  save(gset,file=f)   ## 保存到本地
}
load('GSE42872_eSet.Rdata')  ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺(tái)，所以下載到的是一個(gè)含有一個(gè)元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
# (平均表達(dá)量/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 
##這個(gè)矩陣聲明诱建，我們要把progres.組跟stable進(jìn)行差異分析比較
##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.我們不生產(chǎn)代碼，我們只是代碼的板運(yùn)工俺猿。
2.理解R包原理后會(huì)修改一些較復(fù)雜圖的參數(shù)茎匠。
3.apply,%in%,match(),字符串操作(base和stringr，拆分用的較多)押袍，sort诵冒，ifelse。
4.普通boxplot谊惭，帶統(tǒng)計(jì)值的boxplot汽馋，生存分析。
5.與網(wǎng)頁(yè)小工具的交替使用圈盔，可提高效率豹芯。

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

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