R語言練習(xí)題10道，有答案代碼，還有視頻

來源：http://www.bio-info-trainee.com/3750.html
答案：https://github.com/jmzeng1314/R_bilibili
視頻：【生信技能樹】生信人應(yīng)該這樣學(xué)R語言 https://www.bilibili.com/video/av25643438/?p=13
備注：最好自己下載答案，因?yàn)镴immy老師已經(jīng)把需要的數(shù)據(jù)、代碼打包了，可以完美演示（太貼心了）！前五題主要針對一些基本操作，后五題主要是針對GEO數(shù)據(jù)的一些操作。

作業(yè) 1

根據(jù)R包org.Hs.eg.db找到下面ensembl 基因ID 對應(yīng)的基因名(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)

答案如下：
重點(diǎn)：org.Hs.eg.db包、lapply函數(shù)、merge函數(shù)

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('e1.txt')
head(a)
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)
head(g2e)
library(stringr)
a$ensembl_id=unlist(lapply(a$V1,function(x){
  strsplit(as.character(x),'[.]')[[1]][1]
})
)
tmp1=merge(a,g2e,by='ensembl_id')
tmp2=merge(tmp1,g2s,by='gene_id')
#或者采用下列方法
data <- merge(g2s, g2e, by='gene_id');head(data)
tmp3=merge(a,data,by = c(" V1"= "ensembl_id"))[,-2]

作業(yè) 2

根據(jù)R包hgu133a.db找到下面探針對應(yīng)的基因名(symbol)

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

答案如下

重點(diǎn)：hgu133a.db包的用法

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('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),]

作業(yè) 3

找到R包CLL內(nèi)置的數(shù)據(jù)集的表達(dá)矩陣?yán)锩娴腡P53基因的表達(dá)量，并且繪制在 progres.-stable分組的boxplot圖

提示：

suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex)   
library(hgu95av2.db)

想想如何通過 ggpubr 進(jìn)行美化。
答案如下：

rm(list = ls())
options(stringsAsFactors = F)
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex)
pd=pData(sCLLex)
library(hgu95av2.db) 
ids=toTable(hgu95av2SYMBOL)
head(ids)
boxplot(exprSet['1939_at',] ~ pd$Disease) ## sig
boxplot(exprSet['1974_s_at',] ~ pd$Disease)
boxplot(exprSet['31618_at',] ~ pd$Disease)
#ggpubr
d<-cbind(as.data.frame(exprSet['1939_at',]),as.data.frame(pd$Disease))
head(d)
library(ggpubr)
colnames(d)<-c("e","f")
p<-ggboxplot(d, x=("f"), y =("e") ,
          color = "f",
          palette = "jco",
          add = "jitter")
opar<-par(no.readonly=T)
#  Add p-value
p + stat_compare_means()
# Change method
p + stat_compare_means(method = "t.test")
#發(fā)現(xiàn)自己寫的代碼又丑又長

可視化參考：http://www.sthda.com/english/articles/32-r-graphics-essentials/

image.png

作業(yè) 4

找到BRCA1基因在TCGA數(shù)據(jù)庫的乳腺癌數(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)
a=read.table('e4-plot.txt',sep = '\t',fill = T,header = T)
View(a)
colnames(a)=c('id','subtype','expression','mut')
dat=a
library(ggstatsplot)
ggbetweenstats(data =dat, x = subtype,  y = expression)
library(ggplot2)
ggsave('plot-again-BRCA1-TCGA-BRCA-cbioportal.png')
image.png

作業(yè) 5

找到TP53基因在TCGA數(shù)據(jù)庫的乳腺癌數(shù)據(jù)集的表達(dá)量分組看其是否影響生存

提示使用：http://www.oncolnc.org/
答案如下：

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('BRCA_7157_50_50.csv',sep = ',',fill = T,header = T)
#View(a)
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')
image.png

作業(yè)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è)配置，一般來說自動化的配置是足夠的。
#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')  ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺，所以下載到的是一個(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'
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

作業(yè)7

下載數(shù)據(jù)集GSE24673的表達(dá)矩陣計(jì)算樣本的相關(guān)性并且繪制熱圖，需要標(biāo)記上樣本分組信息
答案如下：

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

library(GEOquery)
# 這個(gè)包需要注意兩個(gè)配置，一般來說自動化的配置是足夠的。
#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')  ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺，所以下載到的是一個(gè)含有一個(gè)元素的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)
image.png

作業(yè)8

找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 對應(yīng)的R的bioconductor注釋包，并且安裝它！

options()$reposoptions()$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()$reposoptions()$BioC_mirror

答案如下：可以參考
用R獲取芯片探針與基因的對應(yīng)關(guān)系三部曲-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

作業(yè)9

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

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

library(GEOquery)
# 這個(gè)包需要注意兩個(gè)配置，一般來說自動化的配置是足夠的。
#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')  ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺，所以下載到的是一個(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]
image.png

作業(yè)10

下載數(shù)據(jù)集GSE42872的表達(dá)矩陣，并且根據(jù)分組使用limma做差異分析，得到差異結(jié)果矩陣
答案如下：

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

library(GEOquery)
# 這個(gè)包需要注意兩個(gè)配置，一般來說自動化的配置是足夠的。
#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')  ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺，所以下載到的是一個(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 
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)

注意：代碼中"progres.-stable"這一行會報(bào)錯(cuò)，奇怪的是Jimmy老師視頻中沒有報(bào)錯(cuò)：

Error in eval(ej, envir = levelsenv) : object 'progres.' not found

原因是沒有這個(gè)分組名稱，Jimmy老師應(yīng)該想表達(dá)的是"Control-Vemurafenib"這樣一個(gè)比較。

Attention！

生信基礎(chǔ)知識大全系列：生信基礎(chǔ)知識100講

史上最強(qiáng)的生信自學(xué)環(huán)境準(zhǔn)備課來啦！！

7次改版，11節(jié)課程，14K的講稿，30個(gè)夜晚打磨，100頁P(yáng)PT的課程。

如果需要組裝自己的服務(wù)器：代辦生物信息學(xué)服務(wù)器

如果需要幫忙下載海外數(shù)據(jù)(GEO/TCGA/GTEx等等)，點(diǎn)我？

如果需要線下輔導(dǎo)及培訓(xùn)，看招學(xué)徒** (培訓(xùn))**

如果需要個(gè)人電腦：個(gè)人計(jì)算機(jī)推薦

如果需要置辦生物信息學(xué)書籍，看：生信人必備書單

如果需要實(shí)習(xí)崗位：實(shí)習(xí)職位發(fā)布

如果需要售后：點(diǎn)我

如果需要入門資料大全：點(diǎn)我

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容