药物的作用机制 （MoA）预测

药物的作用机制 （MoA） 是什么？为什么它很重要？

过去，科学家们从天然产物中提取药物或受到传统疗法的启发。非常常见的药物，如扑热息痛，在美国被称为对乙酰氨基酚，在驱动其药理活性的生物学机制被理解之前几十年就已投入临床使用。今天，随着更强大技术的出现，药物发现已经从过去的偶然方法转变为基于对疾病潜在生物学机制的理解的更具针对性的模型。在这个新框架中，科学家们试图确定与疾病相关的蛋白质靶标，并开发一种可以调节该蛋白质靶标的分子。作为描述给定分子生物活性的简写，科学家们分配了一个称为作用机制或简称 MoA 的标签。

我们如何确定新药的 MoAs？

一种方法是用药物处理人类细胞样本，然后使用算法分析细胞反应，这些算法在大型基因组数据库中搜索与已知模式的相似性，例如具有已知 MoA 的药物的基因表达文库或细胞活力模式。

在本次比赛中，您将可以访问一个结合了基因表达和细胞活力数据的独特数据集。该数据基于一项新技术，该技术可同时测量（在同一样本中）人类细胞对 100 种不同细胞类型中药物的反应（从而解决了事前识别哪些细胞类型更适合给定药物的问题）。此外，您还可以访问此数据集中 5,000 多种药物的 MoA 注释。

按照惯例，数据集已分为测试和训练子集。因此，您的任务是使用训练数据集开发一种算法，该算法可自动将测试集中的每个事例标记为一个或多个 MoA 类。请注意，由于药物可以有多个 MoA 注释，因此该任务正式是一个多标签分类问题。

如何评估解决方案的准确性？

基于 MoA 注释，将根据应用于每个药物-MoA 注释对的对数损失函数的平均值来评估解决方案的准确性。

如果成功，您将帮助开发一种算法来预测化合物在给定其细胞特征下的 MoA，从而帮助科学家推进药物发现过程。

1介绍

这是一个与tidy R和ggplot2竞争的全面的探索性数据分析。它应该为你提供开始工作所需的所有信息，甚至可能更多。

这项挑战的目的是“根据生物活性对药物进行分类”。药物发现的目的是识别与特定疾病相关的特定蛋白质，然后开发可以靶向这些蛋白质的分子。分子的MoA编码其生物活性。该数据集描述了100种不同类型的人体细胞对各种药物的反应。这些响应模式将用于对MoA响应进行分类。

这是一个多标签分类问题。药物可以有多个MoA注释，这些注释以不同的方式描述来自不同细胞类型的二元反应。评价指标是列式平均日志损失。

这些数据都是我们熟悉的训练和测试文件。与其他比赛不同，这里我们有两个单独的文件用于训练预测器（train_features.csv）和目标（train_targets_scores .csv）。每一行对应一个特定的处理。此外，我们还给定了一组可选的MoA目标（train_targets_nonscored.csv），我们不需要预测，但可以用于上下文分析。

具体目标是预测测试文件test_features.csv的train_targets_score .csv类概率。请注意，公共排行榜是基于大约25%的测试数据，其中75%分配给私有LB。

让我们开始吧！

2 预处理

2.1 加载库

我们加载了一系列用于通用数据处理和通用可视化的库，以及用于特定工作的更专业的工具。

# general visualisation
library('ggplot2') # visualisation
library('scales') # visualisation
library('patchwork') # visualisation
library('RColorBrewer') # visualisation
library('corrplot') # visualisation
library('ggthemes') # visualisation
library('viridis') # visualisation
library('glue') # visualisation

# general data manipulation
library('dplyr') # data manipulation
library('readr') # input/output
library('vroom') # input/output
library('skimr') # overview
library('tibble') # data wrangling
library('tidyr') # data wrangling
library('stringr') # string manipulation
library('forcats') # factor manipulation
library('janitor') # data cleaning
library('tictoc') # timing

# specific visualisation
library('alluvial') # visualisation
library('ggrepel') # visualisation
library('ggforce') # visualisation
library('ggridges') # visualisation
library('gganimate') # animations
library('GGally') # visualisation
library('wesanderson') # visualisation

# specific data manipulation
library('lazyeval') # data wrangling
library('broom') # data wrangling
library('purrr') # data wrangling
library('reshape2') # data wrangling
library('rlang') # encoding

# dimensionality reduction
library('factoextra')

# modelling
library('recipes')
library('rsample')
library('keras')
library('tfdatasets')

2.2 加载函数

我们使用一个简单的辅助函数来计算二项式置信区间。

# function to extract binomial confidence levels
get_binCI <- function(x,n) as.list(setNames(binom.test(x,n)$conf.int, c("lwr", "upr")))

2.3 加载数据

我们使用vroom包来读取数据：

train <- vroom(str_c(path,'train_features.csv'), col_types = cols())
targets <- vroom(str_c(path, "train_targets_scored.csv"), col_types = cols())
targets_non <- vroom(str_c(path, "train_targets_nonscored.csv"), col_types = cols())
test <- vroom(str_c(path,'test_features.csv'), col_types = cols())
sample_submit <- vroom(str_c(path,'sample_submission.csv'), col_types = cols())

3.文件内容和数据

首先，我们将快速了解数据集及其形状。下面的表是交互式的，涵盖了前50行的所有列：

训练数据：

• 这是一个相当广泛的数据集，几乎有900列。从数据描述中，我们知道以“g-”开头的特征编码基因表达数据（有772个），以“c-”开头的特征（共100个）显示细胞活力数据。
• 此外，我们有3个“cp_”特征：cp_type表示样本处理，而cp_time和cp_dose表示处理的持续时间和剂量。
• sig_id是这个示例的唯一主键。

测试数据：

目标：

• 目标反应也相当广泛；只有200多个不同的二进制输出。正如所料，这张表非常稀疏。
• sig_id可用于将预测器连接到目标。

额外未得分目标：

• 额外的未得分目标包含约400个类别（即列）；是我们得分目标的两倍。
• 这是关于药物本身的很多额外信息。让我们看看能否找到一个好的应用程序。

空缺填充

sum(is.na(train))
sum(is.na(test))

## [1] 0
## [1] 0

目标稀疏

由于目标类别如此之多，moa可能很少见；目标数据框将非常稀疏。数字如下：

non_zero <- targets %>% 
  select(-sig_id) %>% 
  na_if(1) %>% 
  is.na()

non_zero_percent <-  sum(non_zero) / (nrow(non_zero) * ncol(non_zero)) * 100

sprintf("Percentage of non-zero target class values: %.3f%%", non_zero_percent)

## [1] "Percentage of non-zero target class values: 0.343%"

不到0.5%的非零。很不平衡。没有得分的目标呢？

non_zero <- targets_non %>% 
  select(-sig_id) %>% 
  na_if(1) %>% 
  is.na()

non_zero_percent <-  sum(non_zero) / (nrow(non_zero) * ncol(non_zero)) * 100

sprintf("Percentage of non-zero non-scored target class values: %.3f%%", non_zero_percent)

## [1] "Percentage of non-zero non-scored target class values: 0.052%"

这是一个很大的区别。几乎比本来就稀疏的得分目标还要稀疏一个数量级。

质量检查：

首先，一个快速的质量检查，以确保我们的sig_id值确实是唯一的。结果应该是0：

nrow(train) - nrow(train %>% distinct(sig_id))

## [1] 0

nrow(targets) - nrow(targets %>% distinct(sig_id))

## [1] 0

然后，再次检查训练集和目标集的sig_id值是否匹配。Result也应该是0：

train %>% 
  select(sig_id) %>% 
  anti_join(targets, by = "sig_id") %>% 
  nrow()

## [1] 0

4.个体特征可视化

我们首先分别绘制各种预测特征和目标特征的分布，然后再绘制多特征的视觉效果和相关性。在这里，我们逐组处理特征。

4.1治疗功能

这些预测特征更一般地描述了样本是如何被处理的，在剂量，持续时间方面；以及它是“真正的”治疗还是对照。

p1 <- train %>% 
  count(cp_type) %>% 
  add_tally(n, name = "total") %>% 
  mutate(perc = n/total) %>% 
  ggplot(aes(cp_type, perc, fill = cp_type)) +
  geom_col() +
  geom_text(aes(label = sprintf("%s", n)), nudge_y = 0.02) +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_manual(values = c("grey70", "violetred")) +
  theme_hc() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", fill = "State", title = "Sample treatment", subtitle = "(Compound vs Control)")

p2 <- train %>% 
  count(cp_dose) %>% 
  add_tally(n, name = "total") %>% 
  mutate(perc = n/total) %>% 
  ggplot(aes(cp_dose, perc, fill = cp_dose)) +
  geom_col() +
  geom_text(aes(label = sprintf("%s", n)), nudge_y = 0.02) +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_manual(values = c("darkblue", "darkred")) +
  theme_hc() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", fill = "State", title = "Treatment Dose", subtitle = "(high vs low)")

p3 <- train %>% 
  count(cp_time) %>% 
  mutate(cp_time = as.factor(cp_time)) %>% 
  add_tally(n, name = "total") %>% 
  mutate(perc = n/total) %>% 
  ggplot(aes(cp_time, perc, fill = cp_time)) +
  geom_col() +
  geom_text(aes(label = sprintf("%s", n)), nudge_y = 0.01) +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_brewer(type = "seq", palette = "Oranges") +
  theme_hc() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", fill = "State", title = "Treatment duration", subtitle = "(Units of hours)")

p1 + p2 + p3

绝大多数治疗是复合治疗（“trt_cp”），相比之下，约8%的对照摄动治疗（“ctl_vehicle”）。控制没有MoAs。

治疗剂量分为两类，D1和D2，分别表示高剂量和低剂量。这些指标大致均衡，24小时、48小时或72小时的3个治疗时间类别也是如此。

4.2 基因表达特性

这些本质上是匿名特征，标记为从“g-0”到“g-771”。它们的值是数值型的，所以让我们看一下前4个基因特征的密度作为示例：

train %>% 
  select(sig_id, starts_with("g-")) %>% 
  select(seq(1,5)) %>% 
  pivot_longer(starts_with("g-"), names_to = "feature", values_to = "value") %>% 
  ggplot(aes(value, fill = feature)) +
  geom_density() +
  facet_wrap(~ feature) +
  theme_minimal() +
  theme(legend.position = "none") +
  labs(x = "", y = "", fill = "State", title = "Distributions for gene expression features")

这些分布看起来很正态，这很好。它们中有一些是倾斜的，但没有什么值得转换。

在dplyr 1.0版本（及之后）中，我们可以通过在summarise中使用新的across方法优雅地推导出这些分布的一些统计度量（即最小值、最大值、平均值和标准差）。（以前，我们会使用summarise_at或summarise_all。）这些是这些统计数据的分布；注意，每个切面都有自己的、经过裁剪的轴范围：

gene_stats <- train %>% 
  select(starts_with("g-")) %>% 
  summarise(across(everything(), list(min = min, max = max, mean = mean, sd = sd))) %>% 
  pivot_longer(starts_with("g-"), names_to = "features", values_to = "values") %>% 
  separate(features, into = c("features", "stat"), sep = "_")

gene_stats %>% 
  ggplot(aes(values, fill = stat)) +
  geom_density() +
  scale_fill_manual(values = wes_palette("GrandBudapest2")) +
  facet_wrap(~ stat, scales = "free") +
  theme_tufte() +
  theme(legend.position = "none") +
  labs(x = "", y = "", fill = "State", title = "Gene distribution meta statistics")

• 均值很好地分布在0附近；标准差主要在0.5和1.5之间。
• 最小值和最大值是彼此的镜像。在正/负9 - 10的范围内有显著的增长。

4.3 细胞存活特性

与基因特征类似，细胞活力特征是匿名的，标记为“c-0”到“c-99”；100的特性。其分布如下所示：

train %>% 
  select(sig_id, starts_with("c-")) %>% 
  select(seq(1,5)) %>% 
  pivot_longer(starts_with("c-"), names_to = "feature", values_to = "value") %>% 
  ggplot(aes(value, fill = feature)) +
  geom_density() +
  scale_fill_brewer(palette = "Set3") +
  facet_wrap(~ feature) +
  theme_minimal() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", fill = "State", title = "Distributions for cell viability features")

我们发现:

这也很正常，但在-10附近有明显的起伏。这些是真的吗？值得在下文中更详细地研究。

与基因数据相反，这些分布都没有接近正10的值。然而，这可能是一个采样效应，因为面板和基因方面在每个四个特征内都被缩放到全局最小和最大。

让我们放大负尾，并添加另外两个特征：

train %>% 
  select(sig_id, starts_with("c-")) %>% 
  select(seq(1,7)) %>% 
  pivot_longer(starts_with("c-"), names_to = "feature", values_to = "value") %>% 
  filter(value < -4) %>% 
  ggplot(aes(value, fill = feature)) +
  geom_density() +
  scale_fill_brewer(palette = "Set3") +
  facet_wrap(~ feature) +
  theme_minimal() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", title = "Cell viability features - zoom in on negative tail")

我们发现:

• 那些尾巴绝对是重要的。即使在密度没有上升到-10附近的情况下，分布看起来也远不是正态的。
• 有几个功能看起来几乎是多模态的。在预处理步骤中，这可以成为一个有用的考虑因素。

与基因特征类似，我们再一次推导并绘制元特征的分布：

cell_stats <- train %>% 
  select(starts_with("c-")) %>% 
  summarise(across(everything(), list(min = min, max = max, mean = mean, sd = sd))) %>% 
  pivot_longer(starts_with("c-"), names_to = "features", values_to = "values") %>% 
  separate(features, into = c("features", "stat"), sep = "_")

cell_stats %>% 
  ggplot(aes(values, fill = stat)) +
  geom_density() +
  scale_fill_manual(values = wes_palette("GrandBudapest1")) +
  facet_wrap(~ stat, scales = "free") +
  theme_tufte() +
  theme(legend.position = "none") +
  labs(x = "", y = "", fill = "State", title = "Cell distribution meta statistics")

我们发现:

• 这与基因元分布形成了鲜明的对比。最明显的是，最小值几乎都在-9.5以下，一直上升到-10的边界。最大值显示了3和5之间更广泛的分布。
• 这种不平衡的结果是，平均值向负值偏移，约为-0.5。注意，所有的均值都不大于0。与基因数据相比，标准差的分布从1左右移动到2左右，有一个明显的偏向小值的尾巴。

目标：

所有的靶点都是二元列，表明某种细胞类型是否对药物有反应。一些目标类别还测量响应的类型；例如，有肾上腺素受体激动剂类和肾上腺素受体拮抗剂类。对于同一样例行，这些可能不应该是活动的。

我们的挑战是一个多标签分类问题，因此行（即药物样本）可以有多个MoA（即可以激活多个目标类别）。让我们首先看一下一次可以激活的目标类别的分布。

注意，这里我使用了新的dplyr rowwise语法和哲学，这使得将标准dplyr动词应用于行变得容易得多。扩展代码单元，看看它是如何工作的。

rowstats <- targets %>% 
  select(-sig_id) %>% 
  rowwise() %>% 
  mutate(sum = sum(c_across(everything()))) %>% 
  select(sum) %>% 
  ungroup()

rowstats %>% 
  count(sum) %>% 
  add_tally(n, name = "total") %>% 
  mutate(perc = n/total) %>% 
  mutate(sum = as.factor(sum)) %>% 
  ggplot(aes(sum, n, fill = sum)) +
  geom_col() +
  geom_text(aes(label = sprintf("%.2f%%", perc*100)), nudge_y = 500) +
  # scale_y_continuous(labels = scales::percent) +
  scale_fill_brewer(palette = "Set2") +
  theme_minimal() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", title = "Number of Activations per Sample")

我们发现:

大约39%的训练样本根本没有MoA注释（例如，所有目标类的值都为0）。这在一定程度上解释了稀疏目标dataframe，如果其中40%是完全空的。

最大的组（略多于50%的样本）正好有1个MoA注释（例如，一个类别=该行中的一个值“1”）。

对于超过1个MoA注释，我们看到一条尾巴延伸到7个同时MoA（0.03%的病例）。只有2个MoAs略高于5%，3个MoAs略高于1%。其他的情况就更罕见了。注意，没有6个moa的实例。

从不同的角度来看目标：在我们的训练数据中，哪些类别的moa实例最多？即，哪些列在0中有最多的1 ？为此，我们只需要对所有列求和。我们绘制结果分布，并查看顶部和底部类别。注意密度图上的对数x轴：

target_sums <- targets %>% 
  select(-sig_id) %>% 
  summarise(across(everything(), sum)) %>% 
  pivot_longer(everything(), names_to = "target", values_to = "sum")

p1 <- target_sums %>% 
  ggplot(aes(sum)) +
  geom_density(fill = "darkorange") +
  geom_vline(xintercept = 40, linetype = 2) +
  scale_x_log10() +
  theme_minimal() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", title = "MoA count per target class", subtitle = "Dashed line: 40")

p2 <- target_sums %>% 
  arrange(desc(sum)) %>% 
  head(5) %>% 
  mutate(target = str_replace_all(target, "_", " ")) %>% 
  ggplot(aes(reorder(target, sum, FUN = min), sum, fill = sum)) +
  geom_col() +
  coord_flip() +
  scale_fill_gradient(low = "blue1", high = "blue4") +
  scale_x_discrete(labels = function(x) lapply(str_wrap(x, width = 25), paste, collapse="\n")) +
  theme_minimal() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", title = "Classes with most MoAs")

p3 <- target_sums %>% 
  arrange(sum) %>% 
  head(5) %>%  
  mutate(target = str_replace_all(target, "_", " ")) %>% 
  ggplot(aes(reorder(target, sum, FUN = min), sum, fill = sum)) +
  geom_col() +
  coord_flip() +
  scale_fill_gradient(low = "red4", high = "red1") +
  scale_x_discrete(labels = function(x) lapply(str_wrap(x, width = 25), paste, collapse="\n")) +
  theme_minimal() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", title = "Classes with fewest MoAs")

p1 + (p2/p3)

我们发现:

大多数类在24k行训练样本中有大约10到200个moa。虚线标记了40点的峰值。活动行的最大数目刚刚超过800，最小数目是1。

右下图中两个深红色的病例仅检测到1例MoA阳性。得分最高的课程可以在右上角面板中看到。这两个极端都以抑制剂和拮抗剂为主。

通常情况下，最后这一发现开辟了另一种研究途径。尽管类名对我这样的外行来说意义不大，但我不禁注意到有相当多的“抑制剂”、“拮抗剂”、“激动剂”等。让我们试着提取它们并查看它们的频率。并非所有的类名都整齐地归入这些组，因此我们只查看出现多次的类名。这个想法是使用类的名称并提取它的最后一个单词。下面是得到的频率：

target_sums %>% 
  separate(target, into = c("a", "b", "c", "d", "e", "type"), fill = "left") %>% 
  count(type) %>% 
  add_tally(n, name = "total") %>% 
  mutate(perc = n/total) %>% 
  filter(n > 1) %>% 
  ggplot(aes(reorder(type, n, FUN = min), n, fill = n)) +
  geom_col() +
  geom_text(aes(label = sprintf("%.2f%%", perc*100)), nudge_y = 6) +
  coord_flip() +
  scale_fill_viridis() +
  scale_x_discrete(labels = function(x) lapply(str_wrap(x, width = 25), paste, collapse="\n")) +
  theme_minimal() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", title = "Common final terms in class names")

我们发现:

在类名称中有很多共同点：超过50%是“抑制剂”，“拮抗剂”和“激动剂”各自处于15%左右的可比水平，剩下的三个术语要少见得多。

其余的类（约12%）的名称不属于这些简单模式。

5 多功能交互视觉效果

现在我们对各个特征的行为有了更好的了解，让我们研究它们之间的交互。首先，我们将研究同一组中的不同特征，然后将分析扩展到这些组及其组成元素之间的相互作用。我们将使用与上面相同的特征集顺序。

5.1 特征集合内的交互

5.1.1 处理特征

比较这3个处理特征需要一个平面网格。两个特征横跨水平和垂直的网格轴，第3个定义了每个平面内的绘图：

train %>% 
  group_by(cp_type, cp_dose, cp_time) %>% 
  count() %>% 
  mutate(cp_time = as.factor(cp_time)) %>% 
  ggplot(aes(cp_time, n, fill = cp_time)) +
  geom_col() +
  facet_grid(cp_dose ~ cp_type) +
  scale_fill_manual(values = wes_palette("IsleofDogs1")) +
  theme_minimal() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "Treatment Duration", y = "", fill = "State",
       title = "Treatment Feature Interactions", subtitle = "Horizontal: type, Vertical: dose, Bars/Colour: duration")

我们发现:

• 结果图像与我们在上面看到的整体视图一致。对照治疗在不同剂量和持续时间的情况下同样罕见。
• 然而，一个显著的差异是，与分布均匀得多的D2条相比，D1剂量（对照组和复方）48小时处理的百分比略高。

5.1.2 基因表达特性

这是事情可能变得混乱的地方，在这个组中有772个特征。好的一面是，无论哪种方法在这里有效，都应该适用于较小的细胞活力特征集。

让我们试一试。下面是前200个基因特征之间的相关性。没有标签（显然也没有标题），只有高相关（蓝色）和反相关（红色）的颜色编码：

train %>% 
  select(starts_with("g-")) %>% 
  select(seq(1,200)) %>% 
  cor(use="complete.obs", method = "pearson") %>%
  corrplot(type="lower", tl.col = "black",  diag=FALSE, method = "color",
           outline = FALSE, tl.pos = "n", cl.ratio = 0.05)

我们发现:

• 这就是总体情况。它没有告诉我们太多关于单个特征的信息，但肯定有一些模式可以看到。
• 看起来有一个顺序，特征至少是某种程度上相关或反相关的；我们发现，在总列的25%内，一些特征显示出非常小的相关性。

基于以上概述，仔细研究一组特性就足够了。看：仔细看：

train %>% 
  select(starts_with("g-")) %>% 
  select(seq(1,20)) %>% 
  cor(use="complete.obs", method = "pearson") %>%
  corrplot(type="lower", tl.col = "black",  diag=FALSE, method = "color",
           cl.ratio = 0.1)

我们发现:

• 这里我们有前20个基因特征。某些更强的相关性是显而易见的；例如“g-0”vs“g-8”（反相关），或“g-8”vs“g-17”。
• 相关性不大的特征包括“ g-18 ”， “ g-19 ”，还有“ g-2 ”和“ g-3 ”。

让我们通过散点图来详细了解其中一些配对。在这里，我们选择了4对特征样本，并使用简单的线性拟合将它们绘制在一起。每个图标题给出了特征名称及其皮尔逊相关系数值：

p1 <- train %>% 
  janitor::clean_names() %>% 
  ggplot(aes(g_0, g_8)) +
  geom_point(col = "grey40", size = 0.5) +
  geom_smooth(method = "lm", formula = "y~x", col = "darkred") +
  theme_minimal() +
  labs(title = str_c("g-0 vs g-8: coef = ", sprintf("%.2f", cor(train$`g-0`, train$`g-8`))))

p2 <- train %>% 
  janitor::clean_names() %>% 
  ggplot(aes(g_10, g_17)) +
  geom_point(col = "grey40", size = 0.5) +
  geom_smooth(method = "lm", formula = "y~x", col = "darkblue") +
  theme_minimal() +
  labs(title = str_c("g-10 vs g-17: coef = ", sprintf("%.2f", cor(train$`g-10`, train$`g-17`))))

p3 <- train %>% 
  janitor::clean_names() %>% 
  ggplot(aes(g_0, g_3)) +
  geom_point(col = "grey40", size = 0.5) +
  geom_smooth(method = "lm", formula = "y~x", col = "black") +
  theme_minimal() +
  labs(title = str_c("g-0 vs g-3: coef = ", sprintf("%.2f", cor(train$`g-0`, train$`g-3`))))

p4 <- train %>% 
  janitor::clean_names() %>% 
  ggplot(aes(g_10, g_19)) +
  geom_point(col = "grey40", size = 0.5) +
  geom_smooth(method = "lm", formula = "y~x", col = "black") +
  theme_minimal() +
  labs(title = str_c("g-10 vs g-19: coef = ", sprintf("%.2f", cor(train$`g-10`, train$`g-19`))))


(p1 + p2) / (p3 + p4)

我们发现:

• 这些关系当然不是微不足道的。即使对于接近零相关性的实例，点云的形状也被拉长，并显示出多种模式（例如“g-10”vs“g-19”）。
• 相关特征也不直接：“g-0”vs“g-8”显示了至少两类不同的点，浓度为g-8 == -10。在“g-10”vs“g-17”中，我们看到了与“g-10”vs“g-19”类似的图像，只是略有旋转以达到0.5的相关性。

5.1.3 细胞活力特征

现在我们可以对细胞活力特征做同样的分析。在这里，我们一开始的特征较少（“只有”100个）。首先，概述相关图：

train %>% 
  select(starts_with("c-")) %>% 
  cor(use="complete.obs", method = "pearson") %>%
  corrplot(type="lower", tl.col = "black",  diag=FALSE, method = "color",
           outline = FALSE, tl.pos = "n", cl.ratio = 0.05)

我们发现:

• 至少可以这么说，看起来和上面很不一样。然而，鉴于我们在上面看到的分布向-10倾斜的情况，这种情况并不完全出乎意料。
• 为了进一步了解剩余数据点之间的关系，值得考虑删除这些极端负值。

在这里，仔细研究前10个功能就足够了。这里我们选择了一个显示变量，直接显示系数值，颜色编码与上面的相同：

train %>% 
  select(starts_with("c-")) %>% 
  select(seq(1,10)) %>% 
  cor(use="complete.obs", method = "pearson") %>%
  corrplot(type="lower", tl.col = "black",  diag=FALSE, method = "number",
           cl.ratio = 0.1)

我们发现:

数值在0.75和0.9之间变化。这些是相当强的相关性，但我担心负异常值（如果它们确实是异常值）主导了这幅图。

让我们仔细看看散点图：

p1 <- train %>% 
  janitor::clean_names() %>% 
  ggplot(aes(c_1, c_2)) +
  geom_point(col = "grey40", size = 0.5) +
  geom_smooth(method = "lm", formula = "y~x", col = "darkblue") +
  theme_minimal() +
  labs(title = str_c("c-1 vs c-2: coef = ", sprintf("%.2f", cor(train$`c-1`, train$`c-2`))))

p2 <- train %>% 
  janitor::clean_names() %>% 
  ggplot(aes(c_3, c_4)) +
  geom_point(col = "grey40", size = 0.5) +
  geom_smooth(method = "lm", formula = "y~x", col = "darkblue") +
  theme_minimal() +
  labs(title = str_c("c-3 vs c-4: coef = ", sprintf("%.2f", cor(train$`c-3`, train$`c-4`))))

p3 <- train %>% 
  janitor::clean_names() %>% 
  ggplot(aes(c_5, c_9)) +
  geom_point(col = "grey40", size = 0.5) +
  geom_smooth(method = "lm", formula = "y~x", col = "darkblue") +
  theme_minimal() +
  labs(title = str_c("c-5 vs c-9: coef = ", sprintf("%.2f", cor(train$`c-5`, train$`c-9`))))

p4 <- train %>% 
  janitor::clean_names() %>% 
  ggplot(aes(c_0, c_6)) +
  geom_point(col = "grey40", size = 0.5) +
  geom_smooth(method = "lm", formula = "y~x", col = "darkblue") +
  theme_minimal() +
  labs(title = str_c("c-0 vs c-6: coef = ", sprintf("%.2f", cor(train$`c-0`, train$`c-6`))))

(p1 + p2) / (p3 + p4)

我们发现:

• 好吧，我改正了。-10值在这里确实发挥了一定的作用，但即使没有它们，主点簇中仍然存在大量的相关性。
• 同样，我们可以看到多个子簇的实例以及潜在的异常值。

5.2 特征集之间的交互

在上一节中，我们没有研究众多基因表达或细胞活力列之间的每一个单独的特征相互作用，但我们对相关性的强度有了一个很好的概述，并更详细地研究了几个特定的例子。现在，在最终包含目标特征之前，让我们看看不同组的治疗、基因和细胞特征如何相互作用。

5.2.1 预测特性

下面，我们将研究来自treatment、gene和cell三个不同集合的特征之间的交互效应。

让我们从治疗特征开始分解4个完全随机选择的基因和细胞特征的分布。（真的是吗？）在这里，4个特征垂直排列，而治疗时间定义了水平面。治疗剂量用颜色标记。两种分布重叠的区域是两种颜色的混合：

train %>% 
  janitor::clean_names() %>% 
  select(cp_dose, cp_time, g_525, g_666, c_42, c_22) %>% 
  mutate(cp_time = as.factor(str_c("Duration ", cp_time, "h"))) %>% 
  pivot_longer(starts_with(c("g_", "c_")), names_to = "feature", values_to = "value") %>% 
  ggplot(aes(value, fill = cp_dose)) +
  geom_density(alpha = 0.5) +
  facet_grid(feature ~ cp_time) +
  theme_minimal() +
  theme(legend.position = "top") +
  labs(x = "Feature value", fill = "Dose", title = "Treatment features vs example cell & gene distributions")

我们发现:

• 每个层面的剂量几乎没有差别：颜色几乎完全重叠。
• 同样，相同基因或细胞特征的治疗时间之间的任何偏差最多只能是边际的。

大部分分布没有真正的区别，但单元格特征的尾部呢？让我们仔细看看图5的样式。

train %>% 
  select(cp_time, cp_dose, starts_with("c-")) %>% 
  mutate(cp_time = as.factor(str_c("Duration ", cp_time, "h"))) %>% 
  select(seq(1,5)) %>% 
  pivot_longer(starts_with("c-"), names_to = "feature", values_to = "value") %>% 
  filter(value < -4) %>% 
  ggplot(aes(value, fill = cp_dose)) +
  geom_density(alpha = 0.5) +
  facet_grid(feature ~ cp_time) +
  theme_minimal() +
  theme(legend.position = "top", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", fill = "Dose", title = "Treatment vs cell - zoom in on negative tail")

我们发现:

• 这里有更多的方差。然而，请记住，虽然这个负尾很重要，但总体数字仍然相对较低。
• 我们可以看到不同剂量之间的差异，特别是在较短的治疗时间内，但在治疗时间之间仍然存在较大的差异。例如，比较24小时和48小时的“c-0”持续时间，其中（截断）分布的质量在（截断）x轴范围的另一边。

现在，我们还剩下一个治疗特征没有关注：即cp_type：指定我们是有一个真正的复合药物治疗（trt_cp）还是一个对照扰动（ctl_vehicle）。对照组不应该有MoAs，但可以告诉我们一些关于没有药物存在的特征方差的信息。因此，我们预计control和compound会有一些更显著的差异：

train %>% 
  janitor::clean_names() %>% 
  select(cp_type, g_8, g_525, g_666, c_14, c_42, c_22) %>% 
  pivot_longer(starts_with(c("g_", "c_")), names_to = "feature", values_to = "value") %>% 
  ggplot(aes(value, fill = cp_type)) +
  geom_density(alpha = 0.5) +
  facet_wrap(~ feature)+
  theme_minimal() +
  theme(legend.position = "top") +
  labs(x = "Feature value", fill = "Type", title = "Treatment features vs example cell & gene distributions")

我们发现:

• 好吧，现在我们有进展了。并非所有特征都显示出差异，但“g-8”，特别是“g-525”显示出对照与复方的明显转移分布。
• 目前还不清楚这对预测MoAs有多大帮助，但也许它可以告诉我们控制和类似控制（即边际）效应之间的区别。

5.2.2 预测 VS 目标

这就引出了问题的核心：预测器特征对目标类别的影响。这种逐渐的积累可能显得不必要的扩展，但在目标和预测变量对应分布的背景下研究它们之间的相互作用是很重要的。

我们首先将前面导出的行统计量与预测变量特征（治疗组特征和一些选定的细胞和基因特征）连接起来。

stats_all <- train %>% 
  select(starts_with("cp"), num_range(prefix = "g-", c(8, 525)), num_range(prefix = "c-", c(14, 42))) %>% 
  bind_cols(rowstats)

首先，让我们检查控制扰动。这些不应该有任何MoAs：

stats_all %>% 
  group_by(cp_type, sum) %>% 
  summarise(n = n()) %>% 
  add_tally(n, name = "total") %>%
  ungroup() %>% 
  mutate(perc = n/total) %>% 
  mutate(sum = as.factor(sum)) %>% 
  ggplot(aes(sum, n, fill = sum)) +
  geom_col() +
  geom_text(aes(label = sprintf("%.2f%%", perc*100)), nudge_y = 500) +
  scale_fill_brewer(palette = "Set2") +
  facet_grid(~ cp_type, scales = "free_x", space = "free_x") +
  theme_hc() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", title = "Number of Activations per type")

我们发现:

• 控制扰动确实没有MoA，这是预期的，但很好证实。然而，仍然有大约三分之一的剩余复合处理也没有MoAs。这很好；不是每一种试验药物都会有反应。
• 图的其余部分与我们之前在图7中看到的类似。

接下来，我们将研究连续出现的MoAs数量是否会影响细胞或基因分布。在这里，为了样本量的考虑，我们将具有3个或更多moa的实例合并为“3+”组。我们只关注复合治疗；例如：没有控制行。为了一些视觉上的变化，我们将使用小提琴图（即垂直图，镜像密度图）：

stats_all %>% 
  filter(cp_type == "trt_cp") %>% 
  mutate(cp_time = as.factor(str_c("Duration ", cp_time, "h"))) %>% 
  mutate(sum = if_else(sum >= 3, "3+", as.character(sum))) %>% 
  mutate(sum = as.factor(sum)) %>% 
  pivot_longer(starts_with(c("g-", "c-")), names_to = "feature", values_to = "value") %>% 
  ggplot(aes(sum, value, fill = sum)) +
  # geom_violin(draw_quantiles = c(0.25, 0.75)) +
  geom_violin() +
  facet_grid(feature ~ cp_time) +
  theme_minimal() +
  theme(legend.position = "none") +
  labs(x = "Sum of active MoAs per row", y = "Cell or Gene values", fill = "Rowwise sum of MoAs",
       title = "Selected cell & gene distributions for different counts of MoAs per row",
       subtitle = "Facetted by cell/gene vs treatment duration")

我们发现:

• MoAs为0或1与MoAs为2或更多的病例之间存在显著差异。特别是，2种MoA类型似乎占了细胞特征中大多数负尾的原因(以及一些基因；看看“g-8”)。相比之下，“3+”组看起来更类似于“0”和“1”组（除了样本量更小）。
• “2”和“3+”组的cp_time持续时间有一些细微的差异。“0”和“1”组看起来都非常稳定。

这看起来很有希望。现在我们可以对治疗剂量做一些类似的事情。让我们回到密度图，并为剂量特征选择颜色编码，而facet网格由基因和细胞跨越。这里的想法是真正梳理出存在MoA差异时的剂量变化：

stats_all %>% 
  filter(cp_type == "trt_cp") %>% 
  mutate(sum = if_else(sum >= 3, "3+", as.character(sum))) %>% 
  mutate(sum = as.factor(sum)) %>% 
  pivot_longer(starts_with(c("g-", "c-")), names_to = "feature", values_to = "value") %>% 
  ggplot(aes(value, fill = cp_dose)) +
  geom_density(alpha = 0.5) +
  facet_grid(feature ~ sum) +
  theme_minimal() +
  theme(legend.position = "top") +
  labs(y = "", x = "Cell or Gene values", fill = "Dose",
       title = "Selected cell & gene distributions for different counts of MoAs per row",
       subtitle = "Colour-coded treatment dose")

我们发现:

• 至于剂量差异，似乎没有任何差异。几乎所有分布都显示完全重叠。
• 这些密度图从另一个角度突出了“MoA = 2”组与其他组之间的显著差异。

此外，让我们看看各个目标类别及其与预测器特征的相关性。为了给我们一些统计数据，但也为了避免最极端的情况，我将选择“多巴胺受体拮抗剂”使用图8。

foo <- train %>% 
  select(starts_with("cp"), num_range(prefix = "g-", c(8, 525)), num_range(prefix = "c-", c(14, 42))) %>% 
  bind_cols(rowstats) %>% 
  bind_cols(targets %>% select(dopamine_receptor_antagonist)) %>% 
  filter(cp_type == "trt_cp") %>% 
  pivot_longer(starts_with(c("g-", "c-")), names_to = "feature", values_to = "value") %>% 
  select(-cp_type) %>% 
  mutate(dopamine_receptor_antagonist = as.factor(dopamine_receptor_antagonist))

foo %>% 
  ggplot(aes(value, fill = dopamine_receptor_antagonist)) +
  geom_density(alpha = 0.5) +
  facet_wrap(~ feature) +
  theme_minimal() +
  theme(legend.position = "top") +
  labs(y = "", x = "Cell or Gene values", fill = "Class value",
       title = "Selected cell & gene distributions for specific class",
       subtitle = "Example: dopamine receptor antagonist")

我们发现:

• 这很有趣，原因有很多。首先，“g-”特征的分布密度（和均值）有轻微但显著的变化。虽然很微妙，但总的来说，这些转变可能会变得有用。
• 但大多数情况下，在“c-”特征中有一个提示，即极端负尾的值主要为0。

让我们借用图5的方法来更详细地研究这一点。我们将从图8中获取另一个高分类，并放大负面尾部：

foo <- train %>% 
  select(cp_type, num_range(prefix = "c-", seq(1,4))) %>% 
  bind_cols(rowstats) %>% 
  bind_cols(targets %>% select(dopamine_receptor_antagonist, cyclooxygenase_inhibitor)) %>% 
  filter(cp_type == "trt_cp") %>% 
  pivot_longer(starts_with(c("g-", "c-")), names_to = "feature", values_to = "value") %>% 
  select(-cp_type) %>% 
  filter(value < -1) %>% 
  pivot_longer(c(dopamine_receptor_antagonist, cyclooxygenase_inhibitor), names_to = "class_name", values_to = "class_value") %>% 
  mutate(class_value = as.factor(class_value))

foo %>% 
  ggplot(aes(value, fill = class_value)) +
  geom_density(alpha = 0.5) +
  facet_grid(class_name ~ feature) +
  theme_minimal() +
  theme(legend.position = "top") +
  labs(y = "", x = "Cell values", fill = "Class value",
       title = "First cell feature distributions for specific class",
       subtitle = "Example: dopamine receptor antagonist")

我们发现:

• 实际上，这里的分布非常不同。对于“多巴胺”，几乎所有低于-2.5的内容都与值0有关。
• 对于“环氧合酶”，情况要稍微复杂一些，但是“class value == 1”的分布肯定比0值的分布要窄得多。

6个未得分目标

到目前为止，我们在很大程度上忽略了train_targets_nonscored.csv中提供的未得分目标数据，而不是在一开始就观察它并注意到它的稀疏性。现在让我们花点时间仔细看看这个数据集。

我们将以类似于上面的得分目标的方式开始研究它的属性。至于这个笔记本的其余部分，所有的多图布局都由fantastic patchwork包提供支持：

rowstats_non <- targets_non %>% 
  select(-sig_id) %>% 
  rowwise() %>% 
  mutate(sum = sum(c_across(everything()))) %>% 
  select(sum) %>% 
  ungroup()

target_sums_non <- targets_non %>% 
  select(-sig_id) %>% 
  summarise(across(everything(), sum)) %>% 
  pivot_longer(everything(), names_to = "target", values_to = "sum")

p1 <- rowstats_non %>% 
  count(sum) %>% 
  add_tally(n, name = "total") %>% 
  mutate(perc = n/total) %>% 
  mutate(sum = as.factor(sum)) %>% 
  ggplot(aes(sum, n, fill = sum)) +
  geom_col() +
  geom_text(aes(label = sprintf("%.2f%%", perc*100)), nudge_y = 1000) +
  # scale_y_continuous(labels = scales::percent) +
  scale_fill_brewer(palette = "Set2") +
  theme_tufte() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", title = "Number of Activations per Sample Row")

p2 <- target_sums_non %>% 
  ggplot(aes(sum)) +
  geom_density(fill = "darkorange") +
  geom_vline(xintercept = 6, linetype = 2) +
  scale_x_continuous(trans = "log1p", breaks = c(0, 10, 20, 50, 100)) +
  theme_tufte() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", title = "MoA count per target class", subtitle = "Dashed line: 6")

p3 <- target_sums_non %>% 
  arrange(desc(sum)) %>% 
  head(5) %>% 
  mutate(target = str_replace_all(target, "_", " ")) %>% 
  ggplot(aes(reorder(target, sum, FUN = min), sum, fill = sum)) +
  geom_col() +
  coord_flip() +
  scale_fill_gradient(low = "blue1", high = "blue4") +
  scale_x_discrete(labels = function(x) lapply(str_wrap(x, width = 25), paste, collapse="\n")) +
  theme_tufte() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", title = "Classes with most MoAs")

p1 / (p2 + p3) + plot_annotation(title = 'Non-scored target data')

我们发现:

• 正如我们在上面看到的，这个额外的数据框比我们得分的目标数据更加稀疏：80%的行（即药物分子）根本没有MoAs，相比之下，在得分的目标中有40%。
• 有趣的是，MoA计数保持不变：最大值为7，没有出现6的情况。
• 因此，与图8相比，每个目标类的MoA计数分布显著向左偏移。现在有一些类没有任何MoAs；但对于被得分的目标来说，情况并非如此。
• 顶级moa类的计数在几十而不是几百。所有这些都是“抑制剂”。

那么类名中的类型级别部分呢？下面是与图9等价的图形：

target_sums_non %>% 
  separate(target, into = c("a", "b", "c", "d", "e", "f", "g", "type"), fill = "left") %>% 
  count(type) %>% 
  add_tally(n, name = "total") %>% 
  mutate(perc = n/total) %>% 
  filter(n > 1) %>% 
  ggplot(aes(reorder(type, n, FUN = min), n, fill = n)) +
  geom_col() +
  geom_text(aes(label = sprintf("%.2f%%", perc*100)), nudge_y = 12) +
  coord_flip() +
  scale_fill_viridis(option = "cividis") +
  scale_x_discrete(labels = function(x) lapply(str_wrap(x, width = 25), paste, collapse="\n")) +
  theme_minimal() +
  theme(legend.position = "none", plot.subtitle = element_text(size = 10)) +
  labs(x = "", y = "", title = "Non-scored targets: Common final terms in class names")

我们发现:

• 由于类别数量的因子2更大，我们有一些目标类别类型（即目前类名中的最后一个词）出现超过1次的情况。但是，占主导地位的类型保持不变；其中“抑制剂”的贡献超过50%。
• 拮抗剂仍然比激动剂更常见。还有一些奇怪的拼写变化，例如“activators”而不是“activator”。

7. PCA降维

考虑到基因，特别是细胞特征中显著的相关性，让我们测试一些降维方法，看看我们可以将特征空间降低多少。

这里，我们将重点介绍主成分分析（PCA）。PCA本质上是参数空间的旋转，以便新轴（“主成分”又名PC）是正交的，并与最大方差的方向对齐。我们从基因特征开始，然后再看细胞特征。

7.1 基因特性

X <- train %>% 
  select(starts_with("g-"))

pca <- prcomp(X, center = TRUE, scale. = TRUE)

我们的第一张图将显示前5个PC的每个PC解释的方差量（“碎石图”），以及原始特征对前2个PC贡献的方向和大小（颜色）。PCA视觉效果是通过优秀的factoextra包生成的。

p1 <- fviz_eig(pca, title = "Variance", ncp = 5)
p2 <- fviz_pca_var(pca,
             col.var = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE,
             label = "none",
             title = "Variables"
             )

p1 + p2 + plot_annotation(title = 'PCA on gene features: overview')

我们发现:

• 第一个PC的方差约为25%，其他PC的方差小于5%；这是一个显著的下降。
• 变量图显示了PC1的优势，大多数特征的贡献都接近水平轴。贡献的强度在显示的范围内变化相当平稳。

上面右图中的变量图已经表明，许多变量以类似的方式对前2个变量做出贡献。在这里，我们可以仔细看看这些特定的功能。我展示了每个PC维度的前15个变量。（水平线的虚线表示均匀分布的期望值）：

p1 <- fviz_contrib(pca, choice = "var", axes = 1, top = 15)
p2 <- fviz_contrib(pca, choice = "var", axes = 2, top = 15)

p1 / p2 + plot_annotation(title = 'PCA on gene features: variable contributions')

我们发现:

• 这些图对PC1不是特别有用，因为许多不同的特征具有相似的贡献强度。这些特征名称的顺序也不一定揭示任何有趣的东西。
• 在PC2中，我们有“g-173”领先。然而，请记住，这个成分对总方差的贡献小于5%。

现在我们有了主成分，我们可以查看新的参数空间中的数据分布。为了检测有希望的聚类，我们将根据治疗特征（类型、剂量、持续时间）的值以及MoAs（行）之和对数据点进行颜色编码：

p1 <- fviz_pca_ind(pca, label = "none",
                   habillage = train %>% mutate(cp_type = if_else(cp_type == "trt_cp", "Treatment compound", "Treatment control")) %>% pull(cp_type),
                   # habillage = train$cp_type,
                   alpha.ind = 0.3,
              palette = c("#FFCC00", "black"),
             title = "Treatment type") +
  theme(legend.position = "top")

p2 <- fviz_pca_ind(pca, label = "none",
                   habillage = train$cp_dose,
                   alpha.ind = 0.5,
              palette = wes_palette("Cavalcanti1"),
             title = "Treatment dose") +
  theme(legend.position = "top")
        
p3 <- fviz_pca_ind(pca, label = "none",
                   habillage = train$cp_time,
                   alpha.ind = 1,
              palette = "Accent",
             title = "Treatment duration") +
  theme(legend.position = "top")

                        
p4 <- fviz_pca_ind(pca, label = "none",
                   habillage = rowstats %>% mutate(sum = if_else(sum >= 3, "3+", as.character(sum))) %>% pull(sum),
              palette = "Dark2",
              alpha.ind = 0.7,
             title = "Sum of MoAs") +
  theme(legend.position = "top")

(p1 + p2) / (p3 + p4) + plot_annotation(title = 'PCA on gene features by feature groups')

我们发现:

这里有很多分离和聚类；这很有希望。

对照处理在前两个处理中表现出比实际复合处理更紧密的聚类。考虑到它们的天性，这并不奇怪。

剂量1 vs 2的分布非常吻合。这与我们到目前为止看到的关于这个特性的内容是一致的。

另一方面，处理时间在参数空间的一侧显示出一些分离。类似地，MoAs的行求和表明存在某种聚类，特别是sum = 2 vs 0和1。

7.2单元格特性

现在我们可以对单元格特征做同样的事情。鉴于我们在图14和图15中看到的强相关性，我预计PCA看起来与基因特征发现有很大不同。首先概述和可变贡献：

Xc <- train %>% 
  select(starts_with("c-"))

pca_cell <- prcomp(Xc, center = TRUE, scale. = TRUE)

p1 <- fviz_eig(pca_cell, title = "Variance", ncp = 5)
p2 <- fviz_pca_var(pca_cell,
             col.var = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE,
             label = "none",
             title = "Variables"
             )
p3 <- fviz_contrib(pca_cell, choice = "var", axes = 1, top = 8)
p4 <- fviz_contrib(pca_cell, choice = "var", axes = 2, top = 8)

(p1 + p2) / (p3 + p4) + plot_annotation(title = 'PCA on cell features')

我们发现:

• 是的，这些高相关性在第一PC的统治中表现得很明显。
• 贡献图类似于基因特征，因为许多变量对PC1有很大贡献（以几乎一致的方式），而对PC2的贡献强度逐渐下降。

现在我们来看一下新PC空间中的散点图。我将保留之前布局的配色方案，所以不要感到困惑：

p1 <- fviz_pca_ind(pca_cell, label = "none",
                   habillage = train %>% mutate(cp_type = if_else(cp_type == "trt_cp", "Treatment compound", "Treatment control")) %>% pull(cp_type),
                   # habillage = train$cp_type,
                   alpha.ind = 0.3,
              palette = c("#FFCC00", "black"),
             title = "Treatment type") +
  theme(legend.position = "top")

p2 <- fviz_pca_ind(pca_cell, label = "none",
                   habillage = train$cp_dose,
                   alpha.ind = 0.5,
              palette = wes_palette("Cavalcanti1"),
             title = "Treatment dose") +
  theme(legend.position = "top")
        
p3 <- fviz_pca_ind(pca_cell, label = "none",
                   habillage = train$cp_time,
                   alpha.ind = 1,
              palette = "Accent",
             title = "Treatment duration") +
  theme(legend.position = "top")

                        
p4 <- fviz_pca_ind(pca_cell, label = "none",
                   habillage = rowstats %>% mutate(sum = if_else(sum >= 3, "3+", as.character(sum))) %>% pull(sum),
              palette = "Dark2",
              alpha.ind = 0.7,
             title = "Sum of MoAs") +
  theme(legend.position = "top")

(p1 + p2) / (p3 + p4) + plot_annotation(title = 'PCA on cell features by feature groups')

我们发现:

• 特征分布与我们在基因变量中看到的非常相似。复合处理较对照延长得多，D1和D2重叠紧密，在处理时间和MoA总和上分层更明显。
• 就持续时间而言，24小时组尤其突出。48h和72h组在PC空间的边缘划分得稍好一些，但在其他方面重叠得也很紧密。
• 对于MoA和，sum = 2组与其他组最不同。

8 基线模型

在探索的最后一步中，让我们使用我们获得的一些见解来构建一个基准模型，以使我们开始制作排行榜。这个模型本身没有竞争力，但可以很容易地调整和扩展。

到目前为止，神经网络（NN）模型似乎在这场比赛中表现最好。因此，我们将使用Keras/Tensorflow构建我们自己的入门NN。

8.1 Keras + tidymodels

对于主模型，我们将在tidymodels生态系统中进行预处理，特别是great recipes包，然后将预处理后的数据输入Keras。稍后，我计划添加一个额外的部分，直接在Keras中进行预处理。

在这里使用tidymodels的主要原因是，它是一个可以与scikit-learn相媲美的完整的建模生态系统：一致的语法和结构使其易于实验不同的模型，并在相对较少的代码行中构建高效的ML工作流。(注意，tidymodels——流行的caret包的继承者——仍然在积极开发中，到目前为止仍然缺少一些高级功能。然而，基本的ML工作流，特别是recipes包是相当稳定的。)

8.1.1预处理

首先，我们需要做一些预处理以获得正确的数据形状。神经网络是“挑剔”的模型之一，因为它们只能处理没有缺失值的数值输入（除了简单的决策树或catboost）。此外，如果输入是标准化的，它们基于梯度下降的学习风格效果最好。

在这里，我们将使用一个简单的训练/验证分割。要了解更复杂的k折方法（和一些不同的预处理），请查看kxx的这个很棒的内核。

分割将使用tidymodels/rsample函数initial_split。为了方便起见，我们将目标与训练数据连接起来，然后将它们分开，然后再次分开。你想怎么分都行。

分裂由cp_type分层，cp_type当然具有最大的不平衡性和治疗特征的影响。考虑到cp_type == ctl_vehicle没有（也不应该）任何MoAs，就值得考虑是否在建模中包含这些实例。

此外，请注意，目前Keras模型不接受tidyverse tibbles作为输入，只接受矩阵。相应的错误消息并不能很清楚地说明这一点。因此，我们将立即将目标转换为矩阵，并在预处理步骤之后对预测器进行转换。

training <- train %>%
    left_join(targets %>% rename_with(.fn = ~paste0("target_", .), .cols = -sig_id),
              by = "sig_id")

set.seed(4321)
tt_split <- initial_split(training, prop = 0.8, strata = cp_type)

X_train_pre <- training(tt_split) %>% 
    select(-starts_with("target"), -sig_id)
y_train <- training(tt_split) %>% 
    select(starts_with("target")) %>% 
    as.matrix()

X_valid_pre <- testing(tt_split) %>% 
    select(-starts_with("target"), -sig_id)
y_valid <- testing(tt_split) %>% 
    select(starts_with("target")) %>% 
    as.matrix()

rm(training)

现在来定义预处理方法。这包括以下步骤。

• 将cp_type和cp_dose特征编码为整数。因为这两个类别型都只有2层，所以这相当于二进制编码，在这里效果很好。
• cp_time特性的规范化。它描述了一个持续时间，所以我们可以将它视为一个数值特征。
• PCA可以作为本节的另一个预处理步骤。这里我们将其分别应用于基因和细胞特征。我们保持95%的方差。

cl_rec <- X_train_pre %>% 
  recipe() %>% 
  update_role(everything(), new_role = "predictor") %>% 
  step_integer(c(cp_type, cp_dose), zero_based = TRUE) %>%
  step_normalize(cp_time) %>%
  step_pca(starts_with("g-"), threshold = 0.95, prefix = "pcg_") %>%
  step_pca(starts_with("c-"), threshold = 0.95, prefix = "pcc_")

所有预处理转换都只在训练数据上计算。到目前为止，这只是对预处理步骤的说明。为了将它们应用于训练数据，我们首先准备配方，然后使用juice命令。请记住，keras需要矩阵作为输入。

现在，同样的预处理步骤可以应用于验证数据和测试数据，再次准备食谱，并在相应的输入数据集上使用bake方法。这里有不少烹饪的比喻。（我们可以将prep包含在初始配方块中，但这种方法与更通用的tidymodels方法一起工作得更好。）

X_train <- cl_rec %>% prep() %>% juice() %>% as.matrix()
X_valid = cl_rec %>% prep() %>% bake(X_valid_pre) %>% as.matrix()
X_test = cl_rec %>% prep() %>% bake(test) %>% as.matrix()

PCA帮助我们减少了相当大一部分的输入特征：

glue("Number of columns reduced from { ncol(train) } to { ncol(X_train) }")

## Number of columns reduced from 876 to 550

8.1.2模型

首先，让我们使用一个简单的模型：只有两层和一个退出正则化。你会发现它的语法与Python的Keras非常相似。

model <- keras_model_sequential() %>%
  layer_dense(units = 2048, activation = "relu", input_shape = ncol(X_train)) %>% 
  layer_dense(units = 1024, activation = "relu") %>% 
  layer_dropout(0.2) %>%
  layer_dense(units = ncol(y_train), activation = "sigmoid")

model

## Model
## Model: "sequential"
## ________________________________________________________________________________
## Layer (type)                        Output Shape                    Param #     
## ================================================================================
## dense (Dense)                       (None, 2048)                    1128448     
## ________________________________________________________________________________
## dense_1 (Dense)                     (None, 1024)                    2098176     
## ________________________________________________________________________________
## dropout (Dropout)                   (None, 1024)                    0           
## ________________________________________________________________________________
## dense_2 (Dense)                     (None, 206)                     211150      
## ================================================================================
## Total params: 3,437,774
## Trainable params: 3,437,774
## Non-trainable params: 0
## ________________________________________________________________________________

然后使用log-loss和Adam优化器编译模型，并为其选择自定义学习率。如果验证损失在连续5个epoch中没有增加，我们还使用提前停止回调函数来缩短训练时间。注意，restore_best_weights选项默认设置为FALSE。

model %>% compile(
  optimizer = optimizer_adam(lr = 1e-4),
  loss = "binary_crossentropy"
)

cb_early_stopping <- callback_early_stopping(patience = 5, restore_best_weights = TRUE)

现在我们可以训练模型了。我们提供训练和验证数据，以及纪元数和批大小。由于我们使用了提前停止，因此可以选择相对较多的epoch。

由于Jupyter的问题，fit将不会打印建模进度（既不是作为进度条，也不是作为每个epoch的一行），无论详细设置。相反，我们打印训练时长。对于这个简单的模型，我们可以在几分钟内在CPU上进行训练：

tic()
history <- model %>% fit(X_train, y_train,
                         epochs = 20,
                         verbose = 0,
                         batch_size = 32,
                         callbacks = list(cb_early_stopping),
                         validation_data = list(X_valid, y_valid)
                        )
toc()

## 159.851 sec elapsed

当涉及到在提前停止的情况下绘制训练历史时，在撰写本文时，有另一个bug目前使事情变得复杂。幸运的是，这个方法可以在这个代码块中实现，这个代码块后面有一个非常有用的github注释。

#https://github.com/rstudio/keras/issues/1116

# set plot.keras_training_history in your global env to call user functions
plot.keras_training_history <- keras:::plot.keras_training_history
environment(plot.keras_training_history) <- globalenv()

# replace as.dataframe with custom function
as.data.frame.keras_training_history <- function (x, ...){

 if (tensorflow::tf_version() < "2.2") 
  x$metrics <- x$metrics[x$params$metrics]
 values <- x$metrics

 # pad <- x$params$epochs - length(values$loss)
 # pad_data <- list()
 # for (metric in x$params$metrics) pad_data[[metric]] <- rep_len(NA, 
 #                                                                pad)
 # values <- rbind(values, pad_data)

 values[] <- lapply(values, `length<-`, x$params$epochs)
 df <- data.frame(epoch = seq_len(x$params$epochs), value = unlist(values), 
                  metric = rep(sub("^val_", "", names(x$metrics)), each = x$params$epochs),
                  data = rep(grepl("^val_", names(x$metrics)), each = x$params$epochs))
 rownames(df) <- NULL
 df$data <- factor(df$data, c(FALSE, TRUE), c("training", "validation"))
 df$metric <- factor(df$metric, unique(sub("^val_", "", names(x$metrics))))
 df
 
}

我们仍然需要手动调整x轴范围，但至少我们可以比较训练和验证损失历史：

plot(history) +
  coord_cartesian(xlim = c(0, 15), ylim = c(0.01, 0.02)) +
  theme_minimal() +
  labs(title = "Learning curves")

8.1.3预测和完整性检查

现在我们可以使用该模型在测试集上进行预测。请记住，测试数据的预处理方法与训练数据相同。我们还做了一点预处理，将所有控制样本的概率设置为零。然后我们创建提交文件。

pred <- model$predict(X_test)
colnames(pred) <- str_remove_all(colnames(y_valid), "target_")
pred <- pred %>%
    as_tibble()

pred[test$cp_type == "ctl_vehicle",] <- 0

submit <- test %>%
    select(sig_id) %>%
    bind_cols(as_tibble(pred))

这就是我们的提交文件。

在这里，我们添加一些快速的完整性检查，以确保我们的提交文件具有正确的形状、列名和id。这是为了避免因为一个小bug而丢失提交的内容。这三种身份检查的结果都应该是“TRUE”。

c(identical(dim(submit), dim(sample_submit)),
  identical(submit$sig_id, sample_submit$sig_id),
  identical(colnames(submit), colnames(sample_submit)))

## [1] TRUE TRUE TRUE

然后我们写入submit .csv输出。

submit %>%
    write_csv("submission.csv")

9 还有一件事： drug id

为了完整起见，在比赛进入最后几天时，这里是药物id的附加文件（train_drug.csv）的概述，它在11月初添加到正式比赛数据集。

drugs <- vroom(str_c(path, "train_drug.csv"), col_types = cols())

在这个文件中，训练集的每一行都有一个匿名后的drug_id（但测试集没有）：

drugs %>% 
  head(100) %>% 
  DT::datatable()

对于23814行训练数据，有3289个唯一的药物id。让我们看看最频繁的那些以及总体频率分布。

注意，这里我们对密度分布的y轴使用了平方根缩放。当然，这扭曲了对曲线下（橙色）区域的解释，但这种扭曲在这里没有突出低密度实例重要，特别是那些在（对数尺度）x轴较长端。

p1 <- drugs %>% 
  count(drug_id) %>% 
  ggplot(aes(n, fill = "lightorange")) + 
  geom_density() +
  scale_x_continuous(trans = "log10", breaks = c(1, 6, 10, 100, 1000)) +
  scale_y_sqrt() +
  theme_minimal() +
  theme(legend.position = "none", title = element_text(size = 9)) +
  labs(x = "", title = "Distribution of row counts")

p2 <- drugs %>% 
  count(drug_id) %>% 
  add_tally(n, name = "total") %>% 
  mutate(perc = n/total) %>% 
  slice_max(order_by = n, n = 10) %>% 
  ggplot(aes(reorder(drug_id, n, FUN = min), n, fill = n)) +
  geom_col() +
  geom_text(aes(label = sprintf("%.1f%%", perc*100)), nudge_y = 250) +
  coord_flip(ylim = c(0, 2200)) +
  scale_fill_viridis(option = "cividis") +
  scale_x_discrete(labels = function(x) lapply(str_wrap(x, width = 25), paste, collapse="\n")) +
  theme_minimal() +
  theme(legend.position = "none", title = element_text(size = 9)) +
  labs(x = "", y = "", title = "Top 10 drugs")

p1 + p2 + plot_annotation(title = 'Drug IDs')

我们发现:

• 密度分布在每种药物训练6行时达到峰值。还有一些明显的峰值在1以上，超过10，还有一些超过100。
• 超过100的实例表示最常见的药物id；如条形图所示。最常见的是一个平淡的名字“cacb2b860”，几乎占所有行的8%。其次，出现频率第二高的ID是“87d714366”——这个标签在描述上同样有价值。
• 值得注意的是，有一组7种药物的贡献率为0.7% - 1.0%。其他所有行只占0.1%或更少。

现在，大约8%的行听起来像是一个熟悉的数字。如果你一直到这本笔记本的开头你会发现1866的实际行数。是的，这些是对照治疗。因此，“cacb2b860”应该对应于cp_type ==“ctl_vehicle”。我们来确保这是正确的。下表应该只包含1行ctl_vehicle的1866条记录：

有了控制行，我们可以专注于“真正的”药物。这种密度分布看起来相当可疑——它的多个尖峰和普遍缺乏平滑。让我们选择一个稍微不同的角度：我们再次计算每个drug_id的行数，然后计算这些行数。（例如，有3种药物的数量正好是2个实例。）然后我们将这些行计数作为分类特征，并通过barplot显示药物计数。

再次注意y轴上的平方根缩放。通常，不建议对条形图应用这种转换（因为读者应该比较条形高度来判断它们的差异），但这里的重点（再次）是在频繁情况下突出罕见情况。有时候，你必须打破规则：

drugs %>% 
  filter(drug_id != "cacb2b860") %>% 
  count(drug_id, name = "rows") %>% 
  count(rows) %>% 
  add_tally(n, name = "total") %>% 
  mutate(perc = n/total) %>% 
  ggplot(aes(as.factor(rows), perc, fill = n)) +
  geom_col() +
  geom_text(aes(label = sprintf("%i", n)), nudge_y = 0.05) +
  scale_fill_viridis(option = "cividis") +
  scale_x_discrete(labels = function(x) lapply(str_wrap(x, width = 25), paste, collapse="\n")) +
  scale_y_sqrt(labels = scales::percent) +
  theme_minimal() +
  theme(legend.position = "none") +
  labs(x = "", y = "", title = "Drug IDs: Frequency of row counts", subtitle = "Categorical scaling on the x-axis.")

我们发现:

• 确实存在相对较少的distt drug_id计数。此外，只有5种变体出现次数超过10次：绝大多数（> 75%）drug_ids正好出现在6行。第二名分别是7排和1排。另外，5行和12行频率几乎相同。
• 所有出现的数量超过100行的数据都是唯一的。对应于上图33中排名靠前的药物。

通过这个额外的分析，我结束了我的探索性数据分析。当然，人们可以用新的drug_id数据集做更多的事情，但在这里我只做一个概述。

非常感谢您的阅读！我感谢所有留下反馈、鼓励或赞的人。感谢您的支持！