当看完几遍R语言书籍时，可以适当做些题目，而jimmy在生信技能树的这几道题目是很适合练习，每一行代码自己都要着去理解，在Rstudio中可以不断修改和试错，直到弄明白。http://www.bio-info-trainee.com/3750.html 其实这些题目有很多好的答案了。
但佛系笔记嘛，可能在做笔记前花了很多时间才搞懂为什么是这样写那样写，身边有很多小伙伴可以一起学习，终于做出来的感觉挺爽的。

第一题
请根据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)

先甩个生信技能树博客链接—了解常用数据库ID表示方式 http://www.biotrainee.com/thread-411-1-1.html

答1 > rm(list=ls()) #删除已存在变量，可删可不删，看个人习惯 > if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("org.Hs.eg.db")# Bioconductor 命令行下载安装org.Hs.eg.db 包 > suppressMessages(library(org.Hs.eg.db))#去除加载包时的警告，可用可不用 > ls("package:org.Hs.eg.db")#查看包的内容[1] "org.Hs.eg""org.Hs.eg.db" [3] "org.Hs.eg_dbconn""org.Hs.eg_dbfile" [5] "org.Hs.eg_dbInfo""org.Hs.eg_dbschema" [7] "org.Hs.egACCNUM""org.Hs.egACCNUM2EG" [9] "org.Hs.egALIAS2EG""org.Hs.egCHR" [11] "org.Hs.egCHRLENGTHS""org.Hs.egCHRLOC" [13] "org.Hs.egCHRLOCEND""org.Hs.egENSEMBL" [15] "org.Hs.egENSEMBL2EG""org.Hs.egENSEMBLPROT" [17] "org.Hs.egENSEMBLPROT2EG""org.Hs.egENSEMBLTRANS" [19] "org.Hs.egENSEMBLTRANS2EG" "org.Hs.egENZYME" [21] "org.Hs.egENZYME2EG""org.Hs.egGENENAME" [23] "org.Hs.egGO""org.Hs.egGO2ALLEGS" [25] "org.Hs.egGO2EG""org.Hs.egMAP" [27] "org.Hs.egMAP2EG""org.Hs.egMAPCOUNTS" [29] "org.Hs.egOMIM""org.Hs.egOMIM2EG" [31] "org.Hs.egORGANISM""org.Hs.egPATH" [33] "org.Hs.egPATH2EG""org.Hs.egPFAM" [35] "org.Hs.egPMID""org.Hs.egPMID2EG" [37] "org.Hs.egPROSITE""org.Hs.egREFSEQ" [39] "org.Hs.egREFSEQ2EG""org.Hs.egSYMBOL" [41] "org.Hs.egSYMBOL2EG""org.Hs.egUCSCKG" [43] "org.Hs.egUNIGENE""org.Hs.egUNIGENE2EG" [45] "org.Hs.egUNIPROT" > library("org.Hs.eg.db")#加载包 > keytypes(org.Hs.eg.db)# org.Hs.eg.db 包里面的数据有哪些类型，可以看到这题中我们需要的是ENSEMBL，和SYMBOL [1] "ACCNUM""ALIAS""ENSEMBL""ENSEMBLPROT""ENSEMBLTRANS" [6] "ENTREZID""ENZYME""EVIDENCE""EVIDENCEALL""GENENAME" [11] "GO""GOALL""IPI""MAP""OMIM" [16] "ONTOLOGY""ONTOLOGYALL""PATH""PFAM""PMID" [21] "PROSITE""REFSEQ""SYMBOL""UCSCKG""UNIGENE" [26] "UNIPROT"> x<-c("ENSG00000000003","ENSG00000000005","ENSG00000000419","ENSG00000000457","ENSG00000000460","ENSG00000000938") # 将这些ENSEMBL ID 赋值给x > cols <- c("SYMBOL", "GENENAME") > z <- select(org.Hs.eg.db,keys=x,columns=cols,keytype="ENSEMBL") # select 函数，将ensembl ID 对应到 gene name,赋值给z > View(z) #查看z

第二题
根据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)

答2 >if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("hgu95av2.db")#Bioconductor 下载安装hgu95av2.db 包， > suppressMessages(library(hgu133a.db))#去除加载包的警告 > library(hgu133a.db)#加载包 > keytypes(hgu133a.db)#查看包的的内容 [1] "ACCNUM""ALIAS""ENSEMBL""ENSEMBLPROT""ENSEMBLTRANS" [6] "ENTREZID""ENZYME""EVIDENCE""EVIDENCEALL""GENENAME" [11] "GO""GOALL""IPI""MAP""OMIM" [16] "ONTOLOGY""ONTOLOGYALL""PATH""PFAM""PMID" [21] "PROBEID""PROSITE""REFSEQ""SYMBOL""UCSCKG" [26] "UNIGENE""UNIPROT" > x <-c("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") > cols <- c("PROBEID","SYMBOL") > tmp<- select(hgu133a.db,keys=keys(hgu133a.db),columns = cols)#select 函数得到 > tmp2<-tmp[match(x,tmp$PROBEID),]#match 函数 > View(tmp2)

第三题:
找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量，并且绘制在 progres.-stable分组的boxplot图
提示：
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex)
library(hgu95av2.db)
想想如何通过 ggpubr 进行美化。

答3 > if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("CLL")#下载安装CLL包 > suppressMessages(library(CLL)) > library(CLL)#加载CLL包 > data(sCLLex) #调用sCLLex > library(hgu95av2.db)#加载hgu95av2.db包 > exprSet <- exprs(sCLLex) > pdata <- pData(sCLLex) > p2s <- toTable(hgu95av2SYMBOL) > p3 <- filter(p2s,symbol=='TP53')

第四题
找到BRCA1基因在TCGA数据库的乳腺癌数据集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表达情况
提示：使用http://www.cbioportal.org/index.do 定位数据集：http://www.cbioportal.org/datasets

可以在网站中直接查看，也可以用将数据下载后，R代码得到

第五题
找到TP53基因在TCGA数据库的乳腺癌数据集的表达量分组看其是否影响生存
提示使用：http://www.oncolnc.org/

找到官网，[http://www.oncolnc.org/](http://www.oncolnc.org/)

第六题
下载数据集GSE17215的表达矩阵并且提取下面的基因画热图
提示：根据基因名拿到探针ID，缩小表达矩阵绘制热图，没有检查到的基因直接忽略即可。

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

第七题
下载数据集GSE24673的表达矩阵计算样本的相关性并且绘制热图，需要标记上样本分组信息

rm(list=ls())#删除已有变量 options(stringsAsFactors = F) if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install('GEOquery') library(GEOquery)#下载安装GEOquery GSE <- "GSE24673"下载GSE24673 数据集 if(!file.exists(GSE)){ geo <- getGEO(GSE, destdir = '.', getGPL = F, AnnotGPL = F) save(geo, file = paste0(GSE,'.eSet.Rdata')) }#if 函数 得到并保存GSE,'.eSet.Rdata' 命名为pasteo文件，赋值给geo load(paste0(GSE,'.eSet.Rdata'))#加载pasteo expr <- exprs(geo[[1]]) #获得geo的exprimentData数据，并赋值给expr dim(expr) 查看expr有几行几列 expr[1:4,1:4]#查看expr中一至四行，一至四列的内容 pdata <- pData(geo[[1]])#提取phenoData的数据，并赋值给pdata # 自己根据pdata第八列做一个分组信息矩阵 grp <- c('rbc','rbc','rbc', 'rbn','rbn','rbn', 'rbc','rbc','rbc', 'normal','normal')grp_df <- data.frame(group=grp) #将grp 得到一个数据框 rownames(grp_df) <- colnames(expr) #列名等于expr中的名字 new_cor <- cor(expr) pheatmap::pheatmap(new_cor, annotation_col = grp_df) #热图

文章图片

第八题
找到 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/")#Bioconductor 镜像 options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")) #清华镜像 BiocManager::install("hugene10sttranscriptcluster",ask = F,update =T) #R包是"hugene10sttranscriptcluster"，更新到最新版本 options()$repos options()$BioC_mirror #参考芯片与注释包的关系 [https://www.jianshu.com/p/f6906ba703a0](https://www.jianshu.com/p/f6906ba703a0)

第九题
下载数据集GSE42872的表达矩阵，并且分别挑选出 所有样本的(平均表达量/sd/mad/)最大的探针，并且找到它们对应的基因。

library(GEOquery)#加载GEOquery包 GSE<-GSE42872 if(!file.exists(GSE)){ geo <- getGEO(GSE, destdir = '.', getGPL = F, AnnotGPL = F) save(geo, file = paste0(GSE,'.eSet.Rdata')) } exprSet = exprs(gset[[1]])#获得表达矩阵 # exprSet <- exprSet[apply(exprSet, 1, function(x) sum(x>1) > 5),] # exprSet <- log(edgeR::cpm(exprSet) + 1) # 处理中平均值最大的探针 max_mean <- names(sort(apply(exprSet,1, mean), decreasing = T))[1] # 处理中sd最大的探针 max_sd <- names(sort(apply(exprSet,1, sd), decreasing = T))[1] # 处理中mad最大的探针 max_mad <- names(sort(apply(exprSet,1, mad), decreasing = T))[1] # 查询探针ID对应基因的symbol library(hugene10sttranscriptcluster.db) select(hugene10sttranscriptcluster.db, keys = c(max_mean, max_sd, max_mad), columns = c("SYMBOL"), keytype = "PROBEID")

【生信技能树-R语言小作业（中级）】第十题
下载数据集GSE42872的表达矩阵，并且根据分组使用limma做差异分析，得到差异结果矩阵

gset = getGEO('GSE42872', GSEMatrix = TRUE, AnnotGPL = FALSE, getGPL= FALSE) s <- toTable(hugene10sttranscriptclusterSYMBOL) exprSet = exprs(gset[[1]]) group_list <- as.character(c(rep("Control",3),rep("Vemurafenib",3))) dim(exprSet) # 过滤【只保留和注释文件探针id相同的探针】 efilt <- exprSet[rownames(exprSet)%in%s$probe_id,] dim(efilt) ## 整合1【目的：保证一个基因对应一个探针；如果基因和探针一一对应很好说，但如果一个基因对应多个探针：每个探针取一行的均值-》对应同一基因的探针取表达量最大的探针-》按照基因名给他们建索引，因为是按照基因来过滤探针（不用s$probe_id构建索引的原因是，看清楚我们的目的是让注释包的一个基因对应我们自己表达矩阵的一个探针。如果用s$probe_id那么结果就成了让注释包的一个探针对应我们自己表达矩阵的一个探针，当然这样也运行不成功，因为自己表达矩阵的探针过滤后的数量和注释包的探针数量不相等，这样没法一一对应。但基因名数量是不变的，什么是索引？以不变应万变的就是索引）】 maxp = by(efilt,s$symbol,function(x) rownames(x)[which.max(rowMeans(x))]) uniprobes = as.character(maxp) efilt=efilt[rownames(efilt)%in%uniprobes,] ## 整合2【目的：将我们表达矩阵的行名换成刚才一对一的基因名，并且match这个函数保证了表达矩阵和注释包的顺序是一致的】 rownames(efilt)=s[match(rownames(efilt),s$probe_id),2]# 差异分析 suppressMessages(library(limma)) #limma需要三个矩阵：表达矩阵（efilt）、分组矩阵（design）、比较矩阵（contrast） #先做一个分组矩阵～design，说明MAO是哪几个样本，MNO又是哪几个，其中1代表“是” design <- model.matrix(~0+factor(group_list)) colnames(design) <- levels(factor(group_list)) rownames(design) <- colnames(efilt) design #再做一个比较矩阵【一般是case比control】 contrast<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design) contrast DEG <- function(efilt,design,contrast){ ##step1 fit <- lmFit(efilt,design) ##step2 fit2 <- contrasts.fit(fit, contrast) fit2 <- eBayes(fit2) ##step3 mtx = topTable(fit2, coef=1, n=Inf) deg_mtx = na.omit(mtx) return(deg_mtx) } DEG_mtx <- DEG(efilt,design,contrast) #得到全部的差异基因矩阵

推荐阅读

• 

  吐奶|宝宝吐奶怎么办？先别急着送医院，看看这些干货
• 

  cookie底层分析,广告cookie 第三方分析
• 

  新居落成祝福语 盘点新居落成祝福语
• 

  如何登录云端服务器？ 云端服务器怎么登陆
• 

  桔子可以放冰箱保鲜吗？如何正确保存桔子以延长其新鲜度
• 

  涂防晒霜需要卸妆吗 不化妆只涂防晒霜需要卸妆吗
• 

  网盘解析,rar网盘解析
• 

  第一次使用充电宝要全部用尽吗 首次使用充电宝要全部用尽吗
• 

  用于捐赠物品的进项税可以抵扣吗 哪些捐赠允许税前抵扣，税法业务招待费扣除标准
• 

  表格中分析出数据的奇偶,excel表格怎么分析
• 

  大麦若叶青汁过期还能喝吗
• 

  家庭常用锅炉的清洗方法，家用燃气锅炉清洗
• 

  中国国旗简笔画 中国国旗简笔画怎么画
• 

  范成大简介及作品 范成大的作品有哪些
• 

  英文辞职报告模板
• 

  博世壁挂炉报ce怎么处理 郑州市博世壁挂炉总经销处理方法
• 

  安卓手机重装系统软件,手机出现病例怎么办?教你如何自行刷机
• 

  山楂树栽种注意事项 山楂树种植注意事项
• 

  手机投影仪可以给手机充电吗？荣耀手机怎么连接投影仪？
• 

  ug自动判断创建基准面的方法

• CET4听力微技能一
• 【生信技能树】R语言练习题|【生信技能树】R语言练习题 - 中级
• java中如何实现重建二叉树
• 工具|后天就是七夕节，你准备好了吗（送上几个七夕代码，展示你技能的时候到了！）
• 种树郭橐驼传（文言句式+古今异义+词类活用+通假字）
• 白杨树
• 08黑龙江迟淑荣弯柳树网络学院第五期学习赵宗瑞老师主讲的（传统文化与身心健康）教育体系心得体会
• 新媒体时代，你需要掌握的必备技能
• ［原创］能见沂山一棵树，胜读十年无用书！
• 涵养字外功