数据-R语言-图表-决策-Linux-PythonR语言做生信R小推车

R语言练习题10道，有答案代码，还有视频

2019-03-07  本文已影响29人  土豆学生信

来源：http://www.bio-info-trainee.com/3750.html
答案：https://github.com/jmzeng1314/R_bilibili
视频：【生信技能树】生信人应该这样学R语言 https://www.bilibili.com/video/av25643438/?p=13
备注：最好自己下载答案，因为Jimmy老师已经把需要的数据、代码打包了，可以完美演示（太贴心了）！前五题主要针对一些基本操作，后五题主要是针对GEO数据的一些操作。

作业 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)

答案如下：
重点：org.Hs.eg.db包、lapply函数、merge函数

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]

作业 2

根据R包hgu133a.db找到下面探针对应的基因名(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)

答案如下

重点：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),]

作业 3

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

提示：

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

想想如何通过 ggpubr 进行美化。
答案如下：

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")
#发现自己写的代码又丑又长

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

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

作业 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)
#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

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

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

作业8

找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 对应的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获取芯片探针与基因的对应关系三部曲-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]
image.png

作业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)

注意：代码中"progres.-stable"这一行会报错，奇怪的是Jimmy老师视频中没有报错：

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

原因是没有这个分组名称，Jimmy老师应该想表达的是"Control-Vemurafenib"这样一个比较。

Attention！

生信基础知识大全系列：生信基础知识100讲

史上最强的生信自学环境准备课来啦！！

7次改版，11节课程，14K的讲稿，30个夜晚打磨，100页PPT的课程。

如果需要组装自己的服务器：代办生物信息学服务器

如果需要帮忙下载海外数据(GEO/TCGA/GTEx等等)，点我？

如果需要线下辅导及培训，看招学徒** (培训)**

如果需要个人电脑：个人计算机推荐

如果需要置办生物信息学书籍，看：生信人必备书单

如果需要实习岗位：实习职位发布

如果需要售后：点我

如果需要入门资料大全：点我

上一篇下一篇

猜你喜欢

热点阅读