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%以上的观测

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

最后才开始画火山图

  • 发表于:
  • 原文链接https://kuaibao.qq.com/s/20180903G1YFX000?refer=cp_1026
  • 腾讯「云+社区」是腾讯内容开放平台帐号(企鹅号)传播渠道之一,根据《腾讯内容开放平台服务协议》转载发布内容。

扫码关注云+社区

领取腾讯云代金券

玩转腾讯云 有奖征文活动