Tidyverse

Tidyverse

Tidyverse - tidy data

Tidyverse

• Todas as funções do tidyverse tem como primeiro argumento dados.

• Pipe (%>% ou |>)

  # o resultado deste código
  f1(f2(f3(x)))
  
  # é igual ao resultado deste código (porém mais fácil de ler)
  x %>%
    f3() %>%
    f2() %>%
    f1()

Tidyverse

9. {lubridate}

8. {forcats}

7. {readr}

6. {tibble}

5. {purrr}

4. {stringr}

3. {tidyr}

2. {ggplot2}

1. {dplyr}

Importação (e exportação)

• No geral podemos usar read_delim():

  readr::read_delim(file = "arquivo.ext",
                    delim = "\t")

• Caso seu arquivo possua comentarios (Comuns em arquivos de bioinfo) pode-se usar o argumento comment. Exemplo de um arquivo vcf:

  # Exemplo de um arquivo vcf
  readr::read_delim(file = "arquivo.vcf",
                    delim = "\t",
                    comment = "##")

tidy - arrumar o dado

• {tidyr} (Wickham, Vaughan, e Girlich 2023);

• pivot_wider()/pivot_longer();

• separate()/unite();

• nest()/unnest();

tidy - prática

Vamos organizar um arquivo VCF para o formato tidy:

#Importar o VCF
vcf <- read_delim(file = "dados_minicurso/dados_aula/TNF.vcf",
                  delim = "\t",
                  comment = "##")

vcf_tidy <- vcf %>% 
  rename("CHR" = `#CHROM`) %>% #renomear a columa do cromossomo
  select(-ID,-FILTER,-INFO) %>%  #excluir colunas ID e FILTER
  # rename_with(.fn = \(x) paste0("sample_",x),.cols = `1`:`8`) #renomear as amostras
  pivot_longer(cols = `1`:`8`, #converte o arquivo pro formato longo
               names_to = "sample",
               values_to = "genotype") %>% 
  separate(col = genotype,
           into = c("GT","AD","DP","GQ","PL"),
           sep = ':') %>% 
  select(-FORMAT)

tidy - prática

Quais as qualidades (GQ) e cobertura (DP) médias?

vcf_tidy %>% 
  mutate(GQ = as.numeric(GQ),
         DP = as.numeric(DP)) %>% #converter para numerico
  group_by(sample) %>% #agrupar por amostras
  summarise(GQ_media = mean(GQ), #calcular médias
            DP_media = mean(DP))

tidy - prática

vcf_tidy %>% 
  mutate(GQ = as.numeric(GQ),
         DP = as.numeric(DP)) %>% #converter para numerico
  select(POS,sample,DP,GQ) %>% 
  ggplot(aes(x =POS,y = GQ/DP, color = sample)) +
  geom_line() +
  facet_wrap(~sample) +
  theme_bw()

Lidando com multiplos arquivos

• Comumente precisamos lidar com multiplos arquivos

• {purrr} (Wickham e Henry 2023) é desenhado para isso;

• Familia de funções map();

  # Aplica a todos os elementos da lista a função fn1
  map(lista, fn1) 
  # Aplica a todos os elementos do vetor a função fn1
  map(vetor, fn1)

map - prática

• A pasta dados_aula/multiplos contém varios arquivos semelhantes;

• list.files() para retornar o nome deles (caminho absoluto);

• map() com a função read_*() apropriada para importar os dados;

map - prática

(files <- list.files(path = "dados_minicurso/dados_aula/multiplos/",
                    pattern = ".txt",
                    full.names = T))

[1] "dados_minicurso/dados_aula/multiplos/ADP_14d.txt"
[2] "dados_minicurso/dados_aula/multiplos/ADP_7d.txt" 
[3] "dados_minicurso/dados_aula/multiplos/CD_14d.txt" 
[4] "dados_minicurso/dados_aula/multiplos/CD_7d.txt"  
[5] "dados_minicurso/dados_aula/multiplos/OB_14d.txt" 
[6] "dados_minicurso/dados_aula/multiplos/OB_7d.txt"  

dados_mult <- files %>% map(read_tsv)

map - prática

• Vamos nomear a lista!

• {stringr} é um pacote cheio de funções para manipular texto;

  (files_names <-
    str_remove_all(string = files,
                 pattern = "dados_minicurso/dados_aula/multiplos/|.txt"))

  [1] "ADP_14d" "ADP_7d"  "CD_14d"  "CD_7d"   "OB_14d"  "OB_7d"  

  names(dados_mult) <- files_names

O pipeline

expressao <- 
  readxl::read_xlsx(path = "dados_minicurso/dados_aula/pipeline/expressao.xlsx")

glimpse(expressao)

Rows: 6
Columns: 35
$ Group     <chr> "Control1", "Control2", "Control3", "OM1", "OM2", "OM3"
$ `HOXA 1`  <dbl> 0.64, 1.89, 1.10, 5.23, 2.77, 3.04
$ `HOXA 2`  <dbl> 1.950, 1.030, 1.103, 54.480, 50.870, 50.170
$ `HOXA 3`  <dbl> 1.609, 0.707, 1.038, 20.164, 12.421, 19.348
$ `HOXA 4`  <dbl> 1.047301, 0.821207, 1.208541, 131.790600, 102.973600, 101.32~
$ `HOXA 5`  <dbl> 1.06, 0.71, 1.55, 7.76, 14.30, 5.78
$ `HOXA 6`  <dbl> 0.99, 0.77, 1.43, 10.63, 13.80, 7.43
$ `HOXA 7`  <dbl> 1.92, 0.49, 2.39, 1.20, 1.90, 0.86
$ `HOXA 9`  <dbl> 0.97, 1.03, 1.00, 6.59, 9.74, 8.01
$ `HOXA 10` <dbl> 1.11, 0.91, 1.00, 0.49, 0.47, 0.61
$ `HOXA 11` <dbl> 1.26, 1.42, 0.66, 4.45, 1.94, 5.92
$ `HOXA 13` <dbl> 1.81, 0.68, 1.03, 0.80, 1.71, 1.12
$ HOTTIP    <dbl> 1.1380, 0.8880, 1.2005, 114.9360, 189.7810, 234.4220
$ HOTAIR    <dbl> 0.78, 0.85, 1.86, 3.92, 2.40, 1.34
$ ALP       <dbl> 0.68, 1.66, 1.07, 9.18, 5.99, 7.42
$ BSP       <dbl> 0.84, 1.69, 0.95, 0.05, 0.03, 0.04
$ OTX       <dbl> 1.427, 0.891, 1.173, 2.528, 3.123, 1.300
$ RUNX2     <dbl> 0.9423, 1.4230, 1.1690, 2.4660, 5.0030, 3.1090
$ OTC       <dbl> 1.62, 0.95, 1.10, 0.02, 0.05, 0.03
$ OTN       <dbl> 1.01, 0.99, 1.00, 0.08, 0.11, 0.09
$ OPN       <dbl> 0.86, 1.07, 1.11, 1.03, 0.83, 0.87
$ COLA1     <dbl> 0.68, 1.66, 1.07, 9.18, 5.99, 7.42
$ OPG       <dbl> 1.02, 1.24, 0.82, 1.09, 0.79, 1.03
$ BMP2      <dbl> 0.91, 1.11, 1.00, 0.15, 0.34, 0.22
$ BMP7      <dbl> 1.479211, 0.733211, 1.041559, 0.026690, 0.074408, 0.044574
$ BMP9      <dbl> 1.53, 0.72, 1.05, 0.13, 0.24, 0.17
$ BMPRII    <dbl> 1.140000, 0.890000, 1.010000, 0.696667, 1.140000, 0.620000
$ BMPRI     <dbl> 0.93, 1.08, 1.00, 0.63, 0.76, 0.69
$ TGFB      <dbl> 0.77, 1.38, 1.03, 0.63, 0.39, 0.50
$ SMAD2     <dbl> 1.01, 0.62, 0.54, 0.85, 1.20, 0.58
$ SMAD3     <dbl> 1.00, 0.90, 1.10, 0.34, 0.50, 0.41
$ SMAD4     <dbl> 0.75, 1.42, 1.03, 0.44, 0.29, 0.36
$ SMAD5     <dbl> 0.65, 1.81, 1.09, 1.13, 0.66, 0.86
$ SMAD6     <dbl> 0.78, 1.35, 1.03, 100.49, 92.88, 103.00
$ SMAD7     <dbl> 0.70, 1.59, 1.06, 0.34, 0.29, 0.30

O pipeline

Organize os dados (tidy data):

(expressao_tidy <-
  expressao %>% 
  pivot_longer(cols = -Group,
               names_to = "Gene",
               values_to = "Value") %>% 
  rename(Group_name = Group) %>% 
  mutate(Gene = str_remove(Gene," "),
         Group = case_when(
           str_detect(Group_name,"Control") ~ "Control",
           str_detect(Group_name,"OM") ~ "OM"
           ),
         Rep = str_remove_all(Group_name,"Control|OM")
         )
  )

# A tibble: 204 x 5
   Group_name Gene   Value Group   Rep  
   <chr>      <chr>  <dbl> <chr>   <chr>
 1 Control1   HOXA1   0.64 Control 1    
 2 Control1   HOXA2   1.95 Control 1    
 3 Control1   HOXA3   1.61 Control 1    
 4 Control1   HOXA4   1.05 Control 1    
 5 Control1   HOXA5   1.06 Control 1    
 6 Control1   HOXA6   0.99 Control 1    
 7 Control1   HOXA7   1.92 Control 1    
 8 Control1   HOXA9   0.97 Control 1    
 9 Control1   HOXA10  1.11 Control 1    
10 Control1   HOXA11  1.26 Control 1    
# i 194 more rows

O pipeline - nested data

expressao_nested <-
  expressao_tidy %>% 
  arrange(Gene) %>% 
  select(-Group_name) %>% #opcional
  group_by(Gene) %>% 
  nest()

O pipeline - aplicar

Imagine que você queira calcular a média da expressão para cada gene

expressao_nested %>% 
  mutate(media = map(data,~mean(.x$Value)))

# A tibble: 34 x 3
# Groups:   Gene [34]
   Gene   data             media    
   <chr>  <list>           <list>   
 1 ALP    <tibble [6 x 3]> <dbl [1]>
 2 BMP2   <tibble [6 x 3]> <dbl [1]>
 3 BMP7   <tibble [6 x 3]> <dbl [1]>
 4 BMP9   <tibble [6 x 3]> <dbl [1]>
 5 BMPRI  <tibble [6 x 3]> <dbl [1]>
 6 BMPRII <tibble [6 x 3]> <dbl [1]>
 7 BSP    <tibble [6 x 3]> <dbl [1]>
 8 COLA1  <tibble [6 x 3]> <dbl [1]>
 9 HOTAIR <tibble [6 x 3]> <dbl [1]>
10 HOTTIP <tibble [6 x 3]> <dbl [1]>
# i 24 more rows

O pipeline - aplicar

Imagine que você queira calcular a média da expressão para cada gene

expressao_nested %>% 
  mutate(media = map(data,~mean(.x$Value))) %>% 
  unnest(media)

# A tibble: 34 x 3
# Groups:   Gene [34]
   Gene   data              media
   <chr>  <list>            <dbl>
 1 ALP    <tibble [6 x 3]>  4.33 
 2 BMP2   <tibble [6 x 3]>  0.622
 3 BMP7   <tibble [6 x 3]>  0.567
 4 BMP9   <tibble [6 x 3]>  0.64 
 5 BMPRI  <tibble [6 x 3]>  0.848
 6 BMPRII <tibble [6 x 3]>  0.916
 7 BSP    <tibble [6 x 3]>  0.6  
 8 COLA1  <tibble [6 x 3]>  4.33 
 9 HOTAIR <tibble [6 x 3]>  1.86 
10 HOTTIP <tibble [6 x 3]> 90.4  
# i 24 more rows

O pipeline - estatística

Agora como calculamos o teste t para os genes?

(HOXA1 <- expressao_tidy %>% 
  filter(Gene == "HOXA1") %>% 
  select(-Group_name))

# A tibble: 6 x 4
  Gene  Value Group   Rep  
  <chr> <dbl> <chr>   <chr>
1 HOXA1  0.64 Control 1    
2 HOXA1  1.89 Control 2    
3 HOXA1  1.1  Control 3    
4 HOXA1  5.23 OM      1    
5 HOXA1  2.77 OM      2    
6 HOXA1  3.04 OM      3    

O pipeline - estatística

• O pacote {broom} (Robinson, Hayes, e Couch 2023) é feito pra lidar com esses dados;

• Função tidy();

  # A tibble: 1 x 10
    estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high
       <dbl>     <dbl>     <dbl>     <dbl>   <dbl>     <dbl>    <dbl>     <dbl>
  1    -2.47      1.21      3.68     -2.87  0.0683      2.84    -5.30     0.358
  # i 2 more variables: method <chr>, alternative <chr>

O pipeline - estatística

• Outra opção “tidy” é o pacote {infer} (Couch et al. 2021);

• A função t_test() simpifica esse processo;

  HOXA1 %>% 
    infer::t_test(formula = Value ~ Group)

  # A tibble: 1 x 7
    statistic  t_df p_value alternative estimate lower_ci upper_ci
        <dbl> <dbl>   <dbl> <chr>          <dbl>    <dbl>    <dbl>
  1     -2.87  2.84  0.0683 two.sided      -2.47    -5.30    0.358

O pipeline - estatística

• Vamos aplicar a todos?

  expressao_nested %>% 
    mutate(teste_t = map(data,
                       ~infer::t_test(.x, Value ~ Group)))

  # A tibble: 34 x 3
  # Groups:   Gene [34]
     Gene   data             teste_t         
     <chr>  <list>           <list>          
   1 ALP    <tibble [6 x 3]> <tibble [1 x 7]>
   2 BMP2   <tibble [6 x 3]> <tibble [1 x 7]>
   3 BMP7   <tibble [6 x 3]> <tibble [1 x 7]>
   4 BMP9   <tibble [6 x 3]> <tibble [1 x 7]>
   5 BMPRI  <tibble [6 x 3]> <tibble [1 x 7]>
   6 BMPRII <tibble [6 x 3]> <tibble [1 x 7]>
   7 BSP    <tibble [6 x 3]> <tibble [1 x 7]>
   8 COLA1  <tibble [6 x 3]> <tibble [1 x 7]>
   9 HOTAIR <tibble [6 x 3]> <tibble [1 x 7]>
  10 HOTTIP <tibble [6 x 3]> <tibble [1 x 7]>
  # i 24 more rows

O pipeline - estatística

• Vamos aplicar a todos?

  expressao_nested %>% 
    mutate(teste_t = map(data,
                       ~infer::t_test(.x, Value ~ Group))) %>% 
    select(-data) %>% 
    unnest(teste_t)

  # A tibble: 34 x 8
  # Groups:   Gene [34]
     Gene   statistic  t_df  p_value alternative estimate  lower_ci upper_ci
     <chr>      <dbl> <dbl>    <dbl> <chr>          <dbl>     <dbl>    <dbl>
   1 ALP        -6.62  2.38 0.0139   two.sided     -6.39    -9.97     -2.81 
   2 BMP2        9.61  3.99 0.000661 two.sided      0.77     0.547     0.993
   3 BMP7        4.78  2.02 0.0405   two.sided      1.04     0.110     1.96 
   4 BMP9        3.88  2.07 0.0571   two.sided      0.92    -0.0668    1.91 
   5 BMPRI       5.41  3.92 0.00600  two.sided      0.31     0.150     0.470
   6 BMPRII      1.10  2.76 0.359    two.sided      0.194   -0.399     0.787
   7 BSP         4.20  2.00 0.0523   two.sided      1.12    -0.0276    2.27 
   8 COLA1      -6.62  2.38 0.0139   two.sided     -6.39    -9.97     -2.81 
   9 HOTAIR     -1.68  2.83 0.196    two.sided     -1.39    -4.11      1.33 
  10 HOTTIP     -5.12  2.00 0.0360   two.sided   -179.    -329.      -28.7  
  # i 24 more rows

O pipeline - mas e gráficos?

expressao_tidy %>% 
  group_by(Gene,Group) %>% 
  summarise(media = mean(Value),
            desvpad = sd(Value))

# A tibble: 68 x 4
# Groups:   Gene [34]
   Gene  Group    media desvpad
   <chr> <chr>    <dbl>   <dbl>
 1 ALP   Control 1.14    0.493 
 2 ALP   OM      7.53    1.60  
 3 BMP2  Control 1.01    0.100 
 4 BMP2  OM      0.237   0.0961
 5 BMP7  Control 1.08    0.375 
 6 BMP7  OM      0.0486  0.0241
 7 BMP9  Control 1.1     0.407 
 8 BMP9  OM      0.18    0.0557
 9 BMPRI Control 1.00    0.0751
10 BMPRI OM      0.693   0.0651
# i 58 more rows

O pipeline - mas e gráficos?

expressao_tidy %>% 
  group_by(Gene,Group) %>% 
  summarise(media = mean(Value),
            desvpad = sd(Value)) %>%
  nest()

# A tibble: 34 x 2
# Groups:   Gene [34]
   Gene   data            
   <chr>  <list>          
 1 ALP    <tibble [2 x 3]>
 2 BMP2   <tibble [2 x 3]>
 3 BMP7   <tibble [2 x 3]>
 4 BMP9   <tibble [2 x 3]>
 5 BMPRI  <tibble [2 x 3]>
 6 BMPRII <tibble [2 x 3]>
 7 BSP    <tibble [2 x 3]>
 8 COLA1  <tibble [2 x 3]>
 9 HOTAIR <tibble [2 x 3]>
10 HOTTIP <tibble [2 x 3]>
# i 24 more rows

O pipeline - mas e gráficos?

(expressao_plots <-
    expressao_tidy %>% 
  group_by(Gene,Group) %>% 
  summarise(media = mean(Value),
            desvpad = sd(Value)) %>%
  nest() %>% 
  mutate(plots = map(data,
                     \(x) x %>%
                       ggplot(aes(x = Group,
                                  y = media,
                                  fill = Group)) +
                       geom_col(alpha = .8) +
                       geom_errorbar(aes(ymin = media-desvpad,
                                          ymax = media+desvpad,
                                          col = Group),
                                     width = 0.2)
                     )
         )
  )

# A tibble: 34 x 3
# Groups:   Gene [34]
   Gene   data             plots 
   <chr>  <list>           <list>
 1 ALP    <tibble [2 x 3]> <gg>  
 2 BMP2   <tibble [2 x 3]> <gg>  
 3 BMP7   <tibble [2 x 3]> <gg>  
 4 BMP9   <tibble [2 x 3]> <gg>  
 5 BMPRI  <tibble [2 x 3]> <gg>  
 6 BMPRII <tibble [2 x 3]> <gg>  
 7 BSP    <tibble [2 x 3]> <gg>  
 8 COLA1  <tibble [2 x 3]> <gg>  
 9 HOTAIR <tibble [2 x 3]> <gg>  
10 HOTTIP <tibble [2 x 3]> <gg>  
# i 24 more rows

O pipeline - mas e gráficos?

expressao_plots$plots[[1]]

O pipeline - mas e gráficos?

(plots_names <- paste0(expressao_plots$Gene,".png"))

 [1] "ALP.png"    "BMP2.png"   "BMP7.png"   "BMP9.png"   "BMPRI.png" 
 [6] "BMPRII.png" "BSP.png"    "COLA1.png"  "HOTAIR.png" "HOTTIP.png"
[11] "HOXA1.png"  "HOXA10.png" "HOXA11.png" "HOXA13.png" "HOXA2.png" 
[16] "HOXA3.png"  "HOXA4.png"  "HOXA5.png"  "HOXA6.png"  "HOXA7.png" 
[21] "HOXA9.png"  "OPG.png"    "OPN.png"    "OTC.png"    "OTN.png"   
[26] "OTX.png"    "RUNX2.png"  "SMAD2.png"  "SMAD3.png"  "SMAD4.png" 
[31] "SMAD5.png"  "SMAD6.png"  "SMAD7.png"  "TGFB.png"  

walk2(.x = expressao_plots$plots,
      .y = plots_names,
      .f = ~ggsave(filename = .y,
                   plot = .x),
      .progress = T)

O pipeline

graph LR;
  df(tibble) --> f1{group_by};
  f1 --> f2{nest};
  f2 --> f3{mutate};
  f3 --> f4{map};

Aplicações

• O foco será APENAS no código;

• Todos os dados e cenários apresentados foram sintetizados para este curso;

• Os processos estatísticos podem não ser os mais indicados ou insuficientes;

• Nenhum estudante do PPG de Genética foi maltratado durante a construção destes dados;

Cenário 1

Cenário 2

Cenário 3

O que mais posso aprimorar?

• Linux/WSL;

• git/github;

• Quarto markdown;

Desafios para casa