2024年，搞定tidyverse，让数据更出彩！

纯生信操作适合入门，高阶生信分析和数据挖掘仍需要R 或 Python。为此我们不断完善教程，相继推出了纯生信21天，R语言21天和单细胞21天。数据分析的数学基础在20世纪早期就已确立，但直到计算机的出现才使得实际操作成为可能，并得以推广。数据分析是数学与计算机科学相结合的产物。R 语言则是我们数据分析和可视化最常用的工具之一。

R数据科学，大佬博客：https://craig./

在众多R包中，tidyverse包凭借其强大的数据处理和可视化功能，成为许多数据专业人士的宠儿。tidyverse的核心理念是数据整洁性（tidy data）。通过规整一致的数据结构，我们可以更轻松地处理、分析和可视化数据。这种一致性使得协作变得更加高效，同时也降低了出错的可能性。尽管ggplot2也常常被认为是tidyverse包的一部分，但我更倾向于将其看作独立的R包。因此，我们这次主要梳理其数据清洁的功能！

tidyverse由一系列相互关联的R包组成，其中跟数据清洁最相关的R包是dplyr和tidyr。dplyr提供一组动词来解决常见的数据清洁操作。①mutate()添加现有函数的新变量；②select()根据名称选择变量；③filter()根据变量的值选择案例；④summarise()将多个值减少为单个摘要；⑤arrange()更改行的顺序。

这些都自然地结合在一起group_by()，允许“按组”执行任何操作。dplyr 旨在抽象数据的存储方式。这意味着除了使用本地数据框外，还可以使用完全相同的 R 代码处理远程数据库表。而在数据处理，尤其是数据框处理的过程中，dplyr是必学的R包之一。处理上述功能之外，dplyr包还支持管道操作（pipe operation，%>%）。

为了学习的方便，我们以R in Action的内容为蓝本，以GEO数据实例展开。在R in Action的3.11.2中，作者通过实例展示了管道操作的使用。我们先运行实例，然后结合GEO数据实战来分享管道操作的使用。

####----按照R in Action的例子展开介绍_dplyr----####rm(list = ls())manager <- c(1, 2, 3, 4, 5)date <- c('10/24/08', '10/28/08', '10/1/08', '10/12/08', '5/1/09')country <- c('US', 'US', 'UK', 'UK', 'UK')gender <- c('M', 'F', 'F', 'M', 'F')age <- c(32, 45, 25, 39, 99)q1 <- c(5, 3, 3, 3, 2)q2 <- c(4, 5, 5, 3, 2)q3 <- c(5, 2, 5, 4, 1)q4 <- c(5, 5, 5, NA, 2)q5 <- c(5, 5, 2, NA, 1)leadership <- data.frame(manager, date, country, gender, age, q1, q2, q3, q4, q5, stringsAsFactors=FALSE)

通过上述代码，我们得到5行10列的表达矩阵，并赋值给leadership。

接下来，我们通过dplyr包的函数对leadership进行操作。

library(dplyr)leadership <- mutate(leadership,                     total_score = q1   q2   q3   q4   q5,                     mean_score = total_score / 5)

首先加载dplyr包，通过mutate函数增加两列。

high_potentials <- filter(leadership, total_score > 10)

通过filter函数筛选总分 >10的候选者。

high_potentials <- select(high_potentials, manager, country, mean_score)

通过select函数，选择manager、country和mean_score。

high_potentials <- arrange(high_potentials, country, mean_score)

通过arrange函数，按照country和mean_score从小到大的顺序排列。

上述操作，可以通过管道操作，一步完成。这种操作忽略了中间过程，对操作熟练的人来说，可以节省很多时间；对初学者来说，有点难度。

high_potentials <- filter(leadership, total_score > 10) %>%  select(manager, country, mean_score) %>%  arrange(country, mean_score)

教材实例展示到此结束。

我们仍然以数据集GSE66957为例介绍。GSE表达矩阵的加载，数据Log2转换的判断和数据标准化处理同前面（GEO数据实战，合并merge讨论）。差别在探针ID转换的步骤上。

####----003_ID_transform_not_in_platfomMap----####gpl <- getGEO('GPL15048', destdir = '.')colnames(Table(gpl))# ids <- Table(gpl)[,c(1,4)]ids <- Table(gpl)[,c(seq(1:5))]ids1 <- select(ids,ID,GeneSymbol)ids1 <- ids1ids2 <- rename(ids1,probeID ='ID', Symbol='GeneSymbol')

我们可以通过select函数和rename函数对GSE66957数据集的信息进行选取和重命名，但实际上没有必要。在本次实操中，我们主要利用mutate函数和select函数获取清洁数据。

通过前面的函数，我们获得表达矩阵exprSet，探针为行，样品名为列。

探针矩阵ids1，有的探针没有GeneSymbol。

接下来我们利用mutate函数把表达矩阵和探针矩阵合并。

library(dplyr)exprset1 <-mutate(exprset,ids1)exprset1 <-as.data.frame(exprset1)

探针矩阵会加在表达矩阵的最后两类，ID和GeneSymbol均保留，最终得到第1~69列为样品名，第70列为ID，第71列为GeneSymbol的矩阵。

然后，通过select函数，我们获取GeneSymbol在第一列，ID删除的矩阵。

exprset2 <-select(exprset1,GeneSymbol,seq(1:69))

exprset3= avereps(exprset2[,-1],                          ID = exprset2$GeneSymbol)exprset3 <- as.data.frame(exprset3)

然后探针名称重复的取平均值，表达为数据框。

在该数据集中，对照标本12个，肿瘤标本57个，通过select函数选取任意12个肿瘤标本来进行数据分析。

exprset4 <- select(exprset3,seq(1:12),GSM1634937,GSM1634938,GSM1634939,GSM1634940, GSM1634941,GSM1634942,GSM1634943,GSM1634944,GSM1634945,GSM1634946,                   GSM1634950,GSM1634955)exprset4 =exprset4[-1,]                                     

这样通过上述操作，我们也可以得到行为基因名，列为样品名的清洁数据。上述操作无法通过管道操作，一步完成。重在管道操作的理解和使用。

####----GSE66957_pipe_operation----####library(dplyr)exprset5 <- mutate(exprset,ids1) %>%            select(GeneSymbol,seq(1:12),GSM1634937,GSM1634938,GSM1634939,GSM1634940,                       GSM1634941,GSM1634942,GSM1634943,GSM1634944,GSM1634945,GSM1634946,                       GSM1634950,GSM1634955) group3 <- c(rep('normal',12),rep('cancer',13)) group3 <- factor(group3,levels = c('normal','cancer'))exprset6 = avereps(exprset5[,-1],ID = exprset5$GeneSymbol)%>%          as.data.frame()%>%                    mutate(aver=(GSM1634937 GSM1634938 GSM1634939)/3)%>%         arrange(desc(aver))design <- model.matrix(~0   group3)colnames(design) <- levels(group3)designcontrast.matrix <- makeContrasts( 'cancer - normal', levels = design)fit1 <- lmFit(exprset6,design)fit1 <- contrasts.fit( fit1, contrast.matrix ) fit2 <- eBayes(fit1)allDiff=topTable(fit2,adjust='fdr',coef=1,number=Inf) exprset7 <- select(allDiff,logFC,adj.P.Val) %>%                      arrange(desc(logFC))

通过上述操作，我们得到按照logFC从大到小排序的表达数据。

tidyr包则主要用于数据的separate和extract，常用于GEO数据的基因名转换处理。

install.packages('tidyr')library(tidyr)library(dplyr, warn.conflicts = FALSE)## 示例数据vt_census <- tidycensus::get_decennial( geography = 'block', state = 'VT', county = 'Washington', variables = 'P1_001N', year = 2020)#> Getting data from the 2020 decennial Census#> Using the PL 94-171 Redistricting Data summary file#> Note: 2020 decennial Census data use differential privacy, a technique that#> introduces errors into data to preserve respondent confidentiality.#> ℹ Small counts should be interpreted with caution.#> ℹ See https://www.census.gov/library/fact-sheets/2021/protecting-the-confidentiality-of-the-2020-census-redistricting-data.html for additional guidance.#> This message is displayed once per session.vt_census#> # A tibble: 2,150 × 4#> GEOID NAME variable value#> <chr> <chr> <chr> <dbl>#> 1 500239555021014 Block 1014, Block Group 1, Census Tract 9555.… P1_001N 21#> 2 500239555021015 Block 1015, Block Group 1, Census Tract 9555.… P1_001N 19#> 3 500239555021016 Block 1016, Block Group 1, Census Tract 9555.… P1_001N 0#> 4 500239555021017 Block 1017, Block Group 1, Census Tract 9555.… P1_001N 0#> 5 500239555021018 Block 1018, Block Group 1, Census Tract 9555.… P1_001N 43#> 6 500239555021019 Block 1019, Block Group 1, Census Tract 9555.… P1_001N 68#> 7 500239555021020 Block 1020, Block Group 1, Census Tract 9555.… P1_001N 30#> 8 500239555021021 Block 1021, Block Group 1, Census Tract 9555.… P1_001N 0#> 9 500239555021022 Block 1022, Block Group 1, Census Tract 9555.… P1_001N 18#> 10 500239555021023 Block 1023, Block Group 1, Census Tract 9555.… P1_001N 93#> # … with 2,140 more rows

区分GEOID，获取国家，城市，区街信息。

vt_census |>  select(GEOID) |>   separate_wider_position(    GEOID,    widths = c(state = 2, county = 3, tract = 6, block = 4)  )#> # A tibble: 2,150 × 4#>    state county tract  block#>    <chr> <chr>  <chr>  <chr>#>  1 50    023    955502 1014 #>  2 50    023    955502 1015 #>  3 50    023    955502 1016 #>  4 50    023    955502 1017 #>  5 50    023    955502 1018 #>  6 50    023    955502 1019 #>  7 50    023    955502 1020 #>  8 50    023    955502 1021 #>  9 50    023    955502 1022 #> 10 50    023    955502 1023 #> # … with 2,140 more rows

借助正则表达式。

vt_census |> select(NAME) |> separate_wider_regex( NAME, patterns = c( 'Block ', block = '\\d ', ', ', 'Block Group ', block_group = '\\d ', ', ', 'Census Tract ', tract = '\\d .\\d ', ', ', county = '[^,] ', ', ', state = '.*' ) )#> # A tibble: 2,150 × 5#> block block_group tract county state #> <chr> <chr> <chr> <chr> <chr> #> 1 1014 1 9555.02 Washington County Vermont#> 2 1015 1 9555.02 Washington County Vermont#> 3 1016 1 9555.02 Washington County Vermont#> 4 1017 1 9555.02 Washington County Vermont#> 5 1018 1 9555.02 Washington County Vermont#> 6 1019 1 9555.02 Washington County Vermont#> 7 1020 1 9555.02 Washington County Vermont#> 8 1021 1 9555.02 Washington County Vermont#> 9 1022 1 9555.02 Washington County Vermont#> 10 1023 1 9555.02 Washington County Vermont#> # … with 2,140 more rows

这部分稍微有点难度，但对于理解和处理复杂的字符串数据很有帮助，值得实践应用！