R语言20练习题完整版

今天是生信星球陪你的第116天

你想找辆共享单车，发现满街都是别家车，没有一辆你能骑。

你想学点生信，搜了“初学者教程”，满眼尽是高大上，没有一句能看懂。

终于你跨越茫茫宇宙，来到生信星球，发现了初学者的新大陆

豆豆重新写于2018.9.3

感谢生信技能树jimmy出的这么好一份R练习题，http://www.bio-info-trainee.com/3409.html

在此基础上进行加工、重构，加入了自己的理解

.1 安装R包

.2 了解ExpressionSet对象

CLL包中有data(sCLLex)，是一个表达芯片数据对象，其中包含许多信息

第一行的ExpressionSet就是表达矩阵，查看它使用 ，用 查看矩阵大小

使用str查看对象的结构，使用head查看对象的前6行（默认）

.3 安装并了解hgu95av2.db包

官网：

http://www.bioconductor.org/packages/release/data/annotation/html/hgu95av2.db.html

安装

这个数据库中共有36个包，每个包都可以当成一个列表操作，可以用 函数展示数据

3.1 了解hgu95av2.db 【这是一个关于hgu95av2芯片的注释包】

3.2 探索hgu95av2CHR：探针id与染色体编号对应的关系

3.3 探索hgu95av2SYMBOL：探针id与基因名（缩写）的关系3.3.1 先对探针id进行简单的探索：【输出探针】

3.3.2 再对基因名进行探索【匹配、列表、转数据框、查特定基因、基因数】

为什么有5个基因会分别有8个探针，而大部分6555个基因只对应一个探针？

A：不管是Agilent芯片，还是Affymetrix芯片，上面设计的探针都非常短。最长的如Agilent芯片上的探针，往往都是60bp，但是往往一个基因的长度都好几Kb。因此一般多个探针对应一个基因，取最大表达值探针来作为基因的表达量

.4 过滤、整合表达矩阵

4.1 过滤

4.2 整合

4.3 重塑【reshape2:矩阵=》数据框】

.5 画图探索

作出样本和基因表达量之间的关系图，主要基于ggplot

每种图都做了两种版本，一个初始版，一个调整版

5.1 boxplot

5.2 violin plot

5.3 Histogram

5.4 Density

.6 做一些统计

mean,median,max,min,sd,var,mad，t检验

6.1 利用apply函数进行统计

他需要矩阵，也就是之前得到的efilt，按行进行统计即可

6.2 看表达量top30基因之间的重合部分

用UpSetR包结合之前做的top30基因各种统计，适用于样本数量大于5的情况

其实我们平常做的韦恩图也是这个意思，找交集，但是韦恩图样本数量一般都会控制在

6.3 批量T检验——为了得到pvalue【后续分析重点】

有了pvale就能有padj值；另外还需要对照、处理两组均值，这样就能有log2FC

一般来讲，下游分析使用的差异表达矩阵应该是log2后的结果，它的计算公式是log2FC = log2 (mean(处理组/对照组))

这里为什么可以之间相减？

芯片数据的特点：小样本和大变量，因此数据分布呈偏态、标准差大。而对数转换能使上调、下调的基因连续分布在0的周围，更加符合正态分布，同时对数转换可以使荧光信号强度的标准差减少，方便下游分析

因此我们一直用的e也就是exprs(sCLLex)得到的表达矩阵是log2之后的

我们得到的eavg_2 = log2（mean处理组），eavg_1 = log2（mean处理组），

根据公式就可以算出log2(a/b) = log2(a) - log2(b)

.7 做一些分析 表达量、聚类、PCA、火山图

7.1 按mad指标选表达量前30（top30）的基因,做热图可视化

做热图前需要将矩阵数据中心化+标准化【目的为了向数据中心靠拢，减小数据之间的差别】

中心化：数据减去均值后得到的; 标准化则是在中心化后的数据基础上再除以数据的标准差

7.2 聚类分析图

过滤后的样本聚类

过滤后的基因聚类

使用factoextra包

如果看到聚类分析的结果枝干太长，那么就要换种聚类方法

7.3 PCA分析

【目的：简化变量的个数】本质是降维，本来应该有22个变量，现在我们变成了22个主成分,一般前面的几个主成分就能解释所有的数据。解释：https://wenku.baidu.com/view/c22d1539a31614791711cc7931b765ce05087a6f.html

http://www.bio-info-trainee.com/1232.html

7.3.1 使用一键式ggfortify + prcomp

关于ggfortify的使用：https://wenku.baidu.com/view/e5dc63fb763231126fdb1100.html

最大的优点：一行代码，出ggplot风格的图，不用费时调整，提高效率

结果中：不同颜色代笔不同分组;

坐标轴：能最大反映样本差异性的两个成分（PC1、PC2）

百分数：成分贡献率;

坐标轴刻度：没实际意义（代表相对距离）

7.3.2 使用fviz_pca_ind包

7.4 差异分析火山图【也就是logFC加上-log10(Pvals)的散点图】

首先构建差异分析矩阵

关于limma：

limma是基于R和Bioconductor平台的分析芯片数据的综合软件包，该包功能齐全、教程完善、使用率极高，几乎成为了芯片数据处理流程的代名词；

本质就是对表达量矩阵做一个归一化，然后利用理想的统计分布检验函数来计算差异的显著性；

imma的核心函数是lmFit和eBayes， 前者是用于线性拟合，后者根据前者的拟合结果进行统计推断；

lmFit至少需要两个输入，一个是表达矩阵，一个是分组对象；

表达矩阵必须是matrix类数据结构，每一列都是存放一个样本，每一行是一个探针信息或者是注释后的基因名

关于比较矩阵: 参考https://github.com/bioconductor-china/basic/blob/master/makeContrasts.md

7.4.1 构建一个非差异比较矩阵【只有一组处理和对照，所以可以不用比较矩阵】

7.4.2 构建差异比较矩阵

7.4.3 准备火山图

画火山图第一步，设定阈值，选出UP、DOWN、NOT表达基因

关于阈值的设置：一般情况都是设为2，但具体背景不同还是应该设置不同的阈值。根据网站https://www.bmj.com/about-bmj/resources-readers/publications/statistics-square-one/2-mean-and-standard-deviation 的介绍。mean+2SD可以反映95%以上的观测值，这是比较可信的，如果再严格一点，设为mean+3SD，就可以反映97%以上的观测

画火山图第二步，设定图形的标题

最后才开始画火山图