Last updated: 2024-08-02
Checks: 7 0
Knit directory: muse/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20200712) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version b0eadf1. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Ignored files:
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: r_packages_4.3.3/
Ignored: r_packages_4.4.0/
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown (analysis/split_column.Rmd) and HTML
(docs/split_column.html) files. If you’ve configured a
remote Git repository (see ?wflow_git_remote), click on the
hyperlinks in the table below to view the files as they were in that
past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | b0eadf1 | Dave Tang | 2024-08-02 | Update VCF example |
| html | 5fa80f2 | Dave Tang | 2022-09-01 | Build site. |
| Rmd | c8588b7 | Dave Tang | 2022-09-01 | Split GTF group column |
| html | d6644d2 | Dave Tang | 2021-06-22 | Build site. |
| Rmd | 8ae8e21 | Dave Tang | 2021-06-22 | Post on splitting a single column of annotations |
Two widely used file formats in bioinformatics, VCF and GTF, have
single columns that are packed with annotation information. This makes
them a bit inconvenient to work with in R when using data frames because
the values need to be unpacked, i.e. split. In addition, this violates
one of the conditions for tidy data, which is that every cell is a
single value. In this post, we will use tools from the
tidyverse to split the values into multiple columns to make
the data easier to work with in R.
To get started, install the tidyverse if you haven’t
already.
if(!require("tidyverse")){
install.packages("tidyverse")
}
library(tidyverse)We will load a small portion of a VCF file using
read_tsv; in addition we will rename #CHROM to
CHROM and then change all the column names to lower
case.
vcf_url <- "https://raw.githubusercontent.com/davetang/learning_vcf_file/master/eg/Pfeiffer.vcf"
read_tsv(vcf_url, comment = "##", show_col_types = FALSE, n_max = 1000) |>
dplyr::rename(CHROM = `#CHROM`) |>
dplyr::rename_with(tolower) -> my_vcf
head(my_vcf)# A tibble: 6 × 10
chrom pos id ref alt qual filter info format manuel
<dbl> <dbl> <chr> <chr> <chr> <dbl> <chr> <chr> <chr> <chr>
1 1 866511 rs60722469 C CCCCT 259. PASS AC=2;AF=1.00;A… GT:AD… 1/1:6…
2 1 879317 rs7523549 C T 151. PASS AC=1;AF=0.50;A… GT:AD… 0/1:1…
3 1 879482 . G C 485. PASS AC=1;AF=0.50;A… GT:AD… 0/1:2…
4 1 880390 rs3748593 C A 288. PASS AC=1;AF=0.50;A… GT:AD… 0/1:1…
5 1 881627 rs2272757 G A 486. PASS AC=1;AF=0.50;A… GT:AD… 0/1:1…
6 1 884091 rs7522415 C G 65.5 PASS AC=1;AF=0.50;A… GT:AD… 0/1:6…Note that the info column is packed with all sorts of
information for each variant. Also note the consistent format of the
info column: each annotation is separated by a semi-colon
(;) and annotations are stored as key-value pairs with an
equal sign in between.
my_vcf |>
select(info) |>
head()# A tibble: 6 × 1
info
<chr>
1 AC=2;AF=1.00;AN=2;DB;DP=11;FS=0.000;HRun=0;HaplotypeScore=41.3338;MQ0=0;MQ=61…
2 AC=1;AF=0.50;AN=2;BaseQRankSum=1.455;DB;DP=21;Dels=0.00;FS=1.984;HRun=0;Haplo…
3 AC=1;AF=0.50;AN=2;BaseQRankSum=1.934;DP=48;Dels=0.00;FS=4.452;HRun=0;Haplotyp…
4 AC=1;AF=0.50;AN=2;BaseQRankSum=-4.517;DB;DP=29;Dels=0.00;FS=1.485;HRun=0;Hapl…
5 AC=1;AF=0.50;AN=2;BaseQRankSum=0.199;DB;DP=33;Dels=0.00;FS=0.000;HRun=1;Haplo…
6 AC=1;AF=0.50;AN=2;BaseQRankSum=-0.259;DB;DP=12;Dels=0.00;FS=0.000;HRun=1;Hapl…Firstly, we will use separate_rows to create a new row
for each annotation by using ; as the separator/delimiter;
note that I have included \\s* after ;, which
is a regex for specifying a single whitespace occurring 0 or more times.
By including the regex, whitespace after ; will be removed,
which is good because we do not want whitespaces in our data. In
addition, a mutate call is used prior to calling
separate_rows and it is used to remove a trailing
semicolon, if it exists.
my_vcf |>
mutate(info = sub(pattern = ";$", replacement = "", x = .data$info)) |>
separate_rows(info, sep = ";\\s*") |>
head()# A tibble: 6 × 10
chrom pos id ref alt qual filter info format manuel
<dbl> <dbl> <chr> <chr> <chr> <dbl> <chr> <chr> <chr> <chr>
1 1 866511 rs60722469 C CCCCT 259. PASS AC=2 GT:AD:DP:GQ:… 1/1:6…
2 1 866511 rs60722469 C CCCCT 259. PASS AF=1.00 GT:AD:DP:GQ:… 1/1:6…
3 1 866511 rs60722469 C CCCCT 259. PASS AN=2 GT:AD:DP:GQ:… 1/1:6…
4 1 866511 rs60722469 C CCCCT 259. PASS DB GT:AD:DP:GQ:… 1/1:6…
5 1 866511 rs60722469 C CCCCT 259. PASS DP=11 GT:AD:DP:GQ:… 1/1:6…
6 1 866511 rs60722469 C CCCCT 259. PASS FS=0.000 GT:AD:DP:GQ:… 1/1:6…The next step is to split the key-value pairs and we will use the
separate function to separate the pairs into two columns,
which we will name key and value, using the
equal sign as the separator/delimiter. Sometimes a key is missing a
value and in these cases, the value will be NA.
my_vcf |>
mutate(info = sub(pattern = ";$", replacement = "", x = .data$info)) |>
separate_rows(info, sep = ";\\s*") |>
separate(info, c('key', 'value'), sep = "=") |>
head(10)# A tibble: 10 × 11
chrom pos id ref alt qual filter key value format manuel
<dbl> <dbl> <chr> <chr> <chr> <dbl> <chr> <chr> <chr> <chr> <chr>
1 1 866511 rs60722469 C CCCCT 259. PASS AC 2 GT:AD… 1/1:6…
2 1 866511 rs60722469 C CCCCT 259. PASS AF 1.00 GT:AD… 1/1:6…
3 1 866511 rs60722469 C CCCCT 259. PASS AN 2 GT:AD… 1/1:6…
4 1 866511 rs60722469 C CCCCT 259. PASS DB <NA> GT:AD… 1/1:6…
5 1 866511 rs60722469 C CCCCT 259. PASS DP 11 GT:AD… 1/1:6…
6 1 866511 rs60722469 C CCCCT 259. PASS FS 0.000 GT:AD… 1/1:6…
7 1 866511 rs60722469 C CCCCT 259. PASS HRun 0 GT:AD… 1/1:6…
8 1 866511 rs60722469 C CCCCT 259. PASS Haploty… 41.3… GT:AD… 1/1:6…
9 1 866511 rs60722469 C CCCCT 259. PASS MQ0 0 GT:AD… 1/1:6…
10 1 866511 rs60722469 C CCCCT 259. PASS MQ 61.94 GT:AD… 1/1:6…The current state of the transformation produces a new row for each
variant annotation and two columns containing the key and value. If we
want our data in wide format where each annotation is a column, we can
use the pivot_wider function.
In the code below, I have used the first seven columns
(id_cols = chrom:filter) to specify the columns that
uniquely identifies each variant, i.e. the same variant will have the
same values in these columns. We specify our column names from the
key column and the values for these cells will come from
the value column.
my_vcf |>
mutate(info = sub(pattern = ";$", replacement = "", x = .data$info)) |>
separate_rows(info, sep = ";\\s*") |>
separate(info, c('key', 'value'), sep = "=") |>
distinct() |> # remove duplicated annotations, if any
pivot_wider(id_cols = chrom:filter, names_from = key, values_from = value) -> my_vcf_tidy
head(my_vcf_tidy, 10)# A tibble: 10 × 24
chrom pos id ref alt qual filter AC AF AN DB DP
<dbl> <dbl> <chr> <chr> <chr> <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
1 1 866511 rs60722… C CCCCT 259. PASS 2 1.00 2 <NA> 11
2 1 879317 rs75235… C T 151. PASS 1 0.50 2 <NA> 21
3 1 879482 . G C 485. PASS 1 0.50 2 <NA> 48
4 1 880390 rs37485… C A 288. PASS 1 0.50 2 <NA> 29
5 1 881627 rs22727… G A 486. PASS 1 0.50 2 <NA> 33
6 1 884091 rs75224… C G 65.5 PASS 1 0.50 2 <NA> 12
7 1 884101 rs49704… A C 85.8 PASS 1 0.50 2 <NA> 12
8 1 892460 rs41285… G C 1737. PASS 1 0.50 2 <NA> 152
9 1 897730 rs75496… C T 225. PASS 1 0.50 2 <NA> 21
10 1 909238 rs38297… G C 248. PASS 1 0.50 2 <NA> 19
# ℹ 12 more variables: FS <chr>, HRun <chr>, HaplotypeScore <chr>, MQ0 <chr>,
# MQ <chr>, QD <chr>, set <chr>, BaseQRankSum <chr>, Dels <chr>,
# MQRankSum <chr>, ReadPosRankSum <chr>, DS <chr>Now each row is a single variant and each column is a variable,
making it much easier to work with! We can easily subset variants with a
QD > 50 (after transforming the type of numeric).
my_vcf_tidy |>
dplyr::mutate(QD = as.numeric(QD)) |>
dplyr::filter(QD > 50)# A tibble: 6 × 24
chrom pos id ref alt qual filter AC AF AN DB DP
<dbl> <dbl> <chr> <chr> <chr> <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
1 1 1276973 rs70949… G GACAC 459. PASS 2 1.00 2 <NA> 7
2 1 3396278 . T TGCC… 447. PASS 1 0.50 2 <NA> 8
3 1 19208145 rs11254… G GGGT… 2735. PASS 2 1.00 2 <NA> 25
4 1 22207804 rs66803… ACCC… A 1769. PASS 2 1.00 2 <NA> 27
5 1 31764713 rs11166… C CAAG 1018. PASS 2 1.00 2 <NA> 16
6 1 31905889 rs30504… A ACAG 2172. PASS 2 1.00 2 <NA> 39
# ℹ 12 more variables: FS <chr>, HRun <chr>, HaplotypeScore <chr>, MQ0 <chr>,
# MQ <chr>, QD <dbl>, set <chr>, BaseQRankSum <chr>, Dels <chr>,
# MQRankSum <chr>, ReadPosRankSum <chr>, DS <chr>The GTF also has a column packed with key-value pairs.
my_gtf <- read_tsv(
file = "https://github.com/davetang/importbio/raw/master/inst/extdata/gencode.v38.annotation.sample.gtf.gz",
comment = "#",
show_col_types = FALSE,
col_names = c('chr', 'src', 'feat', 'start', 'end', 'score', 'strand', 'frame', 'group')
)
my_gtf |>
select(group) |>
head()# A tibble: 6 × 1
group
<chr>
1 "gene_id \"ENSG00000223972.5\"; gene_type \"transcribed_unprocessed_pseudogen…
2 "gene_id \"ENSG00000223972.5\"; transcript_id \"ENST00000456328.2\"; gene_typ…
3 "gene_id \"ENSG00000223972.5\"; transcript_id \"ENST00000456328.2\"; gene_typ…
4 "gene_id \"ENSG00000223972.5\"; transcript_id \"ENST00000456328.2\"; gene_typ…
5 "gene_id \"ENSG00000223972.5\"; transcript_id \"ENST00000456328.2\"; gene_typ…
6 "gene_id \"ENSG00000223972.5\"; transcript_id \"ENST00000450305.2\"; gene_typ…We can use the same strategy (but with some additional formatting steps) to split the column up.
my_gtf |>
mutate(group = sub(pattern = ";$", replacement = "", x = .data$group)) |>
mutate(group = gsub(pattern = '"', replacement = "", x = .data$group)) |>
separate_rows(group, sep = ";\\s*") |>
separate(group, c('key', 'value'), sep = "\\s") |>
distinct() |> # remove duplicated annotations, if any
pivot_wider(id_cols = chr:frame, names_from = key, values_from = value) -> my_gtf_splitWarning: Values from `value` are not uniquely identified; output will contain list-cols.
• Use `values_fn = list` to suppress this warning.
• Use `values_fn = {summary_fun}` to summarise duplicates.
• Use the following dplyr code to identify duplicates.
{data} |>
dplyr::summarise(n = dplyr::n(), .by = c(chr, src, feat, start, end, score,
strand, frame, key)) |>
dplyr::filter(n > 1L)head(my_gtf_split, 10)# A tibble: 10 × 24
chr src feat start end score strand frame gene_id gene_type gene_name
<chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <list> <list> <list>
1 chr1 HAVANA gene 11869 14409 . + . <chr> <chr [1]> <chr [1]>
2 chr1 HAVANA tran… 11869 14409 . + . <chr> <chr [1]> <chr [1]>
3 chr1 HAVANA exon 11869 12227 . + . <chr> <chr [1]> <chr [1]>
4 chr1 HAVANA exon 12613 12721 . + . <chr> <chr [1]> <chr [1]>
5 chr1 HAVANA exon 13221 14409 . + . <chr> <chr [1]> <chr [1]>
6 chr1 HAVANA tran… 12010 13670 . + . <chr> <chr [1]> <chr [1]>
7 chr1 HAVANA exon 12010 12057 . + . <chr> <chr [1]> <chr [1]>
8 chr1 HAVANA exon 12179 12227 . + . <chr> <chr [1]> <chr [1]>
9 chr1 HAVANA exon 12613 12697 . + . <chr> <chr [1]> <chr [1]>
10 chr1 HAVANA exon 12975 13052 . + . <chr> <chr [1]> <chr [1]>
# ℹ 13 more variables: level <list>, hgnc_id <list>, havana_gene <list>,
# transcript_id <list>, transcript_type <list>, transcript_name <list>,
# transcript_support_level <list>, tag <list>, havana_transcript <list>,
# exon_number <list>, exon_id <list>, ont <list>, protein_id <list>However, the split columns are lists because there were some cases
where there were multiple annotations with the same key and a list is
needed to store multiple values (which was what the warning above was
about). For example the tag key was repeated more than once
with different unique values for some annotations.
map_lgl(my_gtf_split$tag, function(x) length(x) > 1) [1] FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[13] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[25] TRUE FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[37] FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
[49] TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE
[61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE
[73] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE
[85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUEWe can check which columns have multiple values.
check_column <- function(x){
any(map_lgl(x, function(y) length(y) > 1))
}
my_check <- map_lgl(my_gtf_split, check_column)
my_check[my_check] transcript_id transcript_name tag havana_transcript
TRUE TRUE TRUE TRUE
exon_number ont
TRUE TRUE Therefore despite only a subset of the columns containing multiple
values, all the pivoted columns were turned into lists. However we can
turn the columns back into characters, which makes sense for the
gene_id column which only contained single unique character
values in the first place.
my_gtf_split |>
mutate(gene_id = as.character(gene_id)) |>
head()# A tibble: 6 × 24
chr src feat start end score strand frame gene_id gene_type gene_name
<chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr> <list> <list>
1 chr1 HAVANA gene 11869 14409 . + . ENSG00… <chr [1]> <chr [1]>
2 chr1 HAVANA trans… 11869 14409 . + . ENSG00… <chr [1]> <chr [1]>
3 chr1 HAVANA exon 11869 12227 . + . ENSG00… <chr [1]> <chr [1]>
4 chr1 HAVANA exon 12613 12721 . + . ENSG00… <chr [1]> <chr [1]>
5 chr1 HAVANA exon 13221 14409 . + . ENSG00… <chr [1]> <chr [1]>
6 chr1 HAVANA trans… 12010 13670 . + . ENSG00… <chr [1]> <chr [1]>
# ℹ 13 more variables: level <list>, hgnc_id <list>, havana_gene <list>,
# transcript_id <list>, transcript_type <list>, transcript_name <list>,
# transcript_support_level <list>, tag <list>, havana_transcript <list>,
# exon_number <list>, exon_id <list>, ont <list>, protein_id <list>But we can also do this to the tag column (even if it
needed a list to store the multiple values) and entries with multiple
values get turned into R (character) code!
my_gtf_split |>
mutate(tag = as.character(tag)) |>
select(tag) |>
head()# A tibble: 6 × 1
tag
<chr>
1 "NULL"
2 "basic"
3 "basic"
4 "basic"
5 "basic"
6 "c(\"basic\", \"Ensembl_canonical\")"If you don’t mind having R (character) code in your data, you can perform this transformation across all pivoted columns.
my_gtf_split |>
mutate(across(gene_id:protein_id, as.character))# A tibble: 94 × 24
chr src feat start end score strand frame gene_id gene_type gene_name
<chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
1 chr1 HAVANA gene 11869 14409 . + . ENSG00… transcri… DDX11L1
2 chr1 HAVANA tran… 11869 14409 . + . ENSG00… transcri… DDX11L1
3 chr1 HAVANA exon 11869 12227 . + . ENSG00… transcri… DDX11L1
4 chr1 HAVANA exon 12613 12721 . + . ENSG00… transcri… DDX11L1
5 chr1 HAVANA exon 13221 14409 . + . ENSG00… transcri… DDX11L1
6 chr1 HAVANA tran… 12010 13670 . + . ENSG00… transcri… DDX11L1
7 chr1 HAVANA exon 12010 12057 . + . ENSG00… transcri… DDX11L1
8 chr1 HAVANA exon 12179 12227 . + . ENSG00… transcri… DDX11L1
9 chr1 HAVANA exon 12613 12697 . + . ENSG00… transcri… DDX11L1
10 chr1 HAVANA exon 12975 13052 . + . ENSG00… transcri… DDX11L1
# ℹ 84 more rows
# ℹ 13 more variables: level <chr>, hgnc_id <chr>, havana_gene <chr>,
# transcript_id <chr>, transcript_type <chr>, transcript_name <chr>,
# transcript_support_level <chr>, tag <chr>, havana_transcript <chr>,
# exon_number <chr>, exon_id <chr>, ont <chr>, protein_id <chr>The following steps can be used to split a column containing key-value pairs into separate columns:
separate_rows to split a single column into
rowsseparate to split a key-value pair into two
columnspivot_wider to convert the long format table back
to wide formatHowever, sometimes data is packed into a single column because it cannot be nicely formatted in the first place.
sessionInfo()R version 4.4.0 (2024-04-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.4 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: Etc/UTC
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[5] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
[9] ggplot2_3.5.1 tidyverse_2.0.0 workflowr_1.7.1
loaded via a namespace (and not attached):
[1] sass_0.4.9 utf8_1.2.4 generics_0.1.3 stringi_1.8.4
[5] hms_1.1.3 digest_0.6.35 magrittr_2.0.3 timechange_0.3.0
[9] evaluate_0.24.0 grid_4.4.0 fastmap_1.2.0 rprojroot_2.0.4
[13] jsonlite_1.8.8 processx_3.8.4 whisker_0.4.1 ps_1.7.6
[17] promises_1.3.0 httr_1.4.7 fansi_1.0.6 scales_1.3.0
[21] jquerylib_0.1.4 cli_3.6.2 crayon_1.5.2 rlang_1.1.4
[25] bit64_4.0.5 munsell_0.5.1 withr_3.0.0 cachem_1.1.0
[29] yaml_2.3.8 parallel_4.4.0 tools_4.4.0 tzdb_0.4.0
[33] colorspace_2.1-0 httpuv_1.6.15 curl_5.2.1 vctrs_0.6.5
[37] R6_2.5.1 lifecycle_1.0.4 git2r_0.33.0 bit_4.0.5
[41] fs_1.6.4 vroom_1.6.5 pkgconfig_2.0.3 callr_3.7.6
[45] pillar_1.9.0 bslib_0.7.0 later_1.3.2 gtable_0.3.5
[49] glue_1.7.0 Rcpp_1.0.12 xfun_0.44 tidyselect_1.2.1
[53] rstudioapi_0.16.0 knitr_1.47 htmltools_0.5.8.1 rmarkdown_2.27
[57] compiler_4.4.0 getPass_0.2-4