个性化医疗 - 使用tidy R进行探索性分析

1介绍

这是针对个性化医疗：重新定义癌症治疗挑战的全面探索性数据分析。我将使用ggplot2和tidyverse工具来研究和可视化数据中的结构。

我们面临的挑战是，在不影响肿瘤的突变（“驱动基因”）存在的情况下，自动对有助于癌症肿瘤生长的基因突变（所谓的“驱动”）进行分类。

数据来自4个不同的文件。两个csv文件和两个文本文件：

训练/测试变体：这些是基因突变的csv目录以及目标值类别，目标值类别是对突变的（手动）分类评估。特征变量是Gene（发生突变的特定基因）和Variation（突变的性质）。测试数据当然没有类别值。这是我们要预测的。这两个文件分别通过ID变量链接到另一个文件，即：

训练/测试文本：这些文本包含（专家）用于手动标记突变类的证据的广泛描述。

文本信息是分类问题的关键，必须被很好地理解/建模才能达到有用的准确性。

1.1 加载库和数据文件

library('ggplot2') # visualization
library('ggthemes') # visualization
library('scales') # visualization
library('grid') # visualisation
library('gridExtra') # visualisation
library('corrplot') # visualisation
library('ggraph') # visualisation
library('igraph') # visualisation
library('dplyr') # data manipulation
library('readr') # data input
library('tibble') # data wrangling
library('tidyr') # data wrangling
library('stringr') # string manipulation
library('forcats') # factor manipulation
library('tidytext') # text mining
library('SnowballC') # text analysis
library('wordcloud') # test visualisation

读取变体数据表：

train <- read_csv('../input/training_variants')
test  <- read_csv('../input/test_variants')

目前在R中读取文本文件需要花费更多的精力。

train_txt_dump <- tibble(text = read_lines('../input/training_text', skip = 1))
train_txt <- train_txt_dump %>%
  separate(text, into = c("ID", "txt"), sep = "\|\|")
train_txt <- train_txt %>%
  mutate(ID = as.integer(ID))

test_txt_dump <- tibble(text = read_lines('../input/test_text', skip = 1))
test_txt <- test_txt_dump %>%
  separate(text, into = c("ID", "txt"), sep = "\|\|")
test_txt <- test_txt %>%
  mutate(ID = as.integer(ID))

2 变体数据表

在这个EDA中，我们首先看一下变体数据文件，这些文件通过常见的可视化工具更容易访问。

2.1 第一表数据概述：

我们使用summary和glimpse工具进行第一次概述：

train <- train %>%
  mutate(Gene = factor(Gene),
         Variation = factor(Variation),
         Class = factor(Class))

test <- test %>%
  mutate(Gene = factor(Gene),
         Variation = factor(Variation))

summary(train, maxsum = 9)

glimpse(train)

nrow(train)

nrow(test)

sum(is.na(train))

train %>%
  group_by(Gene) %>%
  summarise(ct = n()) %>%
  arrange(desc(ct))

test %>%
  group_by(Gene) %>%
  summarise(ct = n()) %>%
  arrange(desc(ct))

train %>%
  group_by(Variation) %>%
  summarise(ct = n()) %>%
  arrange(desc(ct))

test %>%
  group_by(Variation) %>%
  summarise(ct = n()) %>%
  arrange(desc(ct))

我们发现:

在训练集中有3321个不同的id，包含264个不同的基因表达，2996个不同的变异。用整数水平表示有9个不同的类别。

基因和变异特征包含不同长度的字符串。

测试数据比训练数据多70%。数据描述告诉我们，“一些测试数据是机器生成的，以防止手工标记。，这应该可以解释这种奇怪的失衡。

变体数据中没有缺失值。

训练数据和测试数据中最频繁的基因完全不同。此外，测试数据似乎包含比训练数据更多的不同基因和更少的高频基因。在某种程度上，这可能是测试数据中添加的机器生成条目的影响（通过添加许多不同的随机水平）。因此，频率上的差异可能反映了有效测试数据相对于训练数据的真实比例。

相比之下，训练集和测试集最频繁的变化基本相同；尽管，同样，测试数据中对应的频率较低（为5 - 10的倍数）。

2.2个体特征可视化

这是最频繁的基因值的频率分布：

top_gene <- train %>%
  group_by(Gene) %>%
  summarise(ct = n()) %>%
  filter(ct > 40)

top_gene %>%
  ggplot(aes(reorder(Gene, -ct, FUN = min), ct)) +
  geom_point(size = 4) +
  labs(x = "Gene", y = "Frequency") +
  coord_flip()

top_gene_test <- test %>%
  group_by(Gene) %>%
  summarise(ct = n()) %>%
  filter(ct > 40)

top_gene_test %>%
  ggplot(aes(reorder(Gene, -ct, FUN = min), ct)) +
  geom_point(size = 4) +
  labs(x = "Gene", y = "Frequency") +
  coord_flip()

我们发现:

• 在训练数据和测试数据中，相对较小的基因水平组构成了相当大的一部分特征值。
• 测试数据中高频基因较少。

这些是训练数据（蓝色）与测试数据（红色）中最常见的变化；通过比较表中的数据来确认我们已经看到的：

foo <- train %>% mutate(set = factor("train")) %>% select(-Class, -ID)
bar <- test %>% mutate(set = factor("test")) %>% select(-ID)

foo <- full_join(foo, bar)

foo %>%
  group_by(Variation, set) %>%
  summarise(ct = n()) %>%
  filter(ct > 3) %>%
  ggplot(aes(reorder(Variation, -ct, FUN = median), ct, colour = set)) +
  geom_point(size = 4) +
  coord_cartesian(ylim = c(0, 100)) +
  labs(x = "Variation", y = "Frequency")

下面我们看一下类别目标在训练数据中的分布：

train %>%
  ggplot(aes(Class)) +
  geom_bar()

我们发现:

• 3级、8级和9级的学生人数明显不足
• 级别5和6是可比的，中低频
• 级别1、2和4是可比的，中高频率
• 7级显然是最常见的

2.3特征交互

现在我们来研究特征之间以及特征与目标类别变量之间的相互作用。

2.3.1基因vs类

首先，我们将查看不同类别的总体最频繁基因的频率分布。注意对数频率标度。

train %>%
  filter(Gene %in% str_c(top_gene$Gene)) %>%
  ggplot(aes(Gene)) +
  geom_bar() +
  scale_y_log10() +
  theme(axis.text.x  = element_text(angle=90, vjust=0.5, size=7)) +
  facet_wrap(~ Class)

我们立即看到了显著的差异：

• 一些基因，比如“PTEN”，主要存在于一个类别中（这里是4）。
• 其他基因，如“TP53”，主要在两类之间共享（这里：1和4）。
• 第8类和第9类不包含最常见的基因。

下面是按基因排序的类的样子（同样是日志计数）：

train %>%
  filter(Gene %in% str_c(top_gene$Gene)) %>%
  ggplot(aes(Class)) +
  geom_bar() +
  scale_y_log10() +
  facet_wrap(~ Gene)

这种表示强调了我们关于不同类别中相似/主导基因的发现。

2.3.2基因vs变异

接下来，我们在某种程度上重新利用计数图来可视化最常见基因的变异分布。由于有很多不同的变异，我们放弃y轴标签，只是说明数据中存在多少基因变异组合。

首先是训练数据：

foo <- train %>%
  filter(Gene %in% str_c(top_gene$Gene)) %>%
  group_by(Gene, Variation) %>%
  summarise(ct = n())

y_labels <- str_sub(foo$Variation, start = 1, end = 5)
  
foo %>%
  ggplot(aes(reorder(Gene, ct, FUN = median), reorder(Variation, ct, FUN = median))) +
  geom_count() +
  labs(x = "Gene", y = "Variation") +
  theme(axis.text.x  = element_text(angle=90, vjust=0.5, size=7),
        axis.ticks = element_blank(), axis.text.y = element_blank(),
        legend.position = "none")

然后是测试数据：

foo <- test %>%
  filter(Gene %in% str_c(top_gene$Gene)) %>%
  group_by(Gene, Variation) %>%
  summarise(ct = n())

y_labels <- str_sub(foo$Variation, start = 1, end = 5)
  
foo %>%
  ggplot(aes(reorder(Gene, ct, FUN = median), reorder(Variation, ct, FUN = median))) +
  geom_count() +
  labs(x = "Gene", y = "Variation") +
  theme(axis.text.x  = element_text(angle=90, vjust=0.5, size=7),
        axis.ticks = element_blank(), axis.text.y = element_blank(),
        legend.position = "none")

同样，这两个数据集在这个视图中是相当异构的。

3文本文件

3.1概述

第二种数据文件包含大量看起来像科学论文或会议记录的文本。下面是第一个条目的开头：

str_sub(train_txt$txt[1], start = 1, end = 1e3)

当然，我们可以很容易地确认完整条目的第一部分对应于本文，然后切换到这一部分（以及其他相关的部分）。因此，该数据文件似乎是分类所基于的论文的完整出版文本的数据转储（包括图表说明、手稿结构，有时还包括从属关系）。

我怀疑，在确定哪些关键字重要、哪些不重要方面，一些领域知识将大有帮助。这将是一个有趣的练习，看看信息在科学出版物中传达得有多清楚。

3.2数据清洗和准备

在这里，我想收集我在最初探索中注意到的各种文本特征、工件和全局属性。该列表可能会随着内核的增长而扩展。

科学术语和停用词：大多数科学论文都有一种共同的语言风格，在整个文本文件中合理地保持一致。像“结果”或“讨论”这样的词将频繁出现，但不一定包含我们的预测目标的任何信号。因此，下面我定义了我自己的额外停用词列表。

研究领域相关的停用词：我的印象是，停用词列表可以通过包括整个研究领域的特征术语来扩展，这些特征术语是如此普遍，以至于它们的高频率可能会掩盖真正有趣的术语。诸如“突变”、“癌症”或“肿瘤”之类的词似乎太笼统了，在这里没有太多的区分能力。下面的TF-IDF似乎证实了这一点。从具有领域知识的人那里得到一些反馈，告诉他们哪些其他术语可以先验地从文本中删除，这将是很有趣的。

纸张表示法的怪癖：将纸张文本直接转换为ascii会导致许多人工操作。单独考虑这些因素都不会有很大的影响，但如果放在一起，它们可能会降低分析的准确性：

• 引用编号（如自然杂志使用）附加在相应的词
• 有时，在文本中嵌入类似“SectionNext”的网页导航命令
• 作者的名字和从属关系偶尔也包括在内

3.3特征工程

这里我们收集了一些特征工程的想法。

3.3.1文本长度- txt_len

train_txt <- train_txt %>%
  mutate(txt_len = str_length(txt),
         set = "train")

test_txt <- test_txt %>%
  mutate(txt_len = str_length(txt),
         set = "test")

combine_txt <- full_join(train_txt,test_txt)

作为早期探索，我们可以看看文本特征长度的分布。先验地说，我不希望论文的长度与分类结果有关；但可能有些分类只需要一篇论文，而其他分类则需要检查多个论文。

首先，这是训练集和测试集文本输入长度的总体分布：

combine_txt %>%
  ggplot(aes(txt_len, fill = set)) +
#  geom_density(alpha = 0.5, bw = 5e3) +
  geom_histogram(bins = 50) +
  labs(x = "Length of text entry")

分布形状的差异可能也是由于添加到测试样本的机器生成的条目。

现在，让我们看看这个分布是否会随着不同的目标类别而改变。首先，进行面片包装比较：

foo <- train_txt %>%
  select(ID, txt_len)
bar <- train %>%
  select(ID, Class)

full_join(foo, bar, by = "ID") %>%
  ggplot(aes(txt_len)) +
  geom_density(fill = "red", bw = 5e3) +
  labs(x = "Length of text entry") +
  facet_wrap(~ Class)

然后叠加经验累积密度函数：

foo <- train_txt %>%
  select(ID, txt_len)
bar <- train %>%
  select(ID, Class)

full_join(foo, bar, by = "ID") %>%
  ggplot(aes(txt_len)) +
  stat_ecdf(geom = "step") +
  stat_ecdf(aes(txt_len, color = Class), geom = "step") +
  labs(x = "Length of text entry")

以及每个类别长度的中位数：

foo <- train_txt %>%
  select(ID, txt_len)
bar <- train %>%
  select(ID, Class)

full_join(foo, bar, by = "ID") %>%
  group_by(Class) %>%
  summarise(l_med = median(txt_len))

我们发现:

• 在测试长度分布的形状和中位数上存在显著差异。第8类和第9类平均需要更多的文本，而第3类需要的论文最短/最少。
• 考虑到它的价值，很容易推测单个类别的文本长度分布中明显的多个峰值可能对应于构成临床证据的论文数量。

3.3.2文本值缺失

在讨论中，有人指出，一些观察的文本特征中有一个“null”条目。使用txt_len特征，我们可以确认这一发现，并轻松地表明没有其他文本值小于100个字符（只是为了防止使用另一个Null标识符）：

combine_txt %>%
  filter(txt_len < 100)

3.3.3 关键词频率 - 步行法

我想利用这次比赛更多地了解文本挖掘。在深入研究各种工具和技术的应用时，我会记录下我所学到的内容。如果你和我一样是初学者，那么也许这种方法对你有用。如果你是专家，那么请随意跳过所有初级信息（也许可以告诉我如果我犯了严重的错误）。
在使用专业工具之前，这里首先介绍一种基于标准字符串操作方法的初步方法。
分析临床证据内容的显而易见的第一步是查看相应论文的文本中是否频繁提及某些关键词。
我们在本综述文章中使用了“病原体”和“良性”这两个词来命名这5个类别。以下是每个类别中这两个词的出现频率：

train_txt <- train_txt %>%
  mutate(f_pathogenic = str_count(txt, "pathogenic"),
         f_benign = str_count(txt, "benign")
         )

这些是单词“致病性”在我们的9个类别中的频率分布（注意y轴为对数）：

foo <- train_txt %>%
  select(ID, f_benign, f_pathogenic)
bar <- train %>%
  select(ID, Class)

full_join(foo, bar, by = "ID") %>%
  ggplot(aes(f_pathogenic)) +
  geom_bar() +
  scale_y_log10() +
#  scale_x_log10() +
  facet_wrap(~ Class)

在这里，我们绘制了单词“致病”（致病）每个类别的平均出现率与单词“良性”（benign）平均出现率的比率：

foo <- train_txt %>%
  select(ID, f_benign, f_pathogenic)
bar <- train %>%
  select(ID, Class)

full_join(foo, bar, by = "ID") %>%
  group_by(Class) %>%
  summarise(mean_benign = mean(f_benign),
            mean_pathogenic = mean(f_pathogenic),
            path_ben = mean(f_pathogenic)/mean(f_benign)) %>%
  ggplot(aes(reorder(Class, -path_ben, FUN = max), path_ben)) +
  geom_point(colour = "red", size = 3) +
  labs(x = "Class", y = "# occurences 'pathogenic' / # occurences 'benign'")

我们发现:

• facet图显示，单词“致病”在某些类别（如1、4或5）中明显更频繁
• 比率图证实了这一印象，并显示了两个不同的类别组：2、7、8、9与1、3、4。平均而言，后者提到“致病”的比例高于前者。此外，第5类和第6类的“致病”比例甚至高于“良性”。

当然，其中一些出现可能是“非致病”或“非良性”的，这就是为什么我们需要进一步深入文本分析来解决这个谜题。

3.4使用tidytext进行文本分析的第一步

正如tidytext包的作者所说：整洁的文本格式被定义为每行一个标记的表；标记是一个单词或另一个有意义的文本单位（转述）。通过整理文本，我们可以利用tidyverse强大的工具来处理和分析文本文件。我将学习这本优秀的免费在线书籍。

为了使我们的文本数据保持整齐，我们使用unnest_tokens工具。这还可以去掉标点符号，并将所有内容转换为小写：

t1 <- train_txt %>% select(ID, txt) %>% unnest_tokens(word, txt)
head(t1)

tidytext包包含一个停用词字典，如“and”或“next”，我们可以从整洁的文本数据中删除它们。此外，我们将根据科学论文的典型结构化语言来定义我们自己选择的停用词。我们还删除了只有数字或符号的标记。

data("stop_words")
my_stopwords <- data_frame(word = c(as.character(1:100),
                                    "fig", "figure", "et", "al", "table",
                                    "data", "analysis", "analyze", "study",
                                    "method", "result", "conclusion", "author",
                                    "find", "found", "show", "perform",
                                    "demonstrate", "evaluate", "discuss"))
t1 <- t1 %>%
  anti_join(stop_words, by = "word") %>%
  anti_join(my_stopwords, by = "word") %>%
  filter(str_detect(word, "[a-z]"))

作为第一个概述，我们来看一下最流行的单词及其频率。这是我们第一次将tidyverse和ggplot2工具应用于文本数据：

t1 %>%
  count(word) %>%
  filter(n > 5e4) %>%
  mutate(word = reorder(word, n)) %>%
  ggplot(aes(word, n)) +
  geom_col() +
  xlab(NULL) +
  coord_flip()

总的来说，这些词是我们期望在癌症研究和遗传学的出版物中找到的。例如，您将注意到，排名前4的单词实际上是两个基本单词的两个变体。就我们的目的而言，这些单词变体可能会混淆我们感兴趣的信号。我们可以使用词干提取工具将它们简化为基本含义，它们的词干。

据我所知，tidytext目前没有原生的词干提取函数。因此，我们将使用SnowballC包及其wordStem工具：

t1 <- t1 %>%
  mutate(word = wordStem(word))

t1 %>%
  count(word) %>%
  filter(n > 5e4) %>%
  mutate(word = reorder(word, n)) %>%
  ggplot(aes(word, n)) +
  geom_col() +
  xlab(NULL) +
  coord_flip()

结果向我们展示了在整个文本数据中最频繁的基本单词。另一种可视化这些频率的方法是通过词云。就我个人而言，我怀疑词云可能相当于饼图的文本。但知道如何将它们整合到整洁的文本分析中是很有用的：

t1 %>% 
  count(word) %>%
  with(wordcloud(word, n, max.words = 100))

3.5类相关词频

为了使用这些词频进行预测，我们首先需要分别确定每个类别的词频。下面，我们将“文本”数据与“变体”数据集中的类信息连接起来。然后，按类别确定每个词的相对词频。

在本例中，我们将最常见的类别== 7与类别1和类别2进行比较。此外，为了保持概览，我们将只查看每个类别中出现次数超过1000次的单词。在这里，使用dplyr工具的能力开始得到适当的回报：

foo <- train %>%
  select(ID, Class)

t1_class <- full_join(t1, foo, by = "ID")

frequency <-t1_class %>%
  count(Class, word) %>%
  filter(n > 5e2) %>%
  group_by(Class) %>%
  mutate(freq = n / sum(n)) %>% 
  select(-n) %>% 
  spread(Class, freq) %>% 
  gather(Class, freq, `1`:`2`)

然后，为了可视化概述，我们绘制类别7中的单词频率与其他两个类别的对比（注意对数轴）：

ggplot(frequency, aes(x = freq, y = `7`, color = abs(`7` - freq))) +
  geom_abline(color = "gray40", lty = 2) +
  geom_jitter(alpha = 0.1, size = 2.5, width = 0.1, height = 0.1) +
  geom_text(aes(label = word), check_overlap = TRUE, vjust = 1.5) +
  scale_x_log10(labels = percent_format()) +
  scale_y_log10(labels = percent_format()) +
  #scale_color_gradient(limits = c(0, 0.001), low = "darkslategray4", high = "gray95") +
  facet_wrap(~Class, ncol = 2) +
  theme(legend.position="none") +
  labs(y = "Class 7", x = NULL)

在这些图中，接近虚线的单词（频率相等）在相应的类别中具有相似的频率。在特定类别轴线上越靠近的单词（例如类别7和类别1的“抑制剂”）在该类别中出现的频率更高。蓝灰色刻度表示第7类频率与整体频率的不同程度（相对频率越高越轻）。背景中的（轻微抖动的）点代表完整的（高频）单词集合，而显示的单词是为了避免重叠而选择的。

这些图给了我们一个有用的概述。例如，他们认为类别2和7比类别1和7更相似。对于一个更系统的方法，我们计算每个频率集的相关系数（这次是完整的列表，而不是超过1000次的频率）：

frequency <-t1_class %>%
  count(Class, word) %>%
  #filter(n > 1e3) %>%
  group_by(Class) %>%
  mutate(freq = n / sum(n)) %>% 
  select(-n) %>% 
  spread(Class, freq)

frequency %>%
  select(-word) %>%
  cor(use="complete.obs", method="spearman") %>%
  corrplot(type="lower", method="number", diag=FALSE)

我们发现:

• 类别2和7实际上是最相似的类别，其次是1和4（相关系数大于0.9）。

• 总的来说，最不同的类别似乎是第9类，特别是与第3和第5类（它们彼此并不十分相似）相比。让我们看看这些组合的词频扩展是什么样的：

foo <- train %>%
  select(ID, Class)

t1_class <- full_join(t1, foo, by = "ID")

frequency <-t1_class %>%
  count(Class, word) %>%
  filter(n > 2e1) %>%
  group_by(Class) %>%
  mutate(freq = n / sum(n)) %>% 
  select(-n) %>% 
  spread(Class, freq) %>% 
  gather(Class, freq, `3`,`5`)

ggplot(frequency, aes(x = freq, y = `9`, color = abs(`9` - freq))) +
  geom_abline(color = "gray40", lty = 2) +
  geom_jitter(alpha = 0.1, size = 2.5, width = 0.1, height = 0.1) +
  geom_text(aes(label = word), check_overlap = TRUE, vjust = 1.5) +
  scale_x_log10(labels = percent_format()) +
  scale_y_log10(labels = percent_format()) +
  #scale_color_gradient(limits = c(0, 0.001), low = "darkslategray4", high = "gray95") +
  facet_wrap(~Class, ncol = 2) +
  theme(legend.position="none") +
  labs(y = "Class 9", x = NULL)

我们发现:

• 与前一组图相比，这里的散点明显更多；特别是5级和9级。
• 有趣的是，“良性”和“病原”在类别3和类别9中都更常见。

3.6 TF-IDF分析-基础与应用

随着竞争的进行，您可能会在内核和讨论中越来越频繁地看到这种缩写组合。作为一个像我这样的初学者，你可能不知道它意味着什么。让我们从基础开始：

• TF表示词频；也就是一个单词在文本中出现的频率。这是我们在上面测量的。停用词列表可以用来过滤掉那些对我们要回答的问题可能没有影响的频繁单词（例如“and”或“the”）。然而，使用停用词可能并不总是一种优雅的方法。IDF来营救。
• IDF表示逆文档频率。在这里，我们更强调文档集合（在我们的例子中是指整个文本数据）中罕见的单词。
• 这两种度量方法都可以组合成TF-IDF，它是一种启发式索引，告诉我们一个单词在一个更大的文档（所有类）的上下文中（这里是某个类）的频率。你可以将其理解为相对文本频率与整体文档频率的标准化。这将导致突出的单词具有特定类别的特征，这几乎是我们为了训练模型而想要实现的。

Tidytext提供了bind_tf_idf函数，可以从包含单词及其每个类别的数量的整洁数据集中提取这些指标：

frequency <-t1_class %>%
  count(Class, word)

tf_idf <- frequency %>%
  bind_tf_idf(word, Class, n)

让我们可视化最具特征的单词及其类别：

tf_idf %>%
  arrange(desc(tf_idf)) %>%
  mutate(word = factor(word, levels = rev(unique(word)))) %>%
  top_n(20, tf_idf) %>%
  ggplot(aes(word, tf_idf, fill = Class)) +
  geom_col() +
  labs(x = NULL, y = "tf-idf") +
  coord_flip()

我觉得这看起来很有技术含量。对谷歌的快速搜索显示，“dnmt3b7”实际上是“DNA甲基转移酶dnmt3b7的异常剪接形式，几乎在所有癌细胞系中都有表达，但在正常细胞中表达水平非常低。””(引文)。在这里，它似乎与第8类有关。

让我们概述一下每个类别中最具特征的术语：

tf_idf %>%
  arrange(desc(tf_idf)) %>%
  mutate(word = factor(word, levels = rev(unique(word)))) %>%
  group_by(Class) %>%
  top_n(10, tf_idf) %>%
  ungroup() %>%  
  ggplot(aes(word, tf_idf, fill = Class)) +
  geom_col() +
  labs(x = NULL, y = "tf-idf") +
  theme(legend.position = "none") +
  facet_wrap(~ Class, ncol = 3, scales = "free") +
  coord_flip()

这是非常专业的术语。但是，我们注意到，其中一些（如“brct”）出现在多个类别中，但仍然具有较高的tf-idf。

3.7 词对频率：n元语法

与测量单个单词的频率类似，我们还可以研究同时出现的单词组的属性（如“统计分析”）。这让我们了解了特定文档中单词之间的（典型）关系。

Tidytext和其他工具使用n-gram的概念，其中n是我们想要作为一个组学习的相邻单词的数量。例如，二元ram是由两个单词组成的一对。我们可以像提取单个单词一样，提取所有这些对：

t2 <- train_txt %>% select(ID, txt) %>% unnest_tokens(bigram, txt, token = "ngrams", n = 2)
head(t2)

为了过滤掉停用词，我们需要首先分离二元分词，然后在过滤后将它们重新组合在一起。Separate/unite也是对应的dplyr函数的名称：

bi_sep <- t2 %>%
  separate(bigram, c("word1", "word2"), sep = " ")

bi_filt <- bi_sep %>%
  filter(!word1 %in% stop_words$word) %>%
  filter(!word2 %in% stop_words$word) %>%
  filter(!word1 %in% my_stopwords$word) %>%
  filter(!word2 %in% my_stopwords$word)

# for later
bigram_counts <- bi_filt %>%
  count(word1, word2, sort = TRUE)

t2 <- bi_filt %>%
  unite(bigram, word1, word2, sep = " ")

估计tf-idf:

foo <- train %>%
  select(ID, Class)

t2_class <- full_join(t2, foo, by = "ID")

t2_tf_idf <- t2_class %>%
  count(Class, bigram) %>%
  bind_tf_idf(bigram, Class, n) %>%
  arrange(desc(tf_idf))

并将每个类别的二元图绘制为具有最佳tf-idf值：

t2_tf_idf %>%
  arrange(desc(tf_idf)) %>%
  mutate(bigram = factor(bigram, levels = rev(unique(bigram)))) %>%
  group_by(Class) %>%
  top_n(10, tf_idf) %>%
  ungroup() %>%  
  ggplot(aes(bigram, tf_idf, fill = Class)) +
  geom_col() +
  labs(x = NULL, y = "tf-idf") +
  theme(legend.position = "none") +
  facet_wrap(~ Class, ncol = 3, scales = "free") +
  coord_flip()

注意，这里我们没有将相似的单词归约为相同的词干，这导致了类内相似的出现（例如Class == 8中的“dnmt3b7表达”和“dnmt3b7表达”）。尽管如此，总的来说，类的内容看起来足够不同，对预测有用。

3.8二元分词网络

一旦我们有了二元分词，即相邻单词的序列，我们还可以通过构建网络来可视化它们与其他单词的联系。单词网络是连接节点的组合。这里我们使用igraph包来构建网络，并使用ggraph包在tidyverse的上下文中将其可视化：

bigram_graph <- bigram_counts %>%
  filter(n > 4e3) %>%
  graph_from_data_frame()

set.seed(1234)

a <- grid::arrow(type = "closed", length = unit(.1, "inches"))

ggraph(bigram_graph, layout = "fr") +
  geom_edge_link(aes(edge_alpha = n), show.legend = FALSE,
                 arrow = a, end_cap = circle(.07, 'inches')) +
  geom_node_point(color = "lightblue", size = 3) +
  geom_node_text(aes(label = name), vjust = 1, hjust = 1) +
  theme_void()

也许这些网络对于解决这个特定的问题不是那么重要。但它们可以让像我这样的非生物学家更多地了解各种技术概念是如何联系在一起的。

这里的箭头表示关系一词的方向（例如“基因表达”而不是“表达基因”）。根据这些链接箭头出现的频率，将其应用于透明度（罕见的箭头更透明）。

3.9 个体类别网络

让我们为各个类绘制相同的网络图，以研究它们的特定重要性术语。为了使其工作，我们需要分别提取二元ram计数。为此，我们构建了一个简短的辅助函数，并给它分配了提取在绘图中显示多少二元ram组合的灵活性。这里，函数的第一个参数是类别的数量，第二个参数是二元语法单词组合的下限。

# input parameters: Class name [1:9], minimum count for bigram graph
plot_bigram_net_class <- function(clname, bimin){

  foo <- t2_class %>%
    filter(Class == clname)

  bar <- foo %>%
    separate(bigram, c("word1", "word2"), sep = " ")
  
  bi_filt <- bar %>%
    filter(!word1 %in% stop_words$word) %>%
    filter(!word2 %in% stop_words$word) %>%
    filter(!word1 %in% my_stopwords$word) %>%
    filter(!word2 %in% my_stopwords$word)
  
  bigram_graph <- bi_filt %>%
    count(word1, word2, sort = TRUE) %>%
    filter(n > bimin) %>%
    graph_from_data_frame()
  
  set.seed(1234)

  a <- grid::arrow(type = "closed", length = unit(.1, "inches"))

  ggraph(bigram_graph, layout = "fr") +
    geom_edge_link(aes(edge_alpha = n), show.legend = FALSE,
                 arrow = a, end_cap = circle(.07, 'inches')) +
    geom_node_point(color = "lightblue", size = 3) +
    geom_node_text(aes(label = name), vjust = 1, hjust = 1) +
    theme_void()
}

这是9张不同的图。我们尽量保持网络图相对稀疏，以便我们可以更清楚地看到重要的连接。你可以在这里尝试使用更多的二元分词。

在下文中，我还注意到一些我认为具有特点的术语或组合。作为一个绝对的非专家，我可能也会注意到一些与我们的挑战无关的琐碎术语。随着比赛的进行，我希望获得一些关于如何清理输入数据的提示。

plot_bigram_net_class(1,8e2)

class 1：我们看到了“p53”和“brct”的联系。我们还发现了二元组“肿瘤抑制因子”。

plot_bigram_net_class(2, 6.5e2)

Class 2：我们看到“ba”、“f3”和“3t3”与“癌细胞”的关系。

plot_bigram_net_class(3, 2e2)

Class 3：在这里，“ baf3 ”和“ brca1 ”似乎很重要。也许还有“酪氨酸激酶”。

plot_bigram_net_class(4, 1e3)

Class 4：这里我们又看到了“brca1”和“brct”，还有另一个突出的“肿瘤抑制因子”。

plot_bigram_net_class(5, 5e2)

Class 5：我们有“顺铂敏感性”和“brca1”网络。

plot_bigram_net_class(6, 5e2)

Class 6：再次“抑瘤”和“e2相互作用”。

plot_bigram_net_class(7, 2e3)

Class 7：这里，“egfr”似乎对“突变”很重要，可以发现几个孤立的二元分词。

plot_bigram_net_class(8, 5e1)

Class 8：在这里我们看到3个术语相对较多的联系，如“bor”与“ccnb3”、“rara”或“gbm”与“adult”和“pediatric”。

plot_bigram_net_class(9, 1.2e2)

class 9：图中一个较密集的网络显示了连接“idh1”和“u2af1”的关系。