生信人的20个R语言习题
请先看：生信人的linux考试
- 安装一些R包：
数据包： ALL, CLL, pasilla, airway
软件包：limma，DESeq2，clusterProfiler
工具包：reshape2
绘图包：ggplot2
不同领域的R包使用频率不一样，在生物信息学领域，尤其需要掌握bioconductor系列包。
1 安装R包
# rm(list = ls())
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/"))
options()$repos
options()$BioC_mirror
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(c(" ALL","CLL","pasilla","airway" ),ask = F,update = F) #数据包,ask=F表示不询问中间过程和不升级
BiocManager::install(c("limma","DESeq2","clusterProfiler" ),ask = F,update = F) # 软件包
BiocManager::install(c("genefu","org.Hs.eg.db","hgu133plus2.db" ),ask = F,update = F)
BiocManager::install(c("reshape2"),ask = F,update = F) # 工具包
BiocManager::install(c("ggplot2"),ask = F,update = F) # 绘图包
#另外可以检测某个包是否存在，只有不存在时才会安装 if (! require ('CLL')){ options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
BiocInstaller::biocLite('CLL',ask = F, suppressUpdates = T) }
- 了解ExpressionSet对象，比如
CLL包里面就有data(sCLLex)，找到它包含的元素，提取其表达矩阵(使用exprs函数)，查看其大小
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
exprSet=exprs(sCLLex)
##sCLLex是依赖于CLL这个package的一个对象
samples=sampleNames(sCLLex)
pdata=pData(sCLLex)
group_list=as.character(pdata[,2])
dim(exprSet)
# [1] 12625 22
exprSet[1:5,1:5]
# CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL
# 1000_at 5.743132 6.219412 5.523328 5.340477 5.229904
# 1001_at 2.285143 2.291229 2.287986 2.295313 2.662170
# 1002_f_at 3.309294 3.318466 3.354423 3.327130 3.365113
# 1003_s_at 1.085264 1.117288 1.084010 1.103217 1.074243
# 1004_at 7.544884 7.671801 7.474025 7.152482 6.902932
- 了解 str,head,help函数，作用于 第二步提取到的表达矩阵
str(exprSet)
# str: Compactly display the internal structure of an R object, a diagnostic function and an alternative to summary (and to some extent, dput).
head(exprSet)
- 安装并了解
hgu95av2.db包,看看ls("package:hgu95av2.db")后 显示的那些变量
BiocManager::install("hgu95av2.db")
library("hgu95av2.db")
ls("package:hgu95av2.db")
[1] "hgu95av2" "hgu95av2.db" "hgu95av2_dbconn" "hgu95av2_dbfile"
[5] "hgu95av2_dbInfo" "hgu95av2_dbschema" "hgu95av2ACCNUM" "hgu95av2ALIAS2PROBE"
[9] "hgu95av2CHR" "hgu95av2CHRLENGTHS" "hgu95av2CHRLOC" "hgu95av2CHRLOCEND"
[13] "hgu95av2ENSEMBL" "hgu95av2ENSEMBL2PROBE" "hgu95av2ENTREZID" "hgu95av2ENZYME"
[17] "hgu95av2ENZYME2PROBE" "hgu95av2GENENAME" "hgu95av2GO" "hgu95av2GO2ALLPROBES"
[21] "hgu95av2GO2PROBE" "hgu95av2MAP" "hgu95av2MAPCOUNTS" "hgu95av2OMIM"
[25] "hgu95av2ORGANISM" "hgu95av2ORGPKG" "hgu95av2PATH" "hgu95av2PATH2PROBE"
[29] "hgu95av2PFAM" "hgu95av2PMID" "hgu95av2PMID2PROBE" "hgu95av2PROSITE"
[33] "hgu95av2REFSEQ" "hgu95av2SYMBOL" "hgu95av2UNIGENE" "hgu95av2UNIPROT"
- 理解
head(toTable(hgu95av2SYMBOL))的用法，找到 TP53 基因对应的探针ID
- 理解探针与基因的对应关系，总共多少个基因，基因最多对应多少个探针，是哪些基因，是不是因为这些基因很长，所以在其上面设计多个探针呢？
- 第二步提取到的表达矩阵是12625个探针在22个样本的表达量矩阵，找到那些不在
hgu95av2.db包收录的对应着SYMBOL的探针。- 提示：有1165个探针是没有对应基因名字的。
- 过滤表达矩阵，删除那1165个没有对应基因名字的探针。
- 整合表达矩阵，多个探针对应一个基因的情况下，只保留在所有样本里面平均表达量最大的那个探针。
- 提示，理解 tapply,by,aggregate,split 函数 , 首先对每个基因找到最大表达量的探针。
- 然后根据得到探针去过滤原始表达矩阵
- 把过滤后的表达矩阵更改行名为基因的symbol，因为这个时候探针和基因是一对一关系了。
- 对第10步得到的表达矩阵进行探索，先画第一个样本的所有基因的表达量的boxplot,hist,density ， 然后画所有样本的 这些图
- 参考：http://bio-info-trainee.com/tmp/basic_visualization_for_expression_matrix.html
- 理解ggplot2的绘图语法，数据和图形元素的映射关系
- 理解统计学指标mean,median,max,min,sd,var,mad并计算出每个基因在所有样本的这些统计学指标，最后按照mad值排序，取top 50 mad值的基因，得到列表。
- 注意：这个题目出的并不合规，请仔细看。
- 根据第12步骤得到top 50 mad值的基因列表来取表达矩阵的子集，并且热图可视化子表达矩阵。试试看其它5种热图的包的不同效果。
- 取不同统计学指标mean,median,max,mean,sd,var,mad的各top50基因列表，使用UpSetR包来看他们之间的overlap情况。
- 在第二步的基础上面提取
CLL包里面的data(sCLLex)数据对象的样本的表型数据。 - 对所有样本的表达矩阵进行聚类并且绘图，然后添加样本的临床表型数据信息(更改样本名)
- 对所有样本的表达矩阵进行PCA分析并且绘图，同样要添加表型信息。
- 根据表达矩阵及样本分组信息进行批量T检验，得到检验结果表格
- 使用limma包对表达矩阵及样本分组信息进行差异分析，得到差异分析表格，重点看logFC和P值，画个火山图(就是logFC和-log10(P值)的散点图。)。
- 对T检验结果的P值和limma包差异分析的P值画散点图，看看哪些基因相差很大。
