Unfortunately the TSV files have varying column names, named after the primer. Fix the naming using the tidyr package’s gather() function.
Additionally, create a corresponding fragments object to visualize the cut fragments. We will reuse the pp_id column which is unique to each fragment. In our case, creating the fragments by summarising with head() and str_c() are slow, about 20 seconds per file, because there are millions of length 2 pp_id groups. Summary operations are much faster when there are fewer groups. A better fast, though more complicated alternative would be to split into 2 tables and then perform the summary operations. But we instead use cache() from Biobase to speed up these operations.
suppressPackageStartupMessages({
library(Biobase) # cache
library(tidyverse) # %>% gather mutate rename group_by summarise
library(stringr) # str_split str_c
})
read_tsv_sites <- function(file) {
## The restriction site column names are based on the filenames,
## so we need to do standardize the column names we can combine
## the file rows together.
col_sites <- str_split(tools::file_path_sans_ext(file),
"[_.]", simplify = TRUE) %>%
as.character() %>% tail(2)
col_sites_orig <- str_c("loc", col_sites, "site", sep = "_")
## Don't read all the columns; just these 4.
col_names <- c("chr_alias", col_sites_orig, "pp_id")
## We must specify "chr_alias" as a character, otherwise it would
## be assumed to be an integer, and we would lose the X and Y
## chromosomes.
columns <- list(col_character(), col_integer(), col_integer(),
col_integer())
names(columns) <- col_names
names(col_sites_orig) <- col_sites
read_tsv(file, col_types = do.call(cols_only, columns)) %>%
rename_(.dots = col_sites_orig) %>%
gather_(key = "site", value = "loc", gather_cols = col_sites)
}
cache(
sites_ <- lapply(files, read_tsv_sites) %>% bind_rows() %>%
rename_("chrom" = "chr_alias") %>% arrange(chrom, loc) %>%
mutate(chrom = str_c("chr", chrom))
)
cache(
fragments_ <- sites_ %>%
group_by(pp_id) %>%
summarise(chrom = head(chrom, 1),
start = min(loc),
end = max(loc),
name = str_c(site, collapse = "_"))
)
sites_
#> # A tibble: 5,634,098 × 4
#> chrom pp_id site loc
#> <chr> <int> <chr> <int>
#> 1 chr1 10872946 Csp6I 55158
#> 2 chr1 10872946 BfaI 55575
#> 3 chr1 3258 NlaIII 62708
#> 4 chr1 3258 DpnII 63095
#> 5 chr1 12832397 Csp6I 565605
#> 6 chr1 12832397 DpnII 567010
#> 7 chr1 12832401 Csp6I 567565
#> 8 chr1 12832401 DpnII 568209
#> 9 chr1 12833144 DpnII 718860
#> 10 chr1 12833144 Csp6I 719929
#> # ... with 5,634,088 more rows
fragments_
#> # A tibble: 2,611,623 × 5
#> pp_id chrom start end name
#> <int> <chr> <int> <int> <chr>
#> 1 3258 chr1 62708 63095 NlaIII_DpnII
#> 2 8869 chr1 749642 750061 DpnII_NlaIII
#> 3 9015 chr1 768660 769142 DpnII_NlaIII
#> 4 9076 chr1 777269 777656 NlaIII_DpnII
#> 5 9323 chr1 797184 797691 DpnII_NlaIII
#> 6 9586 chr1 820923 821239 NlaIII_DpnII
#> 7 9703 chr1 831435 832016 DpnII_NlaIII
#> 8 9774 chr1 839455 840285 NlaIII_DpnII
#> 9 9803 chr1 845800 846453 DpnII_NlaIII
#> 10 9835 chr1 850402 850704 DpnII_NlaIII
#> # ... with 2,611,613 more rows
