Подгружаем модули:

library(tidyverse)
## Warning in Sys.timezone(): unable to identify current timezone '!0@0B>2 ':
## please set environment variable 'TZ'
library(readxl)
library(plotly)
library(htmlwidgets)
library(corrplot)
library(ggcorrplot)
library(pROC)

Препроцессинг:

Загрузим таблицу, объединим одинаковые гены из tox/antitox в один столбец, преобразуем все категориальные переменные в факторы:

dat <- read_xlsx("Table_main_without_SYN.xlsx")
colnames(dat)[8] <- "vca_0311 / vca_0385"
dat <- dat[-14]
which(colnames(dat) == "vca_0312")
## [1] 9
colnames(dat)[9] <- "vca_0312 / vca_0386"
dat <- dat[-14]
colnames(dat)[12] <- "vca_0348 / vca_0503"
dat <- dat[-18]
colnames(dat)[13] <- "vca_0349 / vca_0504"
dat <- dat[-18]

dat$cas3 <- replace(dat$cas3, dat$cas3 == "n.f.", "-")
dat$cas3 <- replace(dat$cas3, dat$cas3 == "+**", "+")
dat_fct <- dat %>% 
  mutate(across(vc_1765:phage_sens, as.factor))

Добавим информацию о времени и месте выделения:

dat_year <- read_xlsx("Table1_year_place.xlsx", col_names = FALSE)
## New names:
## • `` -> `...1`
## • `` -> `...2`
colnames(dat_year) <- c("Strain", "Place/year")
dat_year <- dat_year %>% 
  separate(`Place/year`, into = c("Country", "City", "Year", "Source"), sep = ",")
## Warning: Expected 4 pieces. Missing pieces filled with `NA` in 4 rows [13, 14,
## 18, 20].
dat_year$Source[c(13, 14, 18, 20)] <- dat_year$Year[c(13, 14, 18, 20)]
dat_year$Source[28] <- "env"
dat_year$Source <- replace(dat_year$Source, dat_year$Source %in% c(" вн. ср.", " вн.ср."), "env")
dat_year$Source <- replace(dat_year$Source, dat_year$Source == " чел.", "patient")
dat_year$Year[c(13, 14, 18, 20)] <- dat_year$City[c(13, 14, 18, 20)]
dat_year$City[c(13, 14, 18, 20)] <- NA
dat_year <- dat_year %>%
  unite(col = "Place", c("Country", "City"), sep = ",", remove = TRUE)
dat_year$Place <- str_remove(dat_year$Place, pattern = ",NA")
dat_year$Source <- replace(dat_year$Source, dat_year$Source == "env", "Environment")
dat_year$Source <- replace(dat_year$Source, dat_year$Source == "patient", "Patient")
dat_un <- bind_cols(c(dat_year, dat_fct))
## New names:
## • `Strain` -> `Strain...1`
## • `Strain` -> `Strain...5`
dat_un <- dat_un[-1]
colnames(dat_un)[4] <- "Strain"
dat_un <- dat_un %>% 
  relocate(Strain, .before = Place)
dat_un <- dat_un %>% 
  mutate(across(Place:Source, as.factor))

Анализ:

Выведем таблицу частот по каждому показателю:

dat_un[2:22] <- lapply(dat_un[2:22], str_replace_all, "[\r\n]", "")
dat_un[2:22] <- lapply(dat_un[2:22], as.factor)
lapply(dat_un[2:22], fct_count, sort = TRUE, prop = TRUE)
## $Place
## # A tibble: 11 × 3
##    f                      n      p
##    <fct>              <int>  <dbl>
##  1 РФ, Элиста            16 0.533 
##  2 Украина                3 0.1   
##  3 РФ, Астрахань          2 0.0667
##  4 РФ, Сочи               2 0.0667
##  5 РФ, Казань             1 0.0333
##  6 РФ, Ростов-на-Дону     1 0.0333
##  7 РФ, Челябинск          1 0.0333
##  8 Туркменистан           1 0.0333
##  9 Украина, Бердянск      1 0.0333
## 10 Украина, Мариуполь     1 0.0333
## 11 Украина, Ялта          1 0.0333
## 
## $Year
## # A tibble: 17 × 3
##    f           n      p
##    <fct>   <int>  <dbl>
##  1 " 2011"     4 0.133 
##  2 " 2012"     4 0.133 
##  3 " 2013"     3 0.1   
##  4 " 2015"     3 0.1   
##  5 " 1999"     2 0.0667
##  6 " 2000"     2 0.0667
##  7 " 2017"     2 0.0667
##  8 " 1972"     1 0.0333
##  9 " 1981"     1 0.0333
## 10 " 1995"     1 0.0333
## 11 " 1996"     1 0.0333
## 12 " 2004"     1 0.0333
## 13 " 2005"     1 0.0333
## 14 " 2006"     1 0.0333
## 15 " 2009"     1 0.0333
## 16 " 2014"     1 0.0333
## 17 " 2018"     1 0.0333
## 
## $Source
## # A tibble: 2 × 3
##   f               n     p
##   <fct>       <int> <dbl>
## 1 Environment    20 0.667
## 2 Patient        10 0.333
## 
## $vc_1765
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 -        21   0.7
## 2 int       9   0.3
## 
## $vc_1769
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 -        21   0.7
## 2 int       9   0.3
## 
## $dncV
## # A tibble: 3 × 3
##   f                  n      p
##   <fct>          <int>  <dbl>
## 1 -                 26 0.867 
## 2 A1003G (S335G)     3 0.1   
## 3 int                1 0.0333
## 
## $capV
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 -        26 0.867
## 2 int       4 0.133
## 
## $vc_0814
## # A tibble: 6 × 3
##   f                                                      n      p
##   <fct>                                              <int>  <dbl>
## 1 int                                                   12 0.4   
## 2 A65C (E22A);T85C (S29P)                                6 0.2   
## 3 A65C (E22A);T85C (S29P);A101G (N34S);G442A (A148T)     5 0.167 
## 4 A65C (E22A);T85C (S29P);A101G (N34S)                   3 0.1   
## 5 -                                                      2 0.0667
## 6 A65C (E22A);G116T (S39I)                               2 0.0667
## 
## $vc_0815
## # A tibble: 12 × 3
##    f                                                                    n      p
##    <fct>                                                            <int>  <dbl>
##  1 int                                                                 12 0.4   
##  2 A10G (M4V);A196T (I66F);G605A (S202N)                                3 0.1   
##  3 -                                                                    2 0.0667
##  4 A196T (I66F);G605A (S202N);C1115G (A372G);G1156T (V386F);G1314A…     2 0.0667
##  5 A196T (I66F);G605A (S202N);T727C (S243P);C1115A (A372E);A1219T …     2 0.0667
##  6 T104A (V35E);A196T (I66F);G605A (S202N)                              2 0.0667
##  7 T5G (L2W);G605A (S202N)                                              2 0.0667
##  8 A196T (I66F);G605A (S202N);A898G (T300A);A925C (N309H)               1 0.0333
##  9 G605A (S202N);A908G (E303G)                                          1 0.0333
## 10 T104A (V35E);A196T (I66F);G605A (S202N);                             1 0.0333
## 11 T104A (V35E);A196T (I66F);G605A (S202N); del GTGTTG p: 1(del VL…     1 0.0333
## 12 T104A (V35E);A196T (I66F);G605A (S202N);del 45 п.н. p: 876-920 …     1 0.0333
## 
## $`vca_0311 / vca_0385`
## # A tibble: 5 × 3
##   f                                           n      p
##   <fct>                                   <int>  <dbl>
## 1 int                                        18 0.6   
## 2 -                                           6 0.2   
## 3 G106A (A36T);A256G (N86D);T260G (F87C);     3 0.1   
## 4 C310T (L104F)                               2 0.0667
## 5 G106A (A36T);A256G (N86D);T260G (F87C)      1 0.0333
## 
## $`vca_0312 / vca_0386`
## # A tibble: 4 × 3
##   f                n      p
##   <fct>        <int>  <dbl>
## 1 int             18 0.6   
## 2 -                6 0.2   
## 3 A223G (T75A)     4 0.133 
## 4 A199G (S67G)     2 0.0667
## 
## $vca_0323
## # A tibble: 3 × 3
##   f                                                n      p
##   <fct>                                        <int>  <dbl>
## 1 int                                             20 0.667 
## 2 -                                                8 0.267 
## 3 del 24 п.н. P: 274-297(del GSHSDLF p: 92-98)     2 0.0667
## 
## $vca_0324
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 int      22 0.733
## 2 -         8 0.267
## 
## $`vca_0348 / vca_0503`
## # A tibble: 7 × 3
##   f                                         n      p
##   <fct>                                 <int>  <dbl>
## 1 int                                      12 0.4   
## 2 A82G (K28E);C141G (I47M)                  8 0.267 
## 3 A82G (K28E);                              4 0.133 
## 4 A82G (K28E)                               2 0.0667
## 5 A82G (K28E);C97G (L33V);C141G (I47M)      2 0.0667
## 6 -                                         1 0.0333
## 7 A82G (K28E);G106A (A36T);C141G (I47M)     1 0.0333
## 
## $`vca_0349 / vca_0504`
## # A tibble: 7 × 3
##   f                                                        n      p
##   <fct>                                                <int>  <dbl>
## 1 int                                                     24 0.8   
## 2 A298G (I100V);CA320TG (A107V)                            1 0.0333
## 3 C10T (P4S)                                               1 0.0333
## 4 C52A (L18I)                                              1 0.0333
## 5 del 258 п.н. (70%) p: 32-289 (del 85 а.к. p: 12-96)      1 0.0333
## 6 del 337 п.н. (91%) p: 33-369(del 113 а.к. p: 11-123)     1 0.0333
## 7 T50C (V17A);C71A (T24K);A298G (I100V);C320T (A107V)      1 0.0333
## 
## $vca_0422
## # A tibble: 9 × 3
##   f                                                       n      p
##   <fct>                                               <int>  <dbl>
## 1 G10A (V4I)                                             10 0.333 
## 2 -                                                       9 0.3   
## 3 int                                                     3 0.1   
## 4 A244G (I82V)                                            2 0.0667
## 5 A244G (I82V);T287C (I96T)                               2 0.0667
## 6 del 151 п.н. (49%) p: 144-294(del 51 а.к. p: 48-98)     1 0.0333
## 7 del 97 п.н. p: 1-97(frameshift c p:  1, * p: 15)        1 0.0333
## 8 G10A (V4I);                                             1 0.0333
## 9 G10A (V4I);A269T (Q90L)                                 1 0.0333
## 
## $vca_0423
## # A tibble: 12 × 3
##    f                                                                    n      p
##    <fct>                                                            <int>  <dbl>
##  1 int                                                                 13 0.433 
##  2 del 245 п.н. (68%) p: 1-245(frameshift с p: 1, * p: 12)              7 0.233 
##  3 A102T (E34D);AGCCA223GGCCT (SH75GL);G256A (E86K)                     1 0.0333
##  4 A102T (E34D);TTA247GTT (L83V);CAAGAGCAGACGG253AACTGTCGTCAAT(QEQ…     1 0.0333
##  5 AATGAA97CATGAT (NE33HD);AGCCA223GGCCT (SH75GL);A249T (L83F)          1 0.0333
##  6 AATGAA97CATGAT (NE33HD);AGCCA223GGCCT (SH75GL);A249T (L83F);GCA…     1 0.0333
##  7 del 246 п.н. (68%) p 1-246;(del 82 а.к. из 119 p: 1-82)              1 0.0333
##  8 del 251 (69%) п.н. p: 1-251(frameshift с p:1, * p:6)                 1 0.0333
##  9 del 272 п.н. (75%) p: 1-272;(frameshift с p: 1, * p: 3)              1 0.0333
## 10 GCA344TCG (CS115FG)                                                  1 0.0333
## 11 T16C (S6P);C802T (A27V);del 142 п.н. (40%) p:  104-245(frameshi…     1 0.0333
## 12 TTACGTCAAGAGCAGACGGTGAATCCGCCT247GTTCGTAACTGTCGTCAATTGCTCACGGAA…     1 0.0333
## 
## $vca_0444
## # A tibble: 3 × 3
##   f               n      p
##   <fct>       <int>  <dbl>
## 1 int            18 0.6   
## 2 -              10 0.333 
## 3 C67A (Q23K)     2 0.0667
## 
## $vca_0445
## # A tibble: 3 × 3
##   f                n      p
##   <fct>        <int>  <dbl>
## 1 int             18 0.6   
## 2 -               10 0.333 
## 3 G143A (S48N)     2 0.0667
## 
## $cas3
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 -        19 0.633
## 2 +        11 0.367
## 
## $phage_sens
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 +        17 0.567
## 2 -        13 0.433

Почистим от дубликатов:

dat_un$`vca_0311 / vca_0385` <- fct_collapse(dat_un$`vca_0311 / vca_0385`, "G106A (A36T);A256G (N86D);T260G (F87C)" = "G106A (A36T);A256G (N86D);T260G (F87C);")
dat_un$vc_0815 <- fct_collapse(dat_un$vc_0815, "T104A (V35E);A196T (I66F);G605A (S202N)" = "T104A (V35E);A196T (I66F);G605A (S202N);")
dat_un$`vca_0348 / vca_0503` <- fct_collapse(dat_un$`vca_0348 / vca_0503`, "A82G (K28E)" = "A82G (K28E);")
dat_un$vca_0422 <- fct_collapse(dat_un$vca_0422, "G10A (V4I)" = "G10A (V4I);")
lapply(dat_un[2:22], fct_count, sort = TRUE, prop = TRUE)
## $Place
## # A tibble: 11 × 3
##    f                      n      p
##    <fct>              <int>  <dbl>
##  1 РФ, Элиста            16 0.533 
##  2 Украина                3 0.1   
##  3 РФ, Астрахань          2 0.0667
##  4 РФ, Сочи               2 0.0667
##  5 РФ, Казань             1 0.0333
##  6 РФ, Ростов-на-Дону     1 0.0333
##  7 РФ, Челябинск          1 0.0333
##  8 Туркменистан           1 0.0333
##  9 Украина, Бердянск      1 0.0333
## 10 Украина, Мариуполь     1 0.0333
## 11 Украина, Ялта          1 0.0333
## 
## $Year
## # A tibble: 17 × 3
##    f           n      p
##    <fct>   <int>  <dbl>
##  1 " 2011"     4 0.133 
##  2 " 2012"     4 0.133 
##  3 " 2013"     3 0.1   
##  4 " 2015"     3 0.1   
##  5 " 1999"     2 0.0667
##  6 " 2000"     2 0.0667
##  7 " 2017"     2 0.0667
##  8 " 1972"     1 0.0333
##  9 " 1981"     1 0.0333
## 10 " 1995"     1 0.0333
## 11 " 1996"     1 0.0333
## 12 " 2004"     1 0.0333
## 13 " 2005"     1 0.0333
## 14 " 2006"     1 0.0333
## 15 " 2009"     1 0.0333
## 16 " 2014"     1 0.0333
## 17 " 2018"     1 0.0333
## 
## $Source
## # A tibble: 2 × 3
##   f               n     p
##   <fct>       <int> <dbl>
## 1 Environment    20 0.667
## 2 Patient        10 0.333
## 
## $vc_1765
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 -        21   0.7
## 2 int       9   0.3
## 
## $vc_1769
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 -        21   0.7
## 2 int       9   0.3
## 
## $dncV
## # A tibble: 3 × 3
##   f                  n      p
##   <fct>          <int>  <dbl>
## 1 -                 26 0.867 
## 2 A1003G (S335G)     3 0.1   
## 3 int                1 0.0333
## 
## $capV
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 -        26 0.867
## 2 int       4 0.133
## 
## $vc_0814
## # A tibble: 6 × 3
##   f                                                      n      p
##   <fct>                                              <int>  <dbl>
## 1 int                                                   12 0.4   
## 2 A65C (E22A);T85C (S29P)                                6 0.2   
## 3 A65C (E22A);T85C (S29P);A101G (N34S);G442A (A148T)     5 0.167 
## 4 A65C (E22A);T85C (S29P);A101G (N34S)                   3 0.1   
## 5 -                                                      2 0.0667
## 6 A65C (E22A);G116T (S39I)                               2 0.0667
## 
## $vc_0815
## # A tibble: 11 × 3
##    f                                                                    n      p
##    <fct>                                                            <int>  <dbl>
##  1 int                                                                 12 0.4   
##  2 A10G (M4V);A196T (I66F);G605A (S202N)                                3 0.1   
##  3 T104A (V35E);A196T (I66F);G605A (S202N)                              3 0.1   
##  4 -                                                                    2 0.0667
##  5 A196T (I66F);G605A (S202N);C1115G (A372G);G1156T (V386F);G1314A…     2 0.0667
##  6 A196T (I66F);G605A (S202N);T727C (S243P);C1115A (A372E);A1219T …     2 0.0667
##  7 T5G (L2W);G605A (S202N)                                              2 0.0667
##  8 A196T (I66F);G605A (S202N);A898G (T300A);A925C (N309H)               1 0.0333
##  9 G605A (S202N);A908G (E303G)                                          1 0.0333
## 10 T104A (V35E);A196T (I66F);G605A (S202N); del GTGTTG p: 1(del VL…     1 0.0333
## 11 T104A (V35E);A196T (I66F);G605A (S202N);del 45 п.н. p: 876-920 …     1 0.0333
## 
## $`vca_0311 / vca_0385`
## # A tibble: 4 × 3
##   f                                          n      p
##   <fct>                                  <int>  <dbl>
## 1 int                                       18 0.6   
## 2 -                                          6 0.2   
## 3 G106A (A36T);A256G (N86D);T260G (F87C)     4 0.133 
## 4 C310T (L104F)                              2 0.0667
## 
## $`vca_0312 / vca_0386`
## # A tibble: 4 × 3
##   f                n      p
##   <fct>        <int>  <dbl>
## 1 int             18 0.6   
## 2 -                6 0.2   
## 3 A223G (T75A)     4 0.133 
## 4 A199G (S67G)     2 0.0667
## 
## $vca_0323
## # A tibble: 3 × 3
##   f                                                n      p
##   <fct>                                        <int>  <dbl>
## 1 int                                             20 0.667 
## 2 -                                                8 0.267 
## 3 del 24 п.н. P: 274-297(del GSHSDLF p: 92-98)     2 0.0667
## 
## $vca_0324
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 int      22 0.733
## 2 -         8 0.267
## 
## $`vca_0348 / vca_0503`
## # A tibble: 6 × 3
##   f                                         n      p
##   <fct>                                 <int>  <dbl>
## 1 int                                      12 0.4   
## 2 A82G (K28E);C141G (I47M)                  8 0.267 
## 3 A82G (K28E)                               6 0.2   
## 4 A82G (K28E);C97G (L33V);C141G (I47M)      2 0.0667
## 5 -                                         1 0.0333
## 6 A82G (K28E);G106A (A36T);C141G (I47M)     1 0.0333
## 
## $`vca_0349 / vca_0504`
## # A tibble: 7 × 3
##   f                                                        n      p
##   <fct>                                                <int>  <dbl>
## 1 int                                                     24 0.8   
## 2 A298G (I100V);CA320TG (A107V)                            1 0.0333
## 3 C10T (P4S)                                               1 0.0333
## 4 C52A (L18I)                                              1 0.0333
## 5 del 258 п.н. (70%) p: 32-289 (del 85 а.к. p: 12-96)      1 0.0333
## 6 del 337 п.н. (91%) p: 33-369(del 113 а.к. p: 11-123)     1 0.0333
## 7 T50C (V17A);C71A (T24K);A298G (I100V);C320T (A107V)      1 0.0333
## 
## $vca_0422
## # A tibble: 8 × 3
##   f                                                       n      p
##   <fct>                                               <int>  <dbl>
## 1 G10A (V4I)                                             11 0.367 
## 2 -                                                       9 0.3   
## 3 int                                                     3 0.1   
## 4 A244G (I82V)                                            2 0.0667
## 5 A244G (I82V);T287C (I96T)                               2 0.0667
## 6 del 151 п.н. (49%) p: 144-294(del 51 а.к. p: 48-98)     1 0.0333
## 7 del 97 п.н. p: 1-97(frameshift c p:  1, * p: 15)        1 0.0333
## 8 G10A (V4I);A269T (Q90L)                                 1 0.0333
## 
## $vca_0423
## # A tibble: 12 × 3
##    f                                                                    n      p
##    <fct>                                                            <int>  <dbl>
##  1 int                                                                 13 0.433 
##  2 del 245 п.н. (68%) p: 1-245(frameshift с p: 1, * p: 12)              7 0.233 
##  3 A102T (E34D);AGCCA223GGCCT (SH75GL);G256A (E86K)                     1 0.0333
##  4 A102T (E34D);TTA247GTT (L83V);CAAGAGCAGACGG253AACTGTCGTCAAT(QEQ…     1 0.0333
##  5 AATGAA97CATGAT (NE33HD);AGCCA223GGCCT (SH75GL);A249T (L83F)          1 0.0333
##  6 AATGAA97CATGAT (NE33HD);AGCCA223GGCCT (SH75GL);A249T (L83F);GCA…     1 0.0333
##  7 del 246 п.н. (68%) p 1-246;(del 82 а.к. из 119 p: 1-82)              1 0.0333
##  8 del 251 (69%) п.н. p: 1-251(frameshift с p:1, * p:6)                 1 0.0333
##  9 del 272 п.н. (75%) p: 1-272;(frameshift с p: 1, * p: 3)              1 0.0333
## 10 GCA344TCG (CS115FG)                                                  1 0.0333
## 11 T16C (S6P);C802T (A27V);del 142 п.н. (40%) p:  104-245(frameshi…     1 0.0333
## 12 TTACGTCAAGAGCAGACGGTGAATCCGCCT247GTTCGTAACTGTCGTCAATTGCTCACGGAA…     1 0.0333
## 
## $vca_0444
## # A tibble: 3 × 3
##   f               n      p
##   <fct>       <int>  <dbl>
## 1 int            18 0.6   
## 2 -              10 0.333 
## 3 C67A (Q23K)     2 0.0667
## 
## $vca_0445
## # A tibble: 3 × 3
##   f                n      p
##   <fct>        <int>  <dbl>
## 1 int             18 0.6   
## 2 -               10 0.333 
## 3 G143A (S48N)     2 0.0667
## 
## $cas3
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 -        19 0.633
## 2 +        11 0.367
## 
## $phage_sens
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 +        17 0.567
## 2 -        13 0.433

Сохраним оригиналы существующих столбцов:

dat_true <- dat_un[5:22]
colnames(dat_true) <- paste("true", colnames(dat_true), sep = "_")
dat_fin <- bind_cols(c(dat_un, dat_true))

Оставим только те значения, которые встречаются, как минимум у 3 штаммов из выборки, остальные пометим как Other:

dat_fin[-c(1:7, 9, 15, 21, 22, 23:40)] <- lapply(dat_fin[-c(1:7, 9, 15, 21, 22, 23:40)], fct_lump_min, min = 3)
dat_fin$vc_0815 <- fct_lump_min(dat_fin$vc_0815, min = 2)
dat_fin$vc_0815 <- fct_other(dat_fin$vc_0815, drop = c("A196T (I66F);G605A (S202N);T727C (S243P);C1115A (A372E);A1219T (M407L);T1225A (L409I);A1229G (K410R)", "A196T (I66F);G605A (S202N);C1115G (A372G);G1156T (V386F);G1314A (M438I)"))
dat_fin$vc_0814 <- fct_other(dat_fin$vc_0814, drop = c("A65C (E22A);G116T (S39I)", "A65C (E22A);T85C (S29P)", "A65C (E22A);T85C (S29P);A101G (N34S)", "A65C (E22A);T85C (S29P);A101G (N34S);G442A (A148T)"))
lapply(dat_fin[2:22], fct_count, sort = TRUE, prop = TRUE)
## $Place
## # A tibble: 11 × 3
##    f                      n      p
##    <fct>              <int>  <dbl>
##  1 РФ, Элиста            16 0.533 
##  2 Украина                3 0.1   
##  3 РФ, Астрахань          2 0.0667
##  4 РФ, Сочи               2 0.0667
##  5 РФ, Казань             1 0.0333
##  6 РФ, Ростов-на-Дону     1 0.0333
##  7 РФ, Челябинск          1 0.0333
##  8 Туркменистан           1 0.0333
##  9 Украина, Бердянск      1 0.0333
## 10 Украина, Мариуполь     1 0.0333
## 11 Украина, Ялта          1 0.0333
## 
## $Year
## # A tibble: 17 × 3
##    f           n      p
##    <fct>   <int>  <dbl>
##  1 " 2011"     4 0.133 
##  2 " 2012"     4 0.133 
##  3 " 2013"     3 0.1   
##  4 " 2015"     3 0.1   
##  5 " 1999"     2 0.0667
##  6 " 2000"     2 0.0667
##  7 " 2017"     2 0.0667
##  8 " 1972"     1 0.0333
##  9 " 1981"     1 0.0333
## 10 " 1995"     1 0.0333
## 11 " 1996"     1 0.0333
## 12 " 2004"     1 0.0333
## 13 " 2005"     1 0.0333
## 14 " 2006"     1 0.0333
## 15 " 2009"     1 0.0333
## 16 " 2014"     1 0.0333
## 17 " 2018"     1 0.0333
## 
## $Source
## # A tibble: 2 × 3
##   f               n     p
##   <fct>       <int> <dbl>
## 1 Environment    20 0.667
## 2 Patient        10 0.333
## 
## $vc_1765
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 -        21   0.7
## 2 int       9   0.3
## 
## $vc_1769
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 -        21   0.7
## 2 int       9   0.3
## 
## $dncV
## # A tibble: 3 × 3
##   f                  n      p
##   <fct>          <int>  <dbl>
## 1 -                 26 0.867 
## 2 A1003G (S335G)     3 0.1   
## 3 int                1 0.0333
## 
## $capV
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 -        26 0.867
## 2 int       4 0.133
## 
## $vc_0814
## # A tibble: 3 × 3
##   f         n      p
##   <fct> <int>  <dbl>
## 1 Other    16 0.533 
## 2 int      12 0.4   
## 3 -         2 0.0667
## 
## $vc_0815
## # A tibble: 4 × 3
##   f                                           n     p
##   <fct>                                   <int> <dbl>
## 1 int                                        12   0.4
## 2 Other                                      12   0.4
## 3 A10G (M4V);A196T (I66F);G605A (S202N)       3   0.1
## 4 T104A (V35E);A196T (I66F);G605A (S202N)     3   0.1
## 
## $`vca_0311 / vca_0385`
## # A tibble: 4 × 3
##   f                                          n      p
##   <fct>                                  <int>  <dbl>
## 1 int                                       18 0.6   
## 2 -                                          6 0.2   
## 3 G106A (A36T);A256G (N86D);T260G (F87C)     4 0.133 
## 4 Other                                      2 0.0667
## 
## $`vca_0312 / vca_0386`
## # A tibble: 4 × 3
##   f                n      p
##   <fct>        <int>  <dbl>
## 1 int             18 0.6   
## 2 -                6 0.2   
## 3 A223G (T75A)     4 0.133 
## 4 Other            2 0.0667
## 
## $vca_0323
## # A tibble: 3 × 3
##   f         n      p
##   <fct> <int>  <dbl>
## 1 int      20 0.667 
## 2 -         8 0.267 
## 3 Other     2 0.0667
## 
## $vca_0324
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 int      22 0.733
## 2 -         8 0.267
## 
## $`vca_0348 / vca_0503`
## # A tibble: 6 × 3
##   f                                         n      p
##   <fct>                                 <int>  <dbl>
## 1 int                                      12 0.4   
## 2 A82G (K28E);C141G (I47M)                  8 0.267 
## 3 A82G (K28E)                               6 0.2   
## 4 A82G (K28E);C97G (L33V);C141G (I47M)      2 0.0667
## 5 -                                         1 0.0333
## 6 A82G (K28E);G106A (A36T);C141G (I47M)     1 0.0333
## 
## $`vca_0349 / vca_0504`
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 int      24   0.8
## 2 Other     6   0.2
## 
## $vca_0422
## # A tibble: 4 × 3
##   f              n     p
##   <fct>      <int> <dbl>
## 1 G10A (V4I)    11 0.367
## 2 -              9 0.3  
## 3 Other          7 0.233
## 4 int            3 0.1  
## 
## $vca_0423
## # A tibble: 3 × 3
##   f                                                           n     p
##   <fct>                                                   <int> <dbl>
## 1 int                                                        13 0.433
## 2 Other                                                      10 0.333
## 3 del 245 п.н. (68%) p: 1-245(frameshift с p: 1, * p: 12)     7 0.233
## 
## $vca_0444
## # A tibble: 3 × 3
##   f         n      p
##   <fct> <int>  <dbl>
## 1 int      18 0.6   
## 2 -        10 0.333 
## 3 Other     2 0.0667
## 
## $vca_0445
## # A tibble: 3 × 3
##   f         n      p
##   <fct> <int>  <dbl>
## 1 int      18 0.6   
## 2 -        10 0.333 
## 3 Other     2 0.0667
## 
## $cas3
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 -        19 0.633
## 2 +        11 0.367
## 
## $phage_sens
## # A tibble: 2 × 3
##   f         n     p
##   <fct> <int> <dbl>
## 1 +        17 0.567
## 2 -        13 0.433
levels(dat_fin$`vca_0348 / vca_0503`)
## [1] "-"                                    
## [2] "A82G (K28E)"                          
## [3] "A82G (K28E);C141G (I47M)"             
## [4] "A82G (K28E);C97G (L33V);C141G (I47M)" 
## [5] "A82G (K28E);G106A (A36T);C141G (I47M)"
## [6] "int"

Преобразуем датафрейм в длинный формат:

dat_long <- as_tibble(dat_fin[1:22])
dat_long <- pivot_longer(dat_long, cols = 5:22, names_to = "Gene", values_to = "Value")
dat_long_2 <- as_tibble(dat_fin[c(1:4, 23:40)])
dat_long_2 <- dat_long_2 %>% 
  relocate(true_cas3:true_phage_sens, .after = true_vca_0445)
dat_long_2 <- pivot_longer(dat_long_2, cols = c(5:22), names_to = "Gene_true", values_to = "Value_true")

Сделаем так, чтоб наш датафрейм содержал одновременно истинные и схлопнутые значения:

dat_unit <- cbind(dat_long, dat_long_2)
dat_unit <- dat_unit[-c(7:10)]
dat_unit$Gene <- as.factor(dat_unit$Gene)
dat_unit$Gene_true <- as.factor(dat_unit$Gene_true)

Добавим числовой столбец, обозначающий порядок расположения штаммов на графике:

dat_unit <- dat_unit %>% 
  mutate(col_fct = rep(1:30, each = 18))

Создадим палитру:

pal <- c(
  "coral", "beige",
  "darkseagreen2",
  "steelblue",
  "lightpink",
  "darkolivegreen4", "gold1",
  "thistle3", "slategrey",
  "palegreen2",
  "skyblue", 
  "#FDBF6F",
  "grey", "khaki2",
  "maroon",  "beige"
)

Создадим интерактивный график: цвета будут соответствовать схлопнутой таблице, а при наведении будет доступна полная информация:

p <- ggplot(data = dat_unit, aes(x = Gene, y = fct_reorder(Strain, rev(col_fct)), fill = as.factor(Value),
                                  text = paste("Strain: ", `Strain`,
                                               "<br>Year: ", `Year`,
                                               "<br>Gene: ", Gene,
                                               "<br>Str: ", Value_true))) +
  geom_tile(color = "black") +
  labs(title = "Structure of antiphage system's genes") +
  xlab(label = "Strain") +
  ylab(label = "Gene") +
  scale_fill_manual(values = pal, guide = guide_legend(title = NULL)) +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 7), axis.text.x = element_text(angle = 45))

ggplotly(p, tooltip = "text")

Создадим функцию, которая принимает на выход название показателя в виде dat_fin$value и выдает результаты теста chi-squared по всем показателям датафрейма dat_fin, кроме Strain. Функция выводит только те показатели, результаты chisq.test с которыми < 0.05:

get_chisq <- function(column){
  name = substitute(column)
  test <- as_tibble(sapply(dat_fin[2:22], chisq.test, y = dat_fin[[name[[3]]]]), rownames = NA)
  test <- as_tibble(t(test), rownames = NA)
  namecols <- rownames(test)
  test <- test %>% 
    mutate(Col = namecols, .before = statistic)
  test <- test[, c(1, 4)]
  test[2] <- unlist(test[2])
  test <- test %>% 
    filter(p.value < 0.05)
  return(test)
}

Поиск взаимосвязей

Протестируем хи-квадрат всех показателей друг с другом (см. график dat_chisq.xlsx):

dat_chisq <- tibble(colnames(dat_fin)[2:22])
colnames(dat_chisq) <- "Col"

test_Place <- get_chisq(dat_fin$Place)
dat_chisq <- left_join(x = dat_chisq, y = test_Place, by = c("Col"))
colnames(dat_chisq)[2] <- "Place"

test_Year <- get_chisq(dat_fin$Year)
dat_chisq <- left_join(x = dat_chisq, y = test_Year, by = c("Col"))
colnames(dat_chisq)[3] <- "Year"

test_Source <- get_chisq(dat_fin$Source)
dat_chisq <- left_join(x = dat_chisq, y = test_Source, by = c("Col"))
colnames(dat_chisq)[4] <- "Source"

test_vc1765 <- get_chisq(dat_fin$vc_1765)
dat_chisq <- left_join(x = dat_chisq, y = test_vc1765, by = c("Col"))
colnames(dat_chisq)[5] <- "vc_1765"

test_vc1769 <- get_chisq(dat_fin$vc_1769)
dat_chisq <- left_join(x = dat_chisq, y = test_vc1769, by = c("Col"))
colnames(dat_chisq)[6] <- "vc_1769"

test_dncV <- get_chisq(dat_fin$dncV)
dat_chisq <- left_join(x = dat_chisq, y = test_dncV, by = c("Col"))
colnames(dat_chisq)[7] <- "dncV"

test_capV <- get_chisq(dat_fin$capV)
dat_chisq <- left_join(x = dat_chisq, y = test_capV, by = c("Col"))
colnames(dat_chisq)[8] <- "capV"

test_vc_0814 <- get_chisq(dat_fin$vc_0814)
dat_chisq <- left_join(x = dat_chisq, y = test_vc_0814, by = c("Col"))
colnames(dat_chisq)[9] <- "vc_0814"

test_vc_0815 <- get_chisq(dat_fin$vc_0815)
dat_chisq <- left_join(x = dat_chisq, y = test_vc_0815, by = c("Col"))
colnames(dat_chisq)[10] <- "vc_0815"

test_vca_0311_vca_0385 <- get_chisq(dat_fin$`vca_0311 / vca_0385`)
dat_chisq <- left_join(dat_chisq, test_vca_0311_vca_0385, by = c("Col"))
colnames(dat_chisq)[11] <- "vca_0311 / vca_0385"

test_vca_0312_vca_0386 <- get_chisq(dat_fin$`vca_0312 / vca_0386`)
dat_chisq <- left_join(dat_chisq, test_vca_0312_vca_0386, by = c("Col"))
colnames(dat_chisq)[12] <- "vca_0312 / vca_0386"

test_vca_0323 <- get_chisq(dat_fin$vca_0323)
dat_chisq <- left_join(dat_chisq, test_vca_0323, by = c("Col"))
colnames(dat_chisq)[13] <- "vca_0323"         

test_vca_0324 <- get_chisq(dat_fin$vca_0324)
dat_chisq <- left_join(dat_chisq, test_vca_0324, by = c("Col"))
colnames(dat_chisq)[14] <- "vca_0324"

test_vca_0348_vca_0503 <- get_chisq(dat_fin$`vca_0348 / vca_0503`)
dat_chisq <- left_join(dat_chisq, test_vca_0348_vca_0503, by = c("Col"))
colnames(dat_chisq)[15] <- "vca_0348 / vca_0503"

test_vca_0349_vca_0504 <- get_chisq(dat_fin$`vca_0349 / vca_0504`)
dat_chisq <- left_join(dat_chisq, test_vca_0349_vca_0504, by = c("Col"))
colnames(dat_chisq)[16] <- "vca_0349 / vca_0504"

test_vca_0422 <- get_chisq(dat_fin$vca_0422)
dat_chisq <- left_join(dat_chisq, test_vca_0422, by = c("Col"))
colnames(dat_chisq)[17] <- "vca_0422"

test_vca_0423 <- get_chisq(dat_fin$vca_0423)
dat_chisq <- left_join(dat_chisq, test_vca_0423, by = c("Col"))
colnames(dat_chisq)[18] <- "vca_0423"

test_vca_0444 <- get_chisq(dat_fin$vca_0444)
dat_chisq <- left_join(dat_chisq, test_vca_0444, by = c("Col"))
colnames(dat_chisq)[19] <- "vca_0444"

test_vca_0445 <- get_chisq(dat_fin$vca_0445)
dat_chisq <- left_join(dat_chisq, test_vca_0445, by = c("Col"))
colnames(dat_chisq)[20] <- "vca_0445"

test_cas3 <- get_chisq(dat_fin$cas3)
dat_chisq <- left_join(dat_chisq, test_cas3, by = c("Col"))
colnames(dat_chisq)[21] <- "cas3"

test_phage_sens <- get_chisq(dat_fin$phage_sens)
dat_chisq <- left_join(dat_chisq, test_phage_sens, by = c("Col"))
colnames(dat_chisq)[22] <- "phage_sens"

dat_chisq2 <- dat_chisq
dat_chisq2$Place[!is.na(dat_chisq2$Place)] <- "+"
dat_chisq2$Year[!is.na(dat_chisq2$Year)] <- "+"
dat_chisq2$Source[!is.na(dat_chisq2$Source)] <- "+"
dat_chisq2$vc_1765[!is.na(dat_chisq2$vc_1765)] <- "+"
dat_chisq2$vc_1769[!is.na(dat_chisq2$vc_1769)] <- "+"
dat_chisq2$dncV[!is.na(dat_chisq2$dncV)] <- "+"
dat_chisq2$capV[!is.na(dat_chisq2$capV)] <- "+"
dat_chisq2$vc_0814[!is.na(dat_chisq2$vc_0814)] <- "+"
dat_chisq2$vc_0815[!is.na(dat_chisq2$vc_0815)] <- "+"
dat_chisq2$`vca_0311 / vca_0385`[!is.na(dat_chisq2$`vca_0311 / vca_0385`)] <- "+"
dat_chisq2$`vca_0312 / vca_0386`[!is.na(dat_chisq2$`vca_0312 / vca_0386`)] <- "+"
dat_chisq2$vca_0323[!is.na(dat_chisq2$vca_0323)] <- "+"
dat_chisq2$vca_0324[!is.na(dat_chisq2$vca_0324)] <- "+"
dat_chisq2$`vca_0348 / vca_0503`[!is.na(dat_chisq2$`vca_0348 / vca_0503`)] <- "+"
dat_chisq2$`vca_0349 / vca_0504`[!is.na(dat_chisq2$`vca_0349 / vca_0504`)] <- "+"
dat_chisq2$vca_0422[!is.na(dat_chisq2$vca_0422)] <- "+"
dat_chisq2$vca_0423[!is.na(dat_chisq2$vca_0423)] <- "+"
dat_chisq2$vca_0444[!is.na(dat_chisq2$vca_0444)] <- "+"
dat_chisq2$vca_0445[!is.na(dat_chisq2$vca_0445)] <- "+"
dat_chisq2$cas3[!is.na(dat_chisq2$cas3)] <- "+"
dat_chisq2$phage_sens[!is.na(dat_chisq2$phage_sens)] <- "+"
dat_chisq2
## # A tibble: 21 × 22
##    Col            Place Year  Source vc_1765 vc_1769 dncV  capV  vc_0814 vc_0815
##    <chr>          <chr> <chr> <chr>  <chr>   <chr>   <chr> <chr> <chr>   <chr>  
##  1 Place          +     +     +      <NA>    <NA>    +     <NA>  <NA>    <NA>   
##  2 Year           +     +     <NA>   <NA>    <NA>    +     +     <NA>    <NA>   
##  3 Source         +     <NA>  +      <NA>    <NA>    <NA>  <NA>  <NA>    <NA>   
##  4 vc_1765        <NA>  <NA>  <NA>   +       +       <NA>  <NA>  +       +      
##  5 vc_1769        <NA>  <NA>  <NA>   +       +       <NA>  <NA>  +       +      
##  6 dncV           +     +     <NA>   <NA>    <NA>    +     +     <NA>    +      
##  7 capV           <NA>  +     <NA>   <NA>    <NA>    +     +     <NA>    +      
##  8 vc_0814        <NA>  <NA>  <NA>   +       +       <NA>  <NA>  +       +      
##  9 vc_0815        <NA>  <NA>  <NA>   +       +       +     +     +       +      
## 10 vca_0311 / vc… <NA>  <NA>  <NA>   +       +       +     +     +       +      
## # ℹ 11 more rows
## # ℹ 12 more variables: `vca_0311 / vca_0385` <chr>,
## #   `vca_0312 / vca_0386` <chr>, vca_0323 <chr>, vca_0324 <chr>,
## #   `vca_0348 / vca_0503` <chr>, `vca_0349 / vca_0504` <chr>, vca_0422 <chr>,
## #   vca_0423 <chr>, vca_0444 <chr>, vca_0445 <chr>, cas3 <chr>,
## #   phage_sens <chr>

Оставим только показатели, которыу статистически значимо связаны с чувствительностью к фагу:

dat_sens <- dat_fin[c(5, 6, 9, 10, 11, 15, 18:22)]
dat_sens
## # A tibble: 30 × 11
##    vc_1765 vc_1769 vc_0814 vc_0815 `vca_0311 / vca_0385` `vca_0348 / vca_0503`
##    <fct>   <fct>   <fct>   <fct>   <fct>                 <fct>                
##  1 -       -       int     int     int                   int                  
##  2 int     int     int     int     int                   int                  
##  3 int     int     int     int     int                   int                  
##  4 int     int     int     int     int                   int                  
##  5 -       -       int     int     int                   int                  
##  6 -       -       int     int     int                   int                  
##  7 int     int     int     int     int                   int                  
##  8 int     int     int     int     int                   int                  
##  9 int     int     int     int     int                   int                  
## 10 int     int     int     int     int                   int                  
## # ℹ 20 more rows
## # ℹ 5 more variables: vca_0423 <fct>, vca_0444 <fct>, vca_0445 <fct>,
## #   cas3 <fct>, phage_sens <fct>

Рассмотрим, как именно разные варианты структуры гена оказывают влияние на увствительность к фагу (см. график ggcorrplot.png):

model.matrix(~ 0 + ., data = dat_sens) %>% 
  cor(use = "pairwise.complete.obs") %>% 
  ggcorrplot(show.diag = FALSE, type = "lower", lab = FALSE, lab_size = 3, ggtheme = ggplot2::theme_minimal, tl.cex = 3)

Построение моделей:

Построим логистическую регрессионную модель для сокращенной таблицы, чтобы оценить вклад каждого показателя в чувствительность к фагу

Сначала добавим в таблицу вспомогательный столбец is_phage = 0 (устойчивый) и 1 (чувствительный):

dat_sens <- dat_sens %>% 
  mutate(is.phage = ifelse(phage_sens == "-", 0, 1), .after = phage_sens)
dat_sens$is.phage <- as.factor(dat_sens$is.phage)

Строим модель:

logmod <- glm(data = dat_sens, formula = dat_sens$is.phage ~ dat_sens$vc_1765 + dat_sens$vc_1769 + dat_sens$vc_0814 + dat_sens$vc_0815 + dat_sens$`vca_0311 / vca_0385` + dat_sens$`vca_0348 / vca_0503` + dat_sens$vca_0423 + dat_sens$vca_0444 + dat_sens$vca_0445 + dat_sens$cas3, family = "binomial")
## Warning: glm.fit: возникли подогнанные вероятности 0 или 1
summary(logmod)
## 
## Call:
## glm(formula = dat_sens$is.phage ~ dat_sens$vc_1765 + dat_sens$vc_1769 + 
##     dat_sens$vc_0814 + dat_sens$vc_0815 + dat_sens$`vca_0311 / vca_0385` + 
##     dat_sens$`vca_0348 / vca_0503` + dat_sens$vca_0423 + dat_sens$vca_0444 + 
##     dat_sens$vca_0445 + dat_sens$cas3, family = "binomial", data = dat_sens)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.17741  -0.00002   0.00000   0.00002   1.17741  
## 
## Coefficients: (7 not defined because of singularities)
##                                                                        Estimate
## (Intercept)                                                           9.290e+01
## dat_sens$vc_1765int                                                   3.785e+01
## dat_sens$vc_1769int                                                          NA
## dat_sens$vc_0814int                                                  -1.102e+02
## dat_sens$vc_0814Other                                                -7.033e+01
## dat_sens$vc_0815int                                                          NA
## dat_sens$vc_0815T104A (V35E);A196T (I66F);G605A (S202N)              -6.770e+01
## dat_sens$vc_0815Other                                                -9.026e+01
## dat_sens$`vca_0311 / vca_0385`G106A (A36T);A256G (N86D);T260G (F87C) -1.128e+02
## dat_sens$`vca_0311 / vca_0385`int                                    -2.257e+01
## dat_sens$`vca_0311 / vca_0385`Other                                  -6.243e+01
## dat_sens$`vca_0348 / vca_0503`A82G (K28E)                             6.770e+01
## dat_sens$`vca_0348 / vca_0503`A82G (K28E);C141G (I47M)               -7.061e-15
## dat_sens$`vca_0348 / vca_0503`A82G (K28E);C97G (L33V);C141G (I47M)           NA
## dat_sens$`vca_0348 / vca_0503`A82G (K28E);G106A (A36T);C141G (I47M)   7.554e-12
## dat_sens$`vca_0348 / vca_0503`int                                            NA
## dat_sens$vca_0423int                                                  6.243e+01
## dat_sens$vca_0423Other                                                2.257e+01
## dat_sens$vca_0444int                                                 -7.565e-12
## dat_sens$vca_0444Other                                                       NA
## dat_sens$vca_0445int                                                         NA
## dat_sens$vca_0445Other                                                       NA
## dat_sens$cas3+                                                        4.513e+01
##                                                                      Std. Error
## (Intercept)                                                           1.709e+05
## dat_sens$vc_1765int                                                   3.554e+04
## dat_sens$vc_1769int                                                          NA
## dat_sens$vc_0814int                                                   2.096e+05
## dat_sens$vc_0814Other                                                 1.138e+05
## dat_sens$vc_0815int                                                          NA
## dat_sens$vc_0815T104A (V35E);A196T (I66F);G605A (S202N)               7.620e+04
## dat_sens$vc_0815Other                                                 9.017e+04
## dat_sens$`vca_0311 / vca_0385`G106A (A36T);A256G (N86D);T260G (F87C)  1.229e+05
## dat_sens$`vca_0311 / vca_0385`int                                     8.348e+04
## dat_sens$`vca_0311 / vca_0385`Other                                   1.039e+05
## dat_sens$`vca_0348 / vca_0503`A82G (K28E)                             1.022e+05
## dat_sens$`vca_0348 / vca_0503`A82G (K28E);C141G (I47M)                6.816e+04
## dat_sens$`vca_0348 / vca_0503`A82G (K28E);C97G (L33V);C141G (I47M)           NA
## dat_sens$`vca_0348 / vca_0503`A82G (K28E);G106A (A36T);C141G (I47M)   1.181e+05
## dat_sens$`vca_0348 / vca_0503`int                                            NA
## dat_sens$vca_0423int                                                  5.154e+04
## dat_sens$vca_0423Other                                                4.820e+04
## dat_sens$vca_0444int                                                  6.816e+04
## dat_sens$vca_0444Other                                                       NA
## dat_sens$vca_0445int                                                         NA
## dat_sens$vca_0445Other                                                       NA
## dat_sens$cas3+                                                        5.903e+04
##                                                                      z value
## (Intercept)                                                            0.001
## dat_sens$vc_1765int                                                    0.001
## dat_sens$vc_1769int                                                       NA
## dat_sens$vc_0814int                                                   -0.001
## dat_sens$vc_0814Other                                                 -0.001
## dat_sens$vc_0815int                                                       NA
## dat_sens$vc_0815T104A (V35E);A196T (I66F);G605A (S202N)               -0.001
## dat_sens$vc_0815Other                                                 -0.001
## dat_sens$`vca_0311 / vca_0385`G106A (A36T);A256G (N86D);T260G (F87C)  -0.001
## dat_sens$`vca_0311 / vca_0385`int                                      0.000
## dat_sens$`vca_0311 / vca_0385`Other                                   -0.001
## dat_sens$`vca_0348 / vca_0503`A82G (K28E)                              0.001
## dat_sens$`vca_0348 / vca_0503`A82G (K28E);C141G (I47M)                 0.000
## dat_sens$`vca_0348 / vca_0503`A82G (K28E);C97G (L33V);C141G (I47M)        NA
## dat_sens$`vca_0348 / vca_0503`A82G (K28E);G106A (A36T);C141G (I47M)    0.000
## dat_sens$`vca_0348 / vca_0503`int                                         NA
## dat_sens$vca_0423int                                                   0.001
## dat_sens$vca_0423Other                                                 0.000
## dat_sens$vca_0444int                                                   0.000
## dat_sens$vca_0444Other                                                    NA
## dat_sens$vca_0445int                                                      NA
## dat_sens$vca_0445Other                                                    NA
## dat_sens$cas3+                                                         0.001
##                                                                      Pr(>|z|)
## (Intercept)                                                             1.000
## dat_sens$vc_1765int                                                     0.999
## dat_sens$vc_1769int                                                        NA
## dat_sens$vc_0814int                                                     1.000
## dat_sens$vc_0814Other                                                   1.000
## dat_sens$vc_0815int                                                        NA
## dat_sens$vc_0815T104A (V35E);A196T (I66F);G605A (S202N)                 0.999
## dat_sens$vc_0815Other                                                   0.999
## dat_sens$`vca_0311 / vca_0385`G106A (A36T);A256G (N86D);T260G (F87C)    0.999
## dat_sens$`vca_0311 / vca_0385`int                                       1.000
## dat_sens$`vca_0311 / vca_0385`Other                                     1.000
## dat_sens$`vca_0348 / vca_0503`A82G (K28E)                               0.999
## dat_sens$`vca_0348 / vca_0503`A82G (K28E);C141G (I47M)                  1.000
## dat_sens$`vca_0348 / vca_0503`A82G (K28E);C97G (L33V);C141G (I47M)         NA
## dat_sens$`vca_0348 / vca_0503`A82G (K28E);G106A (A36T);C141G (I47M)     1.000
## dat_sens$`vca_0348 / vca_0503`int                                          NA
## dat_sens$vca_0423int                                                    0.999
## dat_sens$vca_0423Other                                                  1.000
## dat_sens$vca_0444int                                                    1.000
## dat_sens$vca_0444Other                                                     NA
## dat_sens$vca_0445int                                                       NA
## dat_sens$vca_0445Other                                                     NA
## dat_sens$cas3+                                                          0.999
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 41.0539  on 29  degrees of freedom
## Residual deviance:  5.5452  on 14  degrees of freedom
## AIC: 37.545
## 
## Number of Fisher Scoring iterations: 21

Мы не получили ни одного статистически значимого (на 5%-ном уровне значимости) коэффициента. Это может говорить о том, что ни один из выбранных показателей по одиночке не является ключевым или определяющим чувствительность к фагу. В то же время, мы получили поразительно точную модель. Это можно увидеть, если сравнить реальные данные предсказанными моделью (последние два столбца):

dat_sens <- dat_sens %>% 
  mutate(predicted_values = ifelse(logmod$fitted.values > mean(logmod$fitted.values), 1, 0), .after = is.phage)
dat_sens %>% 
  select(12:13) %>%
  mutate(Strain = dat_fin$Strain, .before = is.phage) %>% 
  print.data.frame()
##     Strain is.phage predicted_values
## 1    M1395        1                1
## 2       56        1                1
## 3      866        1                1
## 4       85        1                1
## 5   P18778        1                1
## 6    M1501        1                1
## 7    M1504        1                1
## 8    M1518        1                1
## 9    M1524        1                1
## 10    2613        1                1
## 11    2687        1                1
## 12     124        1                1
## 13    M988        1                1
## 14     617        0                0
## 15   М1332        0                0
## 16   М1337        0                0
## 17 P-18748        0                0
## 18     102        0                0
## 19   M1457        0                0
## 20    2403        0                0
## 21   M1506        0                0
## 22   M1516        0                0
## 23   M1517        1                0
## 24   M1526        0                0
## 25      29        0                0
## 26     132        1                1
## 27   M1522        1                1
## 28     433        1                0
## 29    3178        0                0
## 30     136        0                0

Вычислим эффективность и специфичность полученной модели:

# Верно предсказанные 0 (True Negative):
TN <- dat_sens %>% 
  filter(is.phage == 0, predicted_values == 0) %>% 
  nrow()
# Ложно предсказанные 1 (False Positive):
FP <- dat_sens %>% 
  filter(is.phage == 0, predicted_values == 1) %>% 
  nrow()
# Ложно предсказанные 0 (False Negative):
FN <- dat_sens %>% 
  filter(is.phage == 1, predicted_values == 0) %>% 
  nrow()
# Верно предсказанные 1 (True Positive):
TP <- dat_sens %>% 
  filter(is.phage == 1, predicted_values == 1) %>% 
  nrow()

# Чувствительность:
sens <- TP / (TP + FN)
sens
## [1] 0.8823529
# Специфичность модели:
specif <- TN / (TN + FP)
specif
## [1] 1
# Чувствительность и специфичность полученной модели = 0.82 и 1. Другими словами наша модель нашла 82% из всех истинно положительных (чувствительность к фагу) исходов. Также в 100% случаев, модель правильно предсказала 0 (устойчивость к фагу) среди общего числа истинно отрицательных исходов.

Построим ROC-кривую, еще один инструмент, который поможет наглядно отобразить качество исследуемой нами модели.

proc <- roc(dat_sens$is.phage ~ logmod$fitted.values)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
proc
## 
## Call:
## roc.formula(formula = dat_sens$is.phage ~ logmod$fitted.values)
## 
## Data: logmod$fitted.values in 13 controls (dat_sens$is.phage 0) < 17 cases (dat_sens$is.phage 1).
## Area under the curve: 0.991
plot(proc)
text(x = 1.2, y = 0.9, "AUC = 0.99")

Площадь под ROC-кривой (AUC) модет принимать значения от 0 до 1. Высокие значения AUC, близкие к 1, соотвествуют высокому качеству регрессионной модели. В нашем случае AUC = 0.99, в связи с чем можно говорить об очень высокой предсказательной силе модели.

Выводы:

colnames(dat_sens[1:10])
##  [1] "vc_1765"             "vc_1769"             "vc_0814"            
##  [4] "vc_0815"             "vca_0311 / vca_0385" "vca_0348 / vca_0503"
##  [7] "vca_0423"            "vca_0444"            "vca_0445"           
## [10] "cas3"