做了一段时间的学徒，走的是周末班的路子，作为一个纯湿实验转生信分析的人，一开始是看不懂曾老师的视频的，现在R语言linux这些有了近一个月的练习了，上周在各种前辈的推文协助下完成了linux20题，这周就把20个R语言习题做一下吧。
参考：
生信人的20个R语言习题
20个R语言习题
R语言20练习题【完整版】
生信人的20个R语言习题及其答案-土豆学习笔记
Day14-生信人的20个R语言习题第11-20题
1.安装一些R包：
数据包： ALL, CLL, pasilla, airway
软件包：limma，DESeq2，clusterProfiler
工具包：reshape2
绘图包：ggplot2
不同领域的R包使用频率不一样，在生物信息学领域，尤其需要掌握bioconductor系列包。
- 解析：已经提示了是
bioconductor系列包,那么首先是了解一下什么是bioconductor，他的包又需要怎么安装？
首先在Bing浏览器输入bioconductor关键词，跳出来第一个就是：
1
进到官网再查看一下相关信息
2
简洁明了：在搜索框输入需要下载的包，查看如何下载
CLL
出来的东西特别多，不懂的时候可能就比较慌，没关系，慢慢探索就好；点击第一个，这个检索对象是个数据包，一般比较大，搜索的人估计也不多，网速不好打开还有点费劲，但是没关系，出来之后就能给我们提供很多信息
CLL
从图中可以看到，安装的代码只有一行BiocManager::install("CLL")但是如果是初次安装需要先装安装这个包的包，就是上面的代码if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")(如果已经安装了就不需要，这里面有个判断语句，如果需要就安装，如果有就不装，还是很人性化的) - 这个网站上的包，安装的方式也都是一样的，前提是你的rstudio版本要对应适当的代码，因为生信发展太快，他会不停的更新，我们学习者也要与时俱进。
了解了以上知识之后，就能很快完成第一题了：
rm(list = ls())
options(stringsAsFactors = F)
# 数据包： ALL, CLL, pasilla, airway
# 软件包：limma，DESeq2，clusterProfiler
# 工具包：reshape2
# 绘图包：ggplot2
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
if (!requireNamespace("ALL", quietly = TRUE))
BiocManager::install("ALL")
if (!requireNamespace("CLL", quietly = TRUE))
BiocManager::install("CLL")
if (!requireNamespace("pasilla", quietly = TRUE))
BiocManager::install("pasilla")
if (!requireNamespace("airway", quietly = TRUE))
BiocManager::install("airway")
if (!requireNamespace("limma", quietly = TRUE))
BiocManager::install("limma")
if (!requireNamespace("DESeq2", quietly = TRUE))
BiocManager::install("DESeq2")
if (!requireNamespace("clusterProfiler", quietly = TRUE))
BiocManager::install("clusterProfiler")
if (!requireNamespace("reshape2", quietly = TRUE))
BiocManager::install("reshape2")
if (!requireNamespace("ggplot2", quietly = TRUE))
BiocManager::install("ggplot2")
- 因为下载安装包需要费时间，有的时候还需要更新，最好是一个一个安装，当然熟练了以后可以一起批量安装，然后不显示更新信息及选择不更新自动处理，初学者可以跟着曾老师B站的视频反复自行练习，慢慢就搞清楚了。
- 有的包不是这个网站上面的，安装的方法也不完全一样，如果碰到了，首先还是Bing一下你就知道了，照着前人的和官网的教程就能解决安装及版本转换问题。所以最重要的是搜索和筛选的能力培养。
在这个步骤中我有的包已经安装的。有的还没有安装，而且由于系统更新等问题，重复询问安装等过程也可能会出一些问题，比如这次的结果出现了一句：Registered S3 method overwritten by 'enrichplot'的红字提示，我也没看到waring,也没看到error，先不管，估计他的意思是：我在安装过程中干了这个事情，至于会出现什么问题我也不知道，先提醒你一下。反正没有error,先不管，不要看到红字就血压飙升，头脑发热。 - 对了，看到这里我也需要提示一下，如果是初次使用的话rstudio需要设置镜像，下面有个下载地址就是设置的镜像，这样下载包的速度才能有保证。如果对这些基础背景都不熟悉，我建议跟我一样报个班系统的学习一下比较好。特别是医学生转学生信的，强烈建议！
> rm(list = ls())
> options(stringsAsFactors = F)
> # 数据包： ALL, CLL, pasilla, airway
> # 软件包：limma，DESeq2，clusterProfiler
> # 工具包：reshape2
> # 绘图包：ggplot2
> if (!requireNamespace("BiocManager", quietly = TRUE))
+ install.packages("BiocManager")
> if (!requireNamespace("ALL", quietly = TRUE))
+ BiocManager::install("ALL")
> if (!requireNamespace("CLL", quietly = TRUE))
+ BiocManager::install("CLL")
> if (!requireNamespace("pasilla", quietly = TRUE))
+ BiocManager::install("pasilla")
Bioconductor version 3.9 (BiocManager 1.30.4), R 3.6.1 (2019-07-05)
Installing package(s) 'pasilla'
installing the source package ‘pasilla’
trying URL 'https://mirrors.ustc.edu.cn/bioc/packages/3.9/data/experiment/src/contrib/pasilla_1.12.0.tar.gz'
Content type 'application/gzip' length 3869162 bytes (3.7 MB)
==================================================
downloaded 3.7 MB
* installing *source* package ‘pasilla’ ...
** using staged installation
** data
** inst
** help
*** installing help indices
** building package indices
** installing vignettes
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (pasilla)
The downloaded source packages are in
‘/private/var/folders/tm/q03dw2_n18v2v6nlsw81yhc80000gn/T/RtmpRtpr89/downloaded_packages’
Update old packages: 'AnnotationDbi', 'backports', 'bayestestR', 'BiocManager',
'BiocParallel', 'biomaRt', 'callr', 'circlize', 'classInt', 'curl', 'data.table',
'DelayedMatrixStats', 'DescTools', 'devtools', 'digest', 'DT', 'ellipsis', 'emmeans',
'ensembldb', 'GenomicRanges', 'ggstatsplot', 'GlobalOptions', 'HDF5Array', 'htmlTable',
'htmltools', 'htmlwidgets', 'httpuv', 'insight', 'IRanges', 'KernSmooth', 'knitr',
'lambda.r', 'later', 'logspline', 'maptools', 'matrixStats', 'metaBMA', 'mgcv',
'org.Mm.eg.db', 'org.Rn.eg.db', 'pairwiseComparisons', 'pkgbuild', 'pkgconfig', 'polspline',
'promises', 'purrr', 'qvcalc', 'rcompanion', 'RcppAnnoy', 'RcppArmadillo', 'RcppHNSW',
'recipes', 'rhdf5', 'Rhdf5lib', 'Rhtslib', 'RJSONIO', 'rmarkdown', 'Rsamtools',
'rstantools', 'rvcheck', 'S4Vectors', 'shiny', 'sjlabelled', 'sjmisc', 'sjstats',
'StanHeaders', 'SummarizedExperiment', 'survminer', 'tidyr', 'tinytex', 'xfun'
Update all/some/none? [a/s/n]: if (!requireNamespace("airway", quietly = TRUE))
Update all/some/none? [a/s/n]:
n
> if (!requireNamespace("airway", quietly = TRUE))
+ BiocManager::install("airway")
> if (!requireNamespace("limma", quietly = TRUE))
+ BiocManager::install("limma")
> if (!requireNamespace("DESeq2", quietly = TRUE))
+ BiocManager::install("DESeq2")
> if (!requireNamespace("clusterProfiler", quietly = TRUE))
+ BiocManager::install("clusterProfiler")
Registered S3 method overwritten by 'enrichplot':
method from
fortify.enrichResult DOSE
> if (!requireNamespace("reshape2", quietly = TRUE))
+ BiocManager::install("reshape2")
> if (!requireNamespace("ggplot2", quietly = TRUE))
+ BiocManager::install("ggplot2")
查看是否安装成功，可以用如下代码，如果不报错，一般没有什么问题
library(ALL)
library(CLL)
library(pasilla)
library(airway)
library(DESeq2)
library(clusterProfiler )
library(reshape2)
library(ggplot2)
2.了解ExpressionSet对象，比如CLL包里面就有data(sCLLex) ，找到它包含的元素，提取其表达矩阵(使用exprs函数)，查看其大小
- 参考：http://www.bio-info-trainee.com/bioconductor_China/software/limma.html
- 参考：https://github.com/bioconductor-china/basic/blob/master/ExpressionSet.md
- 解析：这里已经有提示信息了，先看看曾老师发的2个教程了解到这个表达矩阵的信息是存放在这个
CLL数据包里面，然后我要做的是加载这个数据然后提取有用的信息，提取的函数是exprs(),查看大小可以有很多不同的方式，因为我们要挖掘的信息最终都是要进行过滤和分析的，所以这一步特别的关键！ - 我们找到的数据格式各种各样，但是分析软件要求的格式又有自己的特点，我们想要用差异分析软件来分析数据，就必须先了解他处理数据的格式，而且基本上GEO数据库里面的很多数据都下载方式不同，格式也不一样，后续分析是需要保存为数据包这种对象形式，下游分析需要从这些数据包里面进行提取然后再走流程得到分析结果。因此第一步先认识
ExpressionSet这个对象。 - 之前一直没有理解这个东西，因此跑流程的过程一直出各种问题，所以在此提醒初学者，如果想要系统的入门R，一定要把曾老师的这20个题目吃透，我之前就是没有认真做这些题才发现，踩了好多坑。就光这个
ExpressionSet这个对象就坑了我好多次，而我竟然才意识到！！！ - 完成背景知识的了解，我们就可以照着老师的代码去操作，但是如果不理解即使运行出来了，后面也不会用，这就是我的教训。
- 上面那两个提示的链接最后还是看一下，学习一下，有时间呢最好配合看看B站的视频，真的是良心制作。
做第二题之前先加载一下包，然后发现我上一步遗漏了一个问题这些包的正确运行，需要很多依赖的包，加载的时候就会有提示，如果缺什么就补什么，也不要惊慌。
例如：
> library(CLL)
Loading required package: affy
> BiocManager::install("affy")
Bioconductor version 3.9 (BiocManager 1.30.4), R 3.6.1 (2019-07-05)
Installing package(s) 'affy'
trying URL 'https://mirrors.ustc.edu.cn/bioc/packages/3.9/bioc/bin/macosx/el-capitan/contrib/3.6/affy_1.62.0.tgz'
Content type 'application/octet-stream' length 1648118 bytes (1.6 MB)
==================================================
downloaded 1.6 MB
The downloaded binary packages are in
/var/folders/tm/q03dw2_n18v2v6nlsw81yhc80000gn/T//RtmpRtpr89/downloaded_packages
Update old packages: 'AnnotationDbi', 'backports', 'bayestestR', 'BiocManager',
'BiocParallel', 'biomaRt', 'callr', 'circlize', 'classInt', 'curl', 'data.table',
'DelayedMatrixStats', 'DescTools', 'devtools', 'digest', 'DT', 'ellipsis', 'emmeans',
'ensembldb', 'GenomicRanges', 'ggstatsplot', 'GlobalOptions', 'HDF5Array', 'htmlTable',
'htmltools', 'htmlwidgets', 'httpuv', 'insight', 'IRanges', 'KernSmooth', 'knitr',
'lambda.r', 'later', 'logspline', 'maptools', 'matrixStats', 'metaBMA', 'mgcv',
'org.Mm.eg.db', 'org.Rn.eg.db', 'pairwiseComparisons', 'pkgbuild', 'pkgconfig', 'polspline',
'promises', 'purrr', 'qvcalc', 'rcompanion', 'RcppAnnoy', 'RcppArmadillo', 'RcppHNSW',
'recipes', 'rhdf5', 'Rhdf5lib', 'Rhtslib', 'RJSONIO', 'rmarkdown', 'Rsamtools',
'rstantools', 'rvcheck', 'S4Vectors', 'shiny', 'sjlabelled', 'sjmisc', 'sjstats',
'StanHeaders', 'SummarizedExperiment', 'survminer', 'tidyr', 'tinytex', 'xfun'
Update all/some/none? [a/s/n]:
n
>
- 加载完成这个包就可以开始解题了：
data("sCLLex")
class(sCLLex)
typeof(sCLLex)
dim(sCLLex)
str(sCLLex)
?sCLLex
通过上面的函数就能慢慢去认识这个数据了，这个对象还是蛮复杂的，就当他是个图书馆好了，问号以后就发现有帮助文档协助我们来按图索骥了,帮助文档里面提供了很多函数来协助我们拿到想要的数据，像图书馆的索引一样。
- 查看内容：
> sCLLex
ExpressionSet (storageMode: lockedEnvironment)
assayData: 12625 features, 22 samples
element names: exprs
protocolData: none
phenoData
sampleNames: CLL11.CEL CLL12.CEL ... CLL9.CEL (22 total)
varLabels: SampleID Disease
varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu95av2
这个数据集中有22个样本，12625个基因；使用exprs可以查看这个数据集的表达水平；使用phenoData可以查看这个数据集的分组信息，如下所示：
#查看对象信息
> sCLLex
ExpressionSet (storageMode: lockedEnvironment)
assayData: 12625 features, 22 samples
element names: exprs
protocolData: none
phenoData
sampleNames: CLL11.CEL CLL12.CEL ... CLL9.CEL (22 total)
varLabels: SampleID Disease
varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu95av2
#提取表达矩阵展示前4行前4列
> eset <- exprs(sCLLex)
> eset[1:4,1:4]
CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL
1000_at 5.743132 6.219412 5.523328 5.340477
1001_at 2.285143 2.291229 2.287986 2.295313
1002_f_at 3.309294 3.318466 3.354423 3.327130
1003_s_at 1.085264 1.117288 1.084010 1.103217
#查看表型信息
> grouplist2 <- phenoData(sCLLex)
> head(grouplist2)
An object of class 'AnnotatedDataFrame'
sampleNames: CLL11.CEL CLL12.CEL ... CLL16.CEL (6 total)
varLabels: SampleID Disease
varMetadata: labelDescription
#提取表型信息
> grouplist1 <- pData(sCLLex)
> head(grouplist1)
SampleID Disease
CLL11.CEL CLL11 progres.
CLL12.CEL CLL12 stable
CLL13.CEL CLL13 progres.
CLL14.CEL CLL14 progres.
CLL15.CEL CLL15 progres.
CLL16.CEL CLL16 progres.
#查看样本名称
> sampleNames(sCLLex)
[1] "CLL11.CEL" "CLL12.CEL" "CLL13.CEL" "CLL14.CEL" "CLL15.CEL" "CLL16.CEL" "CLL17.CEL"
[8] "CLL18.CEL" "CLL19.CEL" "CLL20.CEL" "CLL21.CEL" "CLL22.CEL" "CLL23.CEL" "CLL24.CEL"
[15] "CLL2.CEL" "CLL3.CEL" "CLL4.CEL" "CLL5.CEL" "CLL6.CEL" "CLL7.CEL" "CLL8.CEL"
[22] "CLL9.CEL"
#查看所有表型变量
> varMetadata(sCLLex)
labelDescription
SampleID Sample ID
Disease Stable/Progressive
#查看前5行探针信息
> featureNames(sCLLex)[1:5]
[1] "1000_at" "1001_at" "1002_f_at" "1003_s_at" "1004_at"
- 在上面的信息中我们看到
Annotation: hgu95av2指的hgu95av2是一个注释包，后面会用到。 - 然后是查看是否有重复的编码，这里可以用到一个神奇的符号
%>%,加载任意一个tidyverse包就可以调用(在Rstudio，如果是Mac系统请使用快捷键cmd+shit+M,如果是win系统，就改成ctr+shit+M)。
> library(dplyr)#目的是使用 %>%
> featureNames(sCLLex) %>% unique() %>% length()
[1] 12625
关于提取表型信息，参考前人的作业
pdata <- pData(sCLLex)
head(pdata)
group_list=as.character(pdata$Disease)
#从数据框中只要表型信息
table(group_list)
#统计表型信息
- 这个操作对于作图比较重要，我之前就是没有理解这个问题，2，3天都没想明白，我明明做了一个跟人家看起来一模一样的数据就是跑不出类似的图出来，其实看起来一样的东西往往在组成的时候是不一样的，结果当然就不同了。这里的操作把分组信息factor化了，我也不知道怎么解释这个东西，反正就是变成一个变量了。
- 得到表达矩阵和分组信息了就可以绘制芯片数据的质量值，如下所示：
> y <- melt(as.data.frame(eset))
No id variables; using all as measure variables
> p <- ggplot(data=y,aes(x=variable,y=value))
> p <- p + geom_boxplot(aes(fill=variable))
> p <- p + theme(axis.text.x = element_text(angle=45,hjust=1,vjust=1))
> p <- p + xlab("group") + ylab("express") + guides(fill=FALSE)
> p
这里有个提示，但是没有waring,也没有error，就不用管了，我还是遇到的不够多，开始还以为有什么问题，反复测试了，没啥问题，运行结束得到美图如下
2019-9-19
3.了解str,head,help函数，作用于第二步提取到的表达矩阵
这几个函数就是帮助我们学习的，多使用，曾老师说要用到超过100次以上就入门了，我还差得远。上面已经使用过了，这里就不重复。
4.安装并了解hgu95av2.db包,看看ls("package:hgu95av2.db")后显示的那些变量
学习一下db数据格式的相关知识
hgu95av2.db是一个注释包，它为hgu95av2平台的芯片提供注释，这个包中有很多注释文件，如下所示：
先安装包再加载数据，查看帮助文档，找到调用函数：
BiocManager::install("hgu95av2.db")
library(hgu95av2.db)
?hgu95av2.db
ls("package:hgu95av2.db")
> ls("package:hgu95av2.db")
[1] "hgu95av2" "hgu95av2_dbconn" "hgu95av2_dbfile"
[4] "hgu95av2_dbInfo" "hgu95av2_dbschema" "hgu95av2.db"
[7] "hgu95av2ACCNUM" "hgu95av2ALIAS2PROBE" "hgu95av2CHR"
[10] "hgu95av2CHRLENGTHS" "hgu95av2CHRLOC" "hgu95av2CHRLOCEND"
[13] "hgu95av2ENSEMBL" "hgu95av2ENSEMBL2PROBE" "hgu95av2ENTREZID"
[16] "hgu95av2ENZYME" "hgu95av2ENZYME2PROBE" "hgu95av2GENENAME"
[19] "hgu95av2GO" "hgu95av2GO2ALLPROBES" "hgu95av2GO2PROBE"
[22] "hgu95av2MAP" "hgu95av2MAPCOUNTS" "hgu95av2OMIM"
[25] "hgu95av2ORGANISM" "hgu95av2ORGPKG" "hgu95av2PATH"
[28] "hgu95av2PATH2PROBE" "hgu95av2PFAM" "hgu95av2PMID"
[31] "hgu95av2PMID2PROBE" "hgu95av2PROSITE" "hgu95av2REFSEQ"
[34] "hgu95av2SYMBOL" "hgu95av2UNIGENE" "hgu95av2UNIPROT"
5.理解 head(toTable(hgu95av2SYMBOL))的用法，找到 TP53 基因对应的探针ID
> head(toTable(hgu95av2SYMBOL))
probe_id symbol
1 1000_at MAPK3
2 1001_at TIE1
3 1002_f_at CYP2C19
4 1003_s_at CXCR5
5 1004_at CXCR5
6 1005_at DUSP1
> hgu95av2SYMBOL
SYMBOL map for chip hgu95av2 (object of class "ProbeAnnDbBimap")
这里提示hgu95av2SYMBOL是个对象，关于一个chip芯片平台得到的一个东西，那么既然是对象，就需要向前面的第二题一样去找相关资料来理解和提取里面的数据了。
?hgu95av2SYMBOL
# hgu95av2SYMBOL is an R object that provides mappings between manufacturer identifiers and gene abbreviations.
帮助文档帮我们描述了一下这个信息，然后应该也能推理出来，调用里面的数据的函数了，这个toTable()估计就跟上面的exprs()作用类似了，通过类比，我们就可以学习更多的对象了。具体的操作就是看说明书操作了，那么就我这种学了一个月的人来讲，做到现在这几道题已经花费了我超过2个小时的时间，但是在最初学习的时候我就花了一天做了前面的10道题，只是跟着流程走，不太明白具体的东西，但是这个简单的操作背后确实需要积累非常多的知识，需要一步步扎实的走，不然就像我一样，走到一半又回头来重新开始。
- 理解了这个背景之后来解题
# 查看hgu95av2SYMBOL大致信息
summary(hgu95av2SYMBOL)
gene_id <- toTable(hgu95av2SYMBOL)
head(gene_id)
str(gene_id)
#查找一个名字叫"TP53"的gene
filter(gene_id,symbol=="TP53")
- 这里就用到了
filter()过滤函数，问号一下这个函数发现它存在于2个R包，以上stats另一个是上面演示管道符时候加载的一个包dplyr。 - 所有很多只要你想到的功能，前人都已经想好了怎么去化繁为简，基本上每个英语单词都是一个函数，基本上对应它本身的功能。
6. 理解探针与基因的对应关系，总共多少个基因，基因最多对应多少个探针，是哪些基因，是不是因为这些基因很长，所以在其上面设计多个探针呢？
这个就需要去查阅相关的资料了，答案都能搜得到比如这里我就引用前人的答案了，当然自己也需要去理解的。
答：不管是Agilent芯片，还是Affymetrix芯片，上面设计的探针都非常短。最长的如Agilent芯片上的探针，往往都是60bp，但是往往一个基因的长度都好几Kb。因此一般多个探针对应一个基因，取最大表达值探针来作为基因的表达量。（也可以取均值。）
7. 第二步提取到的表达矩阵是12625个探针在22个样本的表达量矩阵，找到那些不在hgu95av2.db包收录的对应着SYMBOL的探针。
提示：有1165个探针是没有对应基因名字的。
解析：这里回到第二题，找到相应的信息
> sCLLex
ExpressionSet (storageMode: lockedEnvironment)
assayData: 12625 features, 22 samples
element names: exprs
protocolData: none
phenoData
sampleNames: CLL11.CEL CLL12.CEL ... CLL9.CEL (22 total)
varLabels: SampleID Disease
varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu95av2
> eset <- exprs(sCLLex)
> eset[1:4,1:4]
CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL
1000_at 5.743132 6.219412 5.523328 5.340477
1001_at 2.285143 2.291229 2.287986 2.295313
1002_f_at 3.309294 3.318466 3.354423 3.327130
1003_s_at 1.085264 1.117288 1.084010 1.103217
> head(toTable(hgu95av2SYMBOL))
probe_id symbol
1 1000_at MAPK3
2 1001_at TIE1
3 1002_f_at CYP2C19
4 1003_s_at CXCR5
5 1004_at CXCR5
6 1005_at DUSP1
> probe_map <- hgu95av2SYMBOL
> length(probe_map)
[1] 12625
> probe_info <- mappedkeys(probe_map)
> length(probe_info)
[1] 11460
> 12625-11460
[1] 1165
查阅帮助文档，了解到如下信息：hug95av2SYMBOL是一个R对象，它提供的是芯片生产厂家与基因缩写之间的映射信息。这个映射的信息主要依据Entrez Gene数据库。可以通过mappedkeys()这个函数，得到映射到基因上的探针信息。
- 可以看到这个数据包里有1165个探针是没有Entrez Gene数据库的映射基因的。
- 现在题目要求是提取不在这个hgu95av2.db包里，但是对应了基因SYMBOL的探针
- 第一步 找到这个包里面所有对应了基因SYMBOL的探针，同时过滤第二步表达矩阵中对应了基因SYMBOL的探针，得到的2个矩阵都是11460行
#法1
probe_map <- hgu95av2SYMBOL
length(probe_map)
probe_info <- mappedkeys(probe_map)
length(probe_info)
gene_info <- as.list(probe_map[probe_info])
# 转化为数据表
length(gene_info)
gene_symbol <- toTable(probe_map[probe_info])
dim(eset)
head(gene_symbol)
#对原始表达矩阵过滤
data_filter <- eset[rownames(eset)%in%probe_info,]
table(rownames(eset)%in%probe_info)
table(rownames(eset) %in% rownames(data_filter))
#法2：
#参考.目的R语言20练习题【完整版】对表达矩阵进行过滤，去掉没有map的探针
#table(rownames(eset)%in%gene_symbol$probe_id) #找到sCLLex表达矩阵(eset)在hgu95av2.db包中没有交叉的探针
#e1 = eset[rownames(eset)%in%probe_info,]
#e2 = eset[match(rownames(eset), probe_info, nomatch = 0),] #使用match过滤
注意区分
%>%是管道符号，相当于linux的|起传递作用，前面已经提到过; 而%in%表示两者求交集
- 第二步 将前面得到的2个矩阵进一步过滤
# 对过滤后的表达值的数据根据探针名称进行排序
data_filter_order <- data_filter[order(rownames(data_filter)),]
# 对gene_symbol按照探针名称进行排序
gene_symbol_order <- gene_symbol[order(gene_symbol$probe_id),]
# 查看地表达值的矩阵的名称是否与gene_symbol中的探针名称是否一一对应
table(gene_symbol_order$probe_id == rownames(data_filter_order))
#由于有的基因名对应了多个探针，用取最大平均值的方式去掉重复的探针
row_mean <- apply(data_filter_order,1,mad)
head(row_mean)
#构建 index_probe 函数，这个函数的功能在于：用gene_symbol$symol（基因名）对row_mean进行分组，然后挑出相同基因中表达值最大的那个值所对应的探针编号
index_probe <- function(x){
names(x)[which.max(x)]
}
#tapply() 函数对矩阵进行计算，lapply()对列表进行计算
max_index <- tapply(row_mean,gene_symbol$symbol,index_probe)
# gene_symbol_order与data_filter_order都按照探针名称进行排序，再挑出max_index中的探针哪些在gene_symbol_order$probe_id中
data_filter_order_unique <- data_filter_order[(gene_symbol_order$probe_id %in% max_index),]
# 将data_filter_order_unique那些与gene_symbol_order$probe_id对应起来的探针与基因名提取出来
unique_map_gene_probe <- gene_symbol_order[(gene_symbol_order$probe_id %in% rownames(data_filter_order_unique)),]
#找到在包里对应了基因symbol的探针有8585个
table(unique_map_gene_probe$probe_id == rownames(data_filter_order_unique))
# 可以发现，unique_map_gene_probe$probe_id中的探针顺序与data_filter_order_unique中的探针顺序是致的
rownames(data_filter_order_unique) <- unique_map_gene_probe$symbol
# 将最终过滤后的矩阵的名称换为基因名
- 按照前人的解题思路最后只筛选出来在这个包里并且有基因symbol的探针，并提取出这些基因symbol对应探针数最多的样本的表达矩阵
- 还是没有能解决不在这个包里面但是又对应了基因symbol的探针的注释，可能是我的理解出现了偏差，既然是不在这个包里面的应该是找不到注释的，如果想找到不在这个包里的探针对应的基因symbol，应该提取那1165个探针的坐标重新注释，然后才能得到不在这个包里但又对应着基因symbol的探针。
- 这时候想到之前有一篇推文（重磅！价值一千元的R代码送给你）芯片探针序列的基因组注释
- 看了一下，都是R语言的操作可以尝试一下。
到这里就是前面10题的总结了，7-10题的解题还需要再练习，特别是第7题还没解完，待我先练习一下再来完成。
