《Modern Statistics for Modern Bi

《Modern Statistics for Modern Biology》Chapter 一： 离散数据模型的预测（1.1 - 1.3）

《Modern Statistics for Modern Biology》Chapter 一： 离散数据模型的预测（1.4 - 1.5）

《Modern Statistics for Modern Biology》Chapter 二: 统计建模（2.1-2.3）

《Modern Statistics for Modern Biology》Chapter 二: 统计建模（2.4-2.5）

《Modern Statistics for Modern Biology》Chapter 二 统计建模（2.6 - 2.7）

《Modern Statistics for Modern Biology》Chapter 二 统计建模（2.8 - 2.9）

《Modern Statistics for Modern Biology》Chapter 二 统计建模（2.10 完结）

从这章开始最开始记录一些markdown等小知识
$\hat{p}=\frac{1}{12}$ ：
掌握R语言中的apply函数族
卡方检验
Hardy-Weinberg equilibrium( 哈迪-温伯格平衡 )
带你理解beta分布
简单介绍一下R中的几种统计分布及常用模型

• 统计分布每一种分布有四个函数：d――density（密度函数），p――分布函数，q――分位数函数，r――随机数函数。比如，正态分布的这四个函数为dnorm，pnorm，qnorm，rnorm。下面我们列出各分布后缀，前面加前缀d、p、q或r就构成函数名：norm：正态，t：t分布，f：F分布，chisq：卡方（包括非中心） unif：均匀，exp：指数，weibull：威布尔，gamma：伽玛，beta：贝塔 lnorm：对数正态，logis：逻辑分布，cauchy：柯西， binom：二项分布，geom：几何分布，hyper：超几何，nbinom：负二项，pois：泊松 signrank：符号秩，wilcox：秩和，tukey：学生化极差
  如何预测一条序列是否含有CpG岛
  图片输出尽量保存为矢量图

3、R语言中的高质量图形

• 哈哈，终于来到最喜欢的环节了。

• 数据可视化分为两类（至少）。第一类能让科学家探索数据并发现工作中的复杂过程；另一类的可视化类型提供了关于他的结果的信息丰富、清晰和具有视觉吸引力的插图，他可以将其展示给他人，并且最终包含在出版物中。
  .....此处省略

3.1 这一章的目标

• 学习怎么通过可视化快速而又灵活的解释数据集
• 画出可发表的高质量图片
• 综述R语言的基础绘图
• ** 理解图形概念语法背后的逻辑**
• 介绍ggplot2绘图功能
• 了解如何在一维，两维甚至三维维度中绘制数据，并探索分面facet
• 为分子生物学数据（或沿着其他序列，例如肽）创建“along-genome”图。
• 讨论我们的一些交互式图形选项

3.2 基础R绘图

• 最基础的绘图是plot函数。在下面的代码中，绘图结果见图 3.2。它用于绘制来自酶联免疫吸附测定（ELISA）测定的数据。该方法用于定量测定脱氧核糖核酸酶(DNase)降解DNA的活性。这些数据被组装在R对象DNase中，该对象DNase是一个数据文件，检测运行；conc: 蛋白质浓度；density: 被测量的光密度。

> head(DNase)
  Run       conc density
1   1 0.04882812   0.017
2   1 0.04882812   0.018
3   1 0.19531250   0.121
4   1 0.19531250   0.124
5   1 0.39062500   0.206
6   1 0.39062500   0.215
> plot(DNase$conc, DNase$density)

图 3.2

• 进一步美化，通过xlab、 ylab来改变xy轴的坐标名，pch改变绘制图形的形状（比如：正方形、圆形）。有关变量的信息存储在DNase中，我们可以使用attr函数访问它。图 3.3
  pch参数对应的图形形状

> plot( DNase$conc, DNase$density, 
+ ylab = attr(DNase, "labels")$y,
+ xlab = paste(attr(DNase, "labels")$x, attr(DNase, "units")$x),
+ pch = 3,
+ col = "blue")
> attr(DNase, "labels")
$`x`
[1] "DNase concentration"

$y
[1] "Optical density"
> attr(DNase, "units")
$`x`
[1] "(ng/ml)"

图 3.3 这里严重没搞懂，明明横坐标12.5的位置就对应两个值，这里怎么会有这么多值
横坐标12.5

► 问题

• 用 “metadata”(如较长的描述、物理单元、来源信息等)对数据文件列进行注释似乎是一个有用的特性。这种存储此类信息的方式(如DNase对象)在R生态系统中是标准化的还是通用的？有没有其他标准化或通用的方法来实现这一点？
  ► 解
• 在常规的R data.frame中没有好的或广泛使用的基础结构，在tidyverse(data_frame，tibble)中也没有。但是看看Bioconductor包S4Vectors中的DataFrame类。除其他外，它用于注释SummarizedExperiment的行和列。
• 除了散点图，我们也可以使用柱形图和箱式图来表示 （图 3.4）

> hist(DNase$density, breaks=25, main = "")
> boxplot(density~Run, data = DNase)

图 3.4

• Boxplots便于在紧致空间中显示彼此相邻的多个分布。我们将在第3.6节中看到更多关于绘制多个单变量分布的内容。
• 基本的R绘图功能对于快速交互式地探索数据非常有用；但是，如果我们想要创建更复杂的显示，我们很快就会遇到它们的局限性。我们将使用一个称为图形语法的可视化框架，该框架在ggplot2包中实现，能够以逻辑和优雅的方式逐步构建高质量的图形。首先，让我们介绍并加载一个示例数据集。

3.3 示例数据集

• 为了正确地测试ggplot2的功能，我们需要一个足够大且具有一定复杂性的数据集，以便可以从许多不同的角度对其进行切片和查看。我们将使用一个基因表达微阵列数据集，报告大约100个小鼠胚胎细胞在早期发育的不同时间点的转录情况。哺乳动物的胚胎最初是单个细胞，即受精卵。通过同步的细胞分裂波，卵细胞增殖成一团细胞，一开始它们之间没有明显的区别。不过，在某些时候，细胞会选择不同的血统。通过进一步的规范，不同的细胞类型和组织出现，这是一个完整的有机体所需要的。实验的目的，由Ohnishi等人解释(2014)，旨在研究胚胎中与第一次symmetry breaking event相关的基因表达变化。我们将进一步解释这些数据。更多细节可以在论文和生物导体数据包 Hiiragi2013的文档中找到。我们首先加载数据：

> BiocManager::install("Hiiragi2013", version = "3.8")
> library("Hiiragi2013")  ## 数据比较大，如果安装失败，重复运行上面代码几次就好
> data("x")
> dim(Biobase::exprs(x))
[1] 45101   101

• 只需在R提示符下键入x，就可以打印出ExpressionSet对象x的更详细摘要。数据矩阵的101列(上面通过Biobase包中的exprs函数访问)对应于样本(每一列都是一个单细胞)，45101行对应于阵列(Affymetrix mouse4302阵列)探测到的基因。使用RMA方法(Irizarry等人)对数据进行归一化处理。(2003年)。原始数据也可在软件包(数据对象a)和EMBL-EBI的ArrayExpress数据库中查阅，登录号为E-MTAB-1681。
• 我们可以看下这个样品中包含什么信息。

> head(pData(x), n = 2)
        File.name Embryonic.day Total.number.of.cells lineage genotype   ScanDate sampleGroup sampleColour
1 E3.25  1_C32_IN         E3.25                    32               WT 2011-03-16       E3.25      #CAB2D6
2 E3.25  2_C32_IN         E3.25                    32               WT 2011-03-16       E3.25      #CAB2D6

• 所提供的信息包括有关细胞的信息(即获得这些细胞的胚胎的年龄、大小和基因型)和技术信息(扫描日期、原始数据文件名)。按照惯例，小鼠胚胎发育的时间以天为单位计算，并报告为，例如，E3.5。此外，在这篇论文中，作者根据年龄、基因型和血统将细胞划分为8个生物组(Samplegroup)，并定义了一个颜色方案来表示这些组(SampleColour)。使用下面的代码(请参阅下面的解释)，我们定义了一个小的数据文件组，其中包含每个组的摘要信息：单元格的数量和首选颜色。

> library("dplyr")
> groups = group_by(pData(x), sampleGroup) %>%
+   summarise(n = n(), color = unique(sampleColour))
> groups
# A tibble: 8 x 3
  sampleGroup         n color  
  <chr>           <int> <chr>  
1 E3.25              36 #CAB2D6
2 E3.25 (FGF4-KO)    17 #FDBF6F
3 E3.5 (EPI)         11 #A6CEE3
4 E3.5 (FGF4-KO)      8 #FF7F00
5 E3.5 (PE)          11 #B2DF8A
6 E4.5 (EPI)          4 #1F78B4
7 E4.5 (FGF4-KO)     10 #E31A1C
8 E4.5 (PE)           4 #33A02C

• 在FGF4-KO基因缺失的胚胎中，FGF4基因是一种重要的细胞分化调节基因。从E3.5开始，野生型细胞(没有FGF4基因敲除)经历第一次对称破缺事件，并分化为不同的细胞系，称为多能上胚层(EPI)和原始内胚层(PE)。

• 由于上面的代码块是我们遇到管道操作符%>%和函数group_by以及包dplyr中的summarise函数的第一个实例，所以让我们分布解析代码。首先，管道%>%。通常，管道对于使嵌套函数调用更易于读取非常有用。在R中下面的两行代码行驶的功能是一样的.

f(x) %>% g(y) %>% h
h(g(f(x), y))

• 也就是说：“求值f(x)，然后将结果作为第一个参数传递给函数g，而将y作为第二个参数传递给g。然后将g的输出传递给函数h。你可以无限地重复这个。特别是如果参数x和y本身是复杂的表达式，或者有相当多的函数链，第一个版本往往更容易阅读。

• group_by函数简单地“标记”数据文件，并指出所有后续操作不应同时应用于整个数据文件，而应用于由sampleGroup因子定义的块。最后，summarise计算汇总统计信息；这可以是求平均值、和；在本例中，我们只计算每个块中的行数n()和流行的颜色。

3.4 ggplot2

• ggplot2是** Hadley Wickham (Wickham 2016)** 创作的。balabala......
• 现在我们使用ggplot2重新绘制图 3.2。

> library(ggplot2)
> ggplot(DNase, aes(x = conc, y = density)) + geom_point()

图 3.6

• balabal..... 现在，让我们转到小鼠单细胞数据，并使用ggplot函数为8个组中的每个组绘制样本数。结果如图3.7所示。

> ggplot(groups, aes(sampleGroup, n)) + geom_bar(stat = "identity")

图 3.7

接下来省略中间话语，主要注重于代码的学习

> groupColor = setNames(groups$color, groups$sampleGroup)
>  groupColor 
          E3.25 E3.25 (FGF4-KO)      E3.5 (EPI)  E3.5 (FGF4-KO)       E3.5 (PE)      E4.5 (EPI)  E4.5 (FGF4-KO) 
      "#CAB2D6"       "#FDBF6F"       "#A6CEE3"       "#FF7F00"       "#B2DF8A"       "#1F78B4"       "#E31A1C" 
      E4.5 (PE) 
      "#33A02C" 
> ggplot(groups, aes(x = sampleGroup, y = n, fill = sampleGroup)) +
+   geom_bar(stat = "identity") +
+   scale_fill_manual(values = groupColor, name = "Groups") +
+   theme(axis.text.x = element_text(angle = 90, hjust = 1))

图 3.8

3.4.1 Data flow

> gg = ggplot(DNase, aes(x = conc, y = density)) + geom_point()
> gg
> print(gg)

3.4.2 保存图片

> ggsave("DNAse-histogram-demo.pdf", plot = gg)

• 存储图形的主要方法有两种：矢量图形和光栅(像素)图形。在矢量图形中，打印存储为一系列几何图元，如点、直线、曲线、形状和印刷字符。R中用于将打印保存为矢量图形格式的首选格式是PDF。在光栅图形中，打印存储在点阵数据结构中。光栅格式的主要限制是分辨率有限，这取决于可用的像素数。在R中，最常用的光栅图形输出设备是PNG。**通常，最好将打印保存为矢量图形格式**，因为以后总是可以将矢量图形文件转换为任何所需分辨率的光栅格式，而反向转换则相当困难。你不希望你的演讲或论文中的数字看起来很糟糕，因为像素化的人工制品！