In arendsee/monadR: A Monadic Pipeline System
This work is funded by the National Science Foundation grant NSF-IOS
1546858.
The problem
Gene Feature Format (GFF) is used to annotate intervals on a genome. Loading
and validating a GFF is a common first step in a bioinformatics analysis.
Mistakes at this step can cause major problems later, so it is important to
validate the GFF and report good diagnostics when problems arise. Here I will
show how such a pipeline can be written using standard methods, and then show
how rmonad can be used to organize and annotate such a pipeline.
GFF files are TAB-delimited where each row corresponds to a single interval.
These intervals, though, may be ontologically related. For example a gene is
the parent of an mRNA, which in turn is a parent of a set of exons and a set of
coding sequences (CDS). These relations are specified in the attribute column
(column 9). Here is an example (with TABs replaced with space) that introduces
the main format and the variants we need to support:
# this is a comment, they can appear anywhere in the GFF
# also note, empty lines can appear anywhere in the file
# This is a simple mono-exonic gene
I . gene 11565 11951 . - . ID=gene0
I . mRNA 11565 11951 . - . ID=mrna0;Parent=gene0
I . exon 11565 11951 . - . Parent=mrna0
I . CDS 11565 11951 . - 0 Parent=mrna0
# this is a gene with two splicing variants:
I . gene 61931 83591 . + . ID=gene1
I . mRNA 61931 83591 . + . ID=rna1;Parent=gene1
I . exon 61931 62344 . + . Parent=rna1
I . exon 81616 82209 . + . Parent=rna1
I . exon 82211 83591 . + . Parent=rna1
I . CDS 61931 62344 . + 0 Parent=rna1
I . CDS 81616 82209 . + 0 Parent=rna1
I . CDS 82211 82681 . + 0 Parent=rna1
I . mRNA 61931 83591 . + . ID=rna2;Parent=gene1
I . exon 61931 62344 . + . Parent=rna2
I . exon 81616 82209 . + . Parent=rna2
I . exon 82211 83591 . + . Parent=rna2
I . CDS 61931 62344 . + 0 Parent=rna2
I . CDS 82211 82681 . + 0 Parent=rna2
# Below are a few variants that occur (unforunately) in the wild
# V1: CDS directly descending from a gene.
b . gene 7235 9016 . - . ID=gene2
b . CDS 7235 9016 . - 0 Parent=gene2
# V2: 'Parent=-' when feature has no parent
I . gene 11565 11951 . - . ID=gene3;Parent=-
# V3: Elements with no tags that need to be treated as IDs
I . gene 11565 11951 . - . gene3
I . mRNA 11565 11951 . - . ID=mrna3;Parent=gene3
I . exon 11565 11951 . - . Parent=mrna3
I . CDS 11565 11951 . - 0 Parent=mrna3
Another issues we need to account for is type synonyms. The feature type
(column 3) is required to be valid Sequence Ontology (SO) terms. For the
purposes of this vignette, I will just handle the following sets of synonyms:
-
gene := gene | SO:0000704
-
mRNA := mRNA | messenger RNA | messenger_RNA | SO:0000234
| transcript | SO:0000673
-
CDS := CDS | coding_sequence | coding sequence | SO:0000316
-
exon := exon | SO:0000147 | coding_exon | coding exon |SO:0000195
For mRNA and exon, I am merging to ontology terms (mRNA and transcript; exon
and coding exon). Formally, this is incorrect, but practically it is probably
the right thing. Since these transformations may be wrong, they need to be
noted.
Solution
To test our solution, I use the gff rmonad dataset.
Conventional approach
Before using rmonad, I will use a more conventional approach.
library(readr)
library(dplyr)
library(tidyr)
library(magrittr)
set.seed(210)
good_gfffile <- system.file('extdata', 'gff', '0.gff3', package='rmonad')
read_gff <- function(file){
readr::read_tsv(
file,
col_names = c(
"seqid",
"source",
"type",
"start",
"stop",
"score",
"strand",
"phase",
"attr"
),
na = ".",
comment = "#",
col_types = "ccciidcic"
)
}
read_gff(good_gfffile)
This is a nice start, also, readr will pick up on any deviations from the
specified column number and type, warning of problems:
# wrong type
bad0 <- system.file('extdata', 'gff', 'bad-0.gff3', package='rmonad')
bad_result0 <- read_gff(bad0)
bad_result0
Notice that the 4th entry has a start value of "NA". This indicates there was
something wrong with the gff file and when we take a look at the original gff,
we see that the entry corresponds to the string "one thousand eight hundred and
seven" instead of the number. A more likely error would be including a comma in
a number (e.g., "1,807"), which R will parse as a string.
```{,eval=FALSE}
incorrect type
a Strex region 1 230218 . + . ID=I
a Strex telomere 1 801 . - . ID=id1
a Strex origin_of_replication 707 776 . + . ID=id2
a Strex gene "one thousand eight hundred and seven" 2169 . - . ID=gene1
Similarly when we load in a gff missing a column.
```r
# wrong column number
bad1 <- system.file('extdata', 'gff', 'bad-1.gff3', package='rmonad')
bad_result1 <- read_gff(bad1)
bad_result1
The read_gff can correctly identify that a column is missing from the following raw file.
```{,eval=FALSE}
missing one column
a region 1 230218 . + . ID=I
a telomere 1 801 . - . ID=id1
a origin_of_replication 707 776 . + . ID=id2
a gene 1807 2169 . - . ID=gene1
R uses `NA` to indicate missing values. The R 'numeric' type corresponds to
Haskell `[Maybe Num]`, i.e. an array of possibly empty values. In a GFF,
columns 2,6,7,8 and 9 may be missing, the others may not. So we need an
additional assertion that these are complete.
```r
g <- read_gff(good_gfffile)
for(col in c("seqid", "type", "start", "stop")){
if(any(is.na(g[[col]]))){
stop("GFFError: Column '", col, "' may not have missing values")
}
}
Now we can account for type synonyms using the following method
gene_synonyms <- 'SO:0000704'
mRNA_synonyms <- c('messenger_RNA', 'messenger RNA', 'SO:0000234')
CDS_synonyms <- c('coding_sequence', 'coding sequence', 'SO:0000316')
exon_synonyms <- 'SO:0000147'
g$type <- ifelse(g$type %in% gene_synonyms, 'gene', g$type)
g$type <- ifelse(g$type %in% mRNA_synonyms, 'mRNA', g$type)
g$type <- ifelse(g$type %in% CDS_synonyms, 'CDS', g$type)
g$type <- ifelse(g$type %in% exon_synonyms, 'exon', g$type)
mRNA_near_synonyms <- c('transcript', 'SO:0000673')
exon_near_synonyms <- c('SO:0000147', 'coding_exon', 'coding exon', 'SO:0000195')
if(any(g$type %in% mRNA_near_synonyms)){
g$type <- ifelse(g$type %in% mRNA_near_synonyms, 'mRNA', g$type)
warning("Substituting transcript types for mRNA types, this is probably OK")
}
if(any(g$type %in% exon_near_synonyms)){
g$type <- ifelse(g$type %in% exon_near_synonyms, 'exon', g$type)
warning("Substituting transcript types for exon types, this is probably OK")
}
Now we need to evaluate the attribute column.
tags <- c("ID", "Parent")
data_frame(
attr = stringr::str_split(g$attr, ";"),
order = 1:nrow(g)
) %>%
dplyr::mutate(ntags = sapply(attr, length)) %>%
tidyr::unnest(attr) %>%
dplyr::mutate(attr = ifelse(grepl('=', attr), attr, paste(".U", attr, sep="="))) %>%
tidyr::separate_(
col = "attr",
into = c("tag", "value"),
sep = "=",
extra = "merge"
) %>%
dplyr::filter(tag %in% c(tags, ".U")) %>%
{
if(nrow(.) > 0){
tidyr::spread(., key="tag", value="value")
} else {
.$tag = NULL
.$value = NULL
.
}
} %>%
{
if("Parent" %in% names(.)){
.$Parent <- ifelse(.$Parent == "-", NA, .$Parent)
}
.
} %>% {
for(tag in c(tags, ".U")){
if(! tag %in% names(.))
.[[tag]] = NA_character_
}
.
} %>%
{
if("ID" %in% names(.))
.$ID <- ifelse(is.na(.$ID) & !is.na(.$.U) & .$ntags == 1, .$.U, .$ID)
.
} %>%
merge(data_frame(order=1:nrow(g)), all=TRUE) %>%
dplyr::arrange(order) %>%
{ cbind(g, .) } %>%
dplyr::select(-.U, -order, -ntags, -attr) %>%
{
if(all(c("ID", "Parent") %in% names(.))){
parents <- subset(., type %in% c("CDS", "exon"))$Parent
parent_types <- subset(., ID %in% parents)$type
if(any(parent_types == "gene"))
warning("Found CDS or exon directly inheriting from a gene, this may be fine.")
if(! all(parent_types %in% c("gene", "mRNA")))
stop("Found CDS or exon with illegal parent")
if( any(is.na(parents)) )
stop("Found CDS or exon with no parent")
if(! any(duplicated(.$ID, incomparables=NA)))
warning("IDs are not unique, this is probably bad")
}
.
}
The beauty of this chain is that it requires few temporary variables (just g,
and tags), it is a pure flow of data. It is an elegant sequence of functions
operating on a single thread of data.
But there are a few problems.
First, it is in dire need of documentation. We could add in comments. But
comments cannot be formatted well. A better approach is some form of literate
programming, such as rewriting the program in Rmarkdown. But this 1) breaks the
pipeline (since we can't pipe between chunks), 2) results in an object we can't
compute on, 3) makes debugging even more difficult, because our code is spread
out.
rmonad approach
With rmonad we can mingle documentation and code in a computable object.
require(readr)
require(dplyr)
require(tidyr)
require(magrittr)
library(rmonad)
set.seed(210)
#' @param file The input GFF file
#' @param tags The GFF tags to keep
#' default these are not stored)
read_gff <- function(file, tags){
"
This function reads a GFF. It builds a data frame from it that with a column
for each value in `tags`.
The `read_gff` function is a nested function, containing an Rmonad pipeline
inside it. This docstring describes the entire wrapped pipeline.
"
raw_gff <- evalwrap(
{
"
Rmonad supports docstrings. If a block begins with a string, this
string is extracted and stored. Python has something similar, where the
first string in a function is cast as documentation.
The `evalwrap` function takes an expression and wraps its result into a
context. It also handles the extraction of this docstring. The result
here is used at more than one place in the pipeline. Rather than
accessing it later as a global, it will be funneled bach in.
"
readr::read_tsv(
file,
col_names = c(
"seqid",
"source",
"type",
"start",
"stop",
"score",
"strand",
"phase",
"attr"
),
na = ".",
comment = "#",
col_types = "ccciidcic"
)
}
)
raw_gff %>% tag('raw_gff') %>>% {
"
The %>>% operator applies the function described in this block to the
input on the left-hand-side. This corresponds to the UNIX '|' or magrittr's
'%>%'. It differs from them in that it is a monadic bind operator, rather
than an application operator. It carries a context along with the
computations. The context can store past values, performance information,
this docstring, and links to the parent chunk. The context is a directed
graph of code chunks and their metadata.
Here I assert that no data is missing in the columns where it is required.
"
for(col in c("seqid", "type", "start", "stop")){
if(any(is.na(.[[col]]))){
stop("GFFError: Column '", col, "' may not have missing values")
}
}
.
} %>>% {
"
According the GFF specification, the value in the type column should be
from the Sequence Ontology (SO). Here I collapse the synonymous terms
for the main features usually present in GFF file (gene, mRNA, CDS, and
exon). In production code, I would probably do something a bit more formal.
"
gene_synonyms <- 'SO:0000704'
mRNA_synonyms <- c('messenger_RNA', 'messenger RNA', 'SO:0000234')
CDS_synonyms <- c('coding_sequence', 'coding sequence', 'SO:0000316')
exon_synonyms <- 'SO:0000147'
.$type <- ifelse(.$type %in% gene_synonyms, 'gene', .$type)
.$type <- ifelse(.$type %in% mRNA_synonyms, 'mRNA', .$type)
.$type <- ifelse(.$type %in% CDS_synonyms, 'CDS', .$type)
.$type <- ifelse(.$type %in% exon_synonyms, 'exon', .$type)
.
} %>_% {
"
The %>_% operator lets this chunk of code be run for its effects, which are
raising warnings if we replace the type with a questionable synonym. We
could alternatively just use %>>% and add a terminal '.' to this chunk.
However, the `%>_%` operator guarantees the input is passed without being
changed, which is safer and also clarifies our intention.
"
mRNA_near_synonyms <- c('transcript', 'SO:0000673')
exon_near_synonyms <- c('SO:0000147', 'coding_exon', 'coding exon', 'SO:0000195')
if(any(.$type %in% mRNA_near_synonyms)){
.$type <- ifelse(.$type %in% mRNA_near_synonyms, 'mRNA', .$type)
warning("Substituting transcript types for mRNA types, this is probably OK")
}
if(any(.$type %in% exon_near_synonyms)){
.$type <- ifelse(.$type %in% exon_near_synonyms, 'exon', .$type)
warning("Substituting transcript types for exon types, this is probably OK")
}
} %>>% {
"
Notice below that I use the magrittr operator '%>%' inside the rmonad
pipeline. When to pipe with rmonad and when to pipe with magrittr is a
matter of granularity. This chunk of code perhaps should form one
documentation unit. And perhaps I don't expect it to fail. If I break this
chunk into several, the failures are more localized, and I can access
intermediate values for debugging. On the other hand, putting every little
operation in a new chunk will clutter the graph and reports derived from
it.
"
data_frame(
attr = stringr::str_split(.$attr, ";"),
order = 1:nrow(.)
) %>%
dplyr::mutate(ntags = sapply(attr, length)) %>%
tidyr::unnest(attr) %>%
dplyr::mutate(attr = ifelse(grepl('=', attr), attr, paste(".U", attr, sep="="))) %>%
tidyr::separate_(
col = "attr",
into = c("tag", "value"),
sep = "=",
extra = "merge"
)
} %>% funnel(raw_gff=raw_gff, tags=tags) %*>% {
"
The %v>% operator stores the input value in memory. We could replace every
%>>% operator with %v>%. This would let us inspect every step of an
analysis at the cost of high memory usage. For brevity, I won't break this
following block down any further.
The `funnel` function packages a list in a monad, merging their histories
and propagating error. That is, if `gff` or `tags` failed upstream, this
function will not be run. `%*>%` takes a list on the left and feeds it into
the function on the right as an argument list. Here `funnel` and `%*>%` are
used together to merge a pipeline (gff) and inject a parameter (tags).
We could not have written
%v>% function(gff=gff, tags=tags)
because this would have brought the monad wrapped gff into scope, not the
value itself.
"
dplyr::filter(., tag %in% c(tags, ".U")) %>%
{
if(nrow(.) > 0){
tidyr::spread(., key="tag", value="value")
} else {
.$tag = NULL
.$value = NULL
.
}
} %>%
{
if("Parent" %in% names(.)){
.$Parent <- ifelse(.$Parent == "-", NA, .$Parent)
}
.
} %>% {
for(tag in c(tags, ".U")){
if(! tag %in% names(.))
.[[tag]] = NA_character_
}
.
} %>%
{
if("ID" %in% names(.))
.$ID <- ifelse(is.na(.$ID) & !is.na(.$.U) & .$ntags == 1, .$.U, .$ID)
.
} %>%
merge(data_frame(order=1:nrow(raw_gff)), all=TRUE) %>%
dplyr::arrange(order) %>%
{ cbind(raw_gff, .) } %>%
dplyr::select(-.U, -order, -ntags, -attr)
} %>_% {
"
And make the last few assertions.
"
if(all(c("ID", "Parent") %in% names(.))){
parents <- subset(., type %in% c("CDS", "exon"))$Parent
parent_types <- subset(., ID %in% parents)$type
if(any(parent_types == "gene"))
warning("Found CDS or exon directly inheriting from a gene, this may be fine.")
if(! all(parent_types %in% c("gene", "mRNA")))
stop("Found CDS or exon with illegal parent")
if(any(is.na(parents)) )
stop("Found CDS or exon with no parent")
if(any(duplicated(.$ID, incomparables=NA)))
warning("IDs are not unique, this is probably bad")
}
} %>_% {
"
I could post some closing comments here. The %>_% operator can be chained
and the output does not affect the output of the main chain. The NULL is
required to distinguish this block from an anonymous function that returns
a string.
"
NULL
}
# End Rmonad chain
}
gff_as_graph <- function(file){
read_gff(file, tags=c("ID", "Parent")) %>>%
dplyr::filter(!is.na(Parent)) %>>%
{
"Make a unique label for each CDS and exon interval"
dplyr::group_by(., Parent) %>%
dplyr::mutate(ID = ifelse(type %in% c("CDS", "exon"),
paste(type, start, sep='-'),
ID)) %>%
dplyr::ungroup()
} %>>% {
"Make a graph of gene models"
dplyr::select(., ID, Parent) %>%
as.matrix %>%
igraph::graph_from_edgelist()
}
}
That is the whole GFF program in an rmonad framework
result <-
system.file('extdata', 'gff', '0.gff3', package='rmonad') %v>% gff_as_graph
result %>% plot
esc will extract the final result, and raise all errors, warnings and
messages that where extracted
esc(result)
Now we see why we might want a little more granularity in our pipeline. To
summarize the results we can use the mtabulate functions.
mtabulate(result)
We can also get a summary of issues
missues(result)
The id column corresponds to a row number in the mtabulate result.
To extract particular values, we can use the get_* family of vectorized
getters
# get a list of every stored value, report uncached values as NULL
get_value(result, warn=FALSE)
# get a list of every docstring
get_doc(result)
plot(result)
A plot of the connections between functions in the pipeline. The node labels
are the rmonad ids. The red arrow indicates a 'nest' relationship. Green nodes
are OK, orange nodes raised a warning, red nodes (none appear in this graph)
represent errors.
If the pipeline fails, the last valid result is saved, for example
m <-
system.file('extdata', 'gff', '15.gff3', package='rmonad') %v>%
gff_as_graph
To diagnose the problem, we might first read the error messages:
get_error(m)
Here we see the 4th node raised an error: we having missing values in the
'type' column.
To get a big picture of the pipeline:
mtabulate(m)
The the last three nodes are in the failing state (see column 'OK'). We can see
all issues:
missues(m)
We see two warnings in node 3. To see what code is wrapped in node 3:
get_code(m, 3)
So there were issues in reading the file. We can get the value of this node's parent:
get_value(m, get_parents(m, 3)[[1]])[[1]]
If we wanted to explore the raw data, we could jump into a Bash shell, navigate
to the file above, and take it apart with UNIX tools (the lines of raw data
were not stored in the Rmonad object).
We might next check out the failing columns:
get_code(m, 4:6)
Now that we know exactly where the problem is, we can start working through
the data that triggered the problem:
get_value(m, get_parents(m, 5)[[1]])
get_value(m, 5)[[1]]
From this it looks like there is raw sequence data contained in the GFF file.
According to the GFF specification, it is legal to have the full sequence data
appended to the tabular GFF data. So we probably should either add a test that
will raise an appropriate error message or extend our function to handle this
case.
rmonad saved the input to the failing function. This makes debugging much
simpler. We did not have to disassemble our pipeline and rerun it step-by-step.
Overall, in rmonad, the output of a pipeline is not just the effluent of the
last node, but the collection of all the nodes along the way. The pipeline
itself becomes data that can be computed upon.
# save data, this is used in testing and examples
gff <- list(good_result=result, bad_result0 = bad_result0, bad_result1 = bad_result1)
save(gff, file="gff.rda", version=2)
arendsee/monadR documentation built on Dec. 16, 2020, 4:26 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
This work is funded by the National Science Foundation grant NSF-IOS 1546858.
The problem
Gene Feature Format (GFF) is used to annotate intervals on a genome. Loading
and validating a GFF is a common first step in a bioinformatics analysis.
Mistakes at this step can cause major problems later, so it is important to
validate the GFF and report good diagnostics when problems arise. Here I will
show how such a pipeline can be written using standard methods, and then show
how rmonad can be used to organize and annotate such a pipeline.
GFF files are TAB-delimited where each row corresponds to a single interval. These intervals, though, may be ontologically related. For example a gene is the parent of an mRNA, which in turn is a parent of a set of exons and a set of coding sequences (CDS). These relations are specified in the attribute column (column 9). Here is an example (with TABs replaced with space) that introduces the main format and the variants we need to support:
# this is a comment, they can appear anywhere in the GFF # also note, empty lines can appear anywhere in the file # This is a simple mono-exonic gene I . gene 11565 11951 . - . ID=gene0 I . mRNA 11565 11951 . - . ID=mrna0;Parent=gene0 I . exon 11565 11951 . - . Parent=mrna0 I . CDS 11565 11951 . - 0 Parent=mrna0 # this is a gene with two splicing variants: I . gene 61931 83591 . + . ID=gene1 I . mRNA 61931 83591 . + . ID=rna1;Parent=gene1 I . exon 61931 62344 . + . Parent=rna1 I . exon 81616 82209 . + . Parent=rna1 I . exon 82211 83591 . + . Parent=rna1 I . CDS 61931 62344 . + 0 Parent=rna1 I . CDS 81616 82209 . + 0 Parent=rna1 I . CDS 82211 82681 . + 0 Parent=rna1 I . mRNA 61931 83591 . + . ID=rna2;Parent=gene1 I . exon 61931 62344 . + . Parent=rna2 I . exon 81616 82209 . + . Parent=rna2 I . exon 82211 83591 . + . Parent=rna2 I . CDS 61931 62344 . + 0 Parent=rna2 I . CDS 82211 82681 . + 0 Parent=rna2 # Below are a few variants that occur (unforunately) in the wild # V1: CDS directly descending from a gene. b . gene 7235 9016 . - . ID=gene2 b . CDS 7235 9016 . - 0 Parent=gene2 # V2: 'Parent=-' when feature has no parent I . gene 11565 11951 . - . ID=gene3;Parent=- # V3: Elements with no tags that need to be treated as IDs I . gene 11565 11951 . - . gene3 I . mRNA 11565 11951 . - . ID=mrna3;Parent=gene3 I . exon 11565 11951 . - . Parent=mrna3 I . CDS 11565 11951 . - 0 Parent=mrna3
Another issues we need to account for is type synonyms. The feature type (column 3) is required to be valid Sequence Ontology (SO) terms. For the purposes of this vignette, I will just handle the following sets of synonyms:
-
gene :=
gene|SO:0000704 -
mRNA :=
mRNA|messenger RNA|messenger_RNA|SO:0000234|transcript|SO:0000673 -
CDS :=
CDS|coding_sequence|coding sequence|SO:0000316 -
exon :=
exon|SO:0000147|coding_exon|coding exon|SO:0000195
For mRNA and exon, I am merging to ontology terms (mRNA and transcript; exon and coding exon). Formally, this is incorrect, but practically it is probably the right thing. Since these transformations may be wrong, they need to be noted.
Solution
To test our solution, I use the gff rmonad dataset.
Conventional approach
Before using rmonad, I will use a more conventional approach.
library(readr) library(dplyr) library(tidyr) library(magrittr) set.seed(210) good_gfffile <- system.file('extdata', 'gff', '0.gff3', package='rmonad') read_gff <- function(file){ readr::read_tsv( file, col_names = c( "seqid", "source", "type", "start", "stop", "score", "strand", "phase", "attr" ), na = ".", comment = "#", col_types = "ccciidcic" ) } read_gff(good_gfffile)
This is a nice start, also, readr will pick up on any deviations from the
specified column number and type, warning of problems:
# wrong type bad0 <- system.file('extdata', 'gff', 'bad-0.gff3', package='rmonad') bad_result0 <- read_gff(bad0) bad_result0
Notice that the 4th entry has a start value of "NA". This indicates there was something wrong with the gff file and when we take a look at the original gff, we see that the entry corresponds to the string "one thousand eight hundred and seven" instead of the number. A more likely error would be including a comma in a number (e.g., "1,807"), which R will parse as a string.
```{,eval=FALSE}
incorrect type
a Strex region 1 230218 . + . ID=I a Strex telomere 1 801 . - . ID=id1 a Strex origin_of_replication 707 776 . + . ID=id2 a Strex gene "one thousand eight hundred and seven" 2169 . - . ID=gene1
Similarly when we load in a gff missing a column. ```r # wrong column number bad1 <- system.file('extdata', 'gff', 'bad-1.gff3', package='rmonad') bad_result1 <- read_gff(bad1) bad_result1
The read_gff can correctly identify that a column is missing from the following raw file. ```{,eval=FALSE}
missing one column
a region 1 230218 . + . ID=I a telomere 1 801 . - . ID=id1 a origin_of_replication 707 776 . + . ID=id2 a gene 1807 2169 . - . ID=gene1
R uses `NA` to indicate missing values. The R 'numeric' type corresponds to Haskell `[Maybe Num]`, i.e. an array of possibly empty values. In a GFF, columns 2,6,7,8 and 9 may be missing, the others may not. So we need an additional assertion that these are complete. ```r g <- read_gff(good_gfffile) for(col in c("seqid", "type", "start", "stop")){ if(any(is.na(g[[col]]))){ stop("GFFError: Column '", col, "' may not have missing values") } }
Now we can account for type synonyms using the following method
gene_synonyms <- 'SO:0000704' mRNA_synonyms <- c('messenger_RNA', 'messenger RNA', 'SO:0000234') CDS_synonyms <- c('coding_sequence', 'coding sequence', 'SO:0000316') exon_synonyms <- 'SO:0000147' g$type <- ifelse(g$type %in% gene_synonyms, 'gene', g$type) g$type <- ifelse(g$type %in% mRNA_synonyms, 'mRNA', g$type) g$type <- ifelse(g$type %in% CDS_synonyms, 'CDS', g$type) g$type <- ifelse(g$type %in% exon_synonyms, 'exon', g$type) mRNA_near_synonyms <- c('transcript', 'SO:0000673') exon_near_synonyms <- c('SO:0000147', 'coding_exon', 'coding exon', 'SO:0000195') if(any(g$type %in% mRNA_near_synonyms)){ g$type <- ifelse(g$type %in% mRNA_near_synonyms, 'mRNA', g$type) warning("Substituting transcript types for mRNA types, this is probably OK") } if(any(g$type %in% exon_near_synonyms)){ g$type <- ifelse(g$type %in% exon_near_synonyms, 'exon', g$type) warning("Substituting transcript types for exon types, this is probably OK") }
Now we need to evaluate the attribute column.
tags <- c("ID", "Parent") data_frame( attr = stringr::str_split(g$attr, ";"), order = 1:nrow(g) ) %>% dplyr::mutate(ntags = sapply(attr, length)) %>% tidyr::unnest(attr) %>% dplyr::mutate(attr = ifelse(grepl('=', attr), attr, paste(".U", attr, sep="="))) %>% tidyr::separate_( col = "attr", into = c("tag", "value"), sep = "=", extra = "merge" ) %>% dplyr::filter(tag %in% c(tags, ".U")) %>% { if(nrow(.) > 0){ tidyr::spread(., key="tag", value="value") } else { .$tag = NULL .$value = NULL . } } %>% { if("Parent" %in% names(.)){ .$Parent <- ifelse(.$Parent == "-", NA, .$Parent) } . } %>% { for(tag in c(tags, ".U")){ if(! tag %in% names(.)) .[[tag]] = NA_character_ } . } %>% { if("ID" %in% names(.)) .$ID <- ifelse(is.na(.$ID) & !is.na(.$.U) & .$ntags == 1, .$.U, .$ID) . } %>% merge(data_frame(order=1:nrow(g)), all=TRUE) %>% dplyr::arrange(order) %>% { cbind(g, .) } %>% dplyr::select(-.U, -order, -ntags, -attr) %>% { if(all(c("ID", "Parent") %in% names(.))){ parents <- subset(., type %in% c("CDS", "exon"))$Parent parent_types <- subset(., ID %in% parents)$type if(any(parent_types == "gene")) warning("Found CDS or exon directly inheriting from a gene, this may be fine.") if(! all(parent_types %in% c("gene", "mRNA"))) stop("Found CDS or exon with illegal parent") if( any(is.na(parents)) ) stop("Found CDS or exon with no parent") if(! any(duplicated(.$ID, incomparables=NA))) warning("IDs are not unique, this is probably bad") } . }
The beauty of this chain is that it requires few temporary variables (just g,
and tags), it is a pure flow of data. It is an elegant sequence of functions
operating on a single thread of data.
But there are a few problems.
First, it is in dire need of documentation. We could add in comments. But comments cannot be formatted well. A better approach is some form of literate programming, such as rewriting the program in Rmarkdown. But this 1) breaks the pipeline (since we can't pipe between chunks), 2) results in an object we can't compute on, 3) makes debugging even more difficult, because our code is spread out.
rmonad approach
With rmonad we can mingle documentation and code in a computable object.
require(readr) require(dplyr) require(tidyr) require(magrittr) library(rmonad) set.seed(210) #' @param file The input GFF file #' @param tags The GFF tags to keep #' default these are not stored) read_gff <- function(file, tags){ " This function reads a GFF. It builds a data frame from it that with a column for each value in `tags`. The `read_gff` function is a nested function, containing an Rmonad pipeline inside it. This docstring describes the entire wrapped pipeline. " raw_gff <- evalwrap( { " Rmonad supports docstrings. If a block begins with a string, this string is extracted and stored. Python has something similar, where the first string in a function is cast as documentation. The `evalwrap` function takes an expression and wraps its result into a context. It also handles the extraction of this docstring. The result here is used at more than one place in the pipeline. Rather than accessing it later as a global, it will be funneled bach in. " readr::read_tsv( file, col_names = c( "seqid", "source", "type", "start", "stop", "score", "strand", "phase", "attr" ), na = ".", comment = "#", col_types = "ccciidcic" ) } ) raw_gff %>% tag('raw_gff') %>>% { " The %>>% operator applies the function described in this block to the input on the left-hand-side. This corresponds to the UNIX '|' or magrittr's '%>%'. It differs from them in that it is a monadic bind operator, rather than an application operator. It carries a context along with the computations. The context can store past values, performance information, this docstring, and links to the parent chunk. The context is a directed graph of code chunks and their metadata. Here I assert that no data is missing in the columns where it is required. " for(col in c("seqid", "type", "start", "stop")){ if(any(is.na(.[[col]]))){ stop("GFFError: Column '", col, "' may not have missing values") } } . } %>>% { " According the GFF specification, the value in the type column should be from the Sequence Ontology (SO). Here I collapse the synonymous terms for the main features usually present in GFF file (gene, mRNA, CDS, and exon). In production code, I would probably do something a bit more formal. " gene_synonyms <- 'SO:0000704' mRNA_synonyms <- c('messenger_RNA', 'messenger RNA', 'SO:0000234') CDS_synonyms <- c('coding_sequence', 'coding sequence', 'SO:0000316') exon_synonyms <- 'SO:0000147' .$type <- ifelse(.$type %in% gene_synonyms, 'gene', .$type) .$type <- ifelse(.$type %in% mRNA_synonyms, 'mRNA', .$type) .$type <- ifelse(.$type %in% CDS_synonyms, 'CDS', .$type) .$type <- ifelse(.$type %in% exon_synonyms, 'exon', .$type) . } %>_% { " The %>_% operator lets this chunk of code be run for its effects, which are raising warnings if we replace the type with a questionable synonym. We could alternatively just use %>>% and add a terminal '.' to this chunk. However, the `%>_%` operator guarantees the input is passed without being changed, which is safer and also clarifies our intention. " mRNA_near_synonyms <- c('transcript', 'SO:0000673') exon_near_synonyms <- c('SO:0000147', 'coding_exon', 'coding exon', 'SO:0000195') if(any(.$type %in% mRNA_near_synonyms)){ .$type <- ifelse(.$type %in% mRNA_near_synonyms, 'mRNA', .$type) warning("Substituting transcript types for mRNA types, this is probably OK") } if(any(.$type %in% exon_near_synonyms)){ .$type <- ifelse(.$type %in% exon_near_synonyms, 'exon', .$type) warning("Substituting transcript types for exon types, this is probably OK") } } %>>% { " Notice below that I use the magrittr operator '%>%' inside the rmonad pipeline. When to pipe with rmonad and when to pipe with magrittr is a matter of granularity. This chunk of code perhaps should form one documentation unit. And perhaps I don't expect it to fail. If I break this chunk into several, the failures are more localized, and I can access intermediate values for debugging. On the other hand, putting every little operation in a new chunk will clutter the graph and reports derived from it. " data_frame( attr = stringr::str_split(.$attr, ";"), order = 1:nrow(.) ) %>% dplyr::mutate(ntags = sapply(attr, length)) %>% tidyr::unnest(attr) %>% dplyr::mutate(attr = ifelse(grepl('=', attr), attr, paste(".U", attr, sep="="))) %>% tidyr::separate_( col = "attr", into = c("tag", "value"), sep = "=", extra = "merge" ) } %>% funnel(raw_gff=raw_gff, tags=tags) %*>% { " The %v>% operator stores the input value in memory. We could replace every %>>% operator with %v>%. This would let us inspect every step of an analysis at the cost of high memory usage. For brevity, I won't break this following block down any further. The `funnel` function packages a list in a monad, merging their histories and propagating error. That is, if `gff` or `tags` failed upstream, this function will not be run. `%*>%` takes a list on the left and feeds it into the function on the right as an argument list. Here `funnel` and `%*>%` are used together to merge a pipeline (gff) and inject a parameter (tags). We could not have written %v>% function(gff=gff, tags=tags) because this would have brought the monad wrapped gff into scope, not the value itself. " dplyr::filter(., tag %in% c(tags, ".U")) %>% { if(nrow(.) > 0){ tidyr::spread(., key="tag", value="value") } else { .$tag = NULL .$value = NULL . } } %>% { if("Parent" %in% names(.)){ .$Parent <- ifelse(.$Parent == "-", NA, .$Parent) } . } %>% { for(tag in c(tags, ".U")){ if(! tag %in% names(.)) .[[tag]] = NA_character_ } . } %>% { if("ID" %in% names(.)) .$ID <- ifelse(is.na(.$ID) & !is.na(.$.U) & .$ntags == 1, .$.U, .$ID) . } %>% merge(data_frame(order=1:nrow(raw_gff)), all=TRUE) %>% dplyr::arrange(order) %>% { cbind(raw_gff, .) } %>% dplyr::select(-.U, -order, -ntags, -attr) } %>_% { " And make the last few assertions. " if(all(c("ID", "Parent") %in% names(.))){ parents <- subset(., type %in% c("CDS", "exon"))$Parent parent_types <- subset(., ID %in% parents)$type if(any(parent_types == "gene")) warning("Found CDS or exon directly inheriting from a gene, this may be fine.") if(! all(parent_types %in% c("gene", "mRNA"))) stop("Found CDS or exon with illegal parent") if(any(is.na(parents)) ) stop("Found CDS or exon with no parent") if(any(duplicated(.$ID, incomparables=NA))) warning("IDs are not unique, this is probably bad") } } %>_% { " I could post some closing comments here. The %>_% operator can be chained and the output does not affect the output of the main chain. The NULL is required to distinguish this block from an anonymous function that returns a string. " NULL } # End Rmonad chain } gff_as_graph <- function(file){ read_gff(file, tags=c("ID", "Parent")) %>>% dplyr::filter(!is.na(Parent)) %>>% { "Make a unique label for each CDS and exon interval" dplyr::group_by(., Parent) %>% dplyr::mutate(ID = ifelse(type %in% c("CDS", "exon"), paste(type, start, sep='-'), ID)) %>% dplyr::ungroup() } %>>% { "Make a graph of gene models" dplyr::select(., ID, Parent) %>% as.matrix %>% igraph::graph_from_edgelist() } }
That is the whole GFF program in an rmonad framework
result <- system.file('extdata', 'gff', '0.gff3', package='rmonad') %v>% gff_as_graph result %>% plot
esc will extract the final result, and raise all errors, warnings and
messages that where extracted
esc(result)
Now we see why we might want a little more granularity in our pipeline. To
summarize the results we can use the mtabulate functions.
mtabulate(result)
We can also get a summary of issues
missues(result)
The id column corresponds to a row number in the mtabulate result.
To extract particular values, we can use the get_* family of vectorized
getters
# get a list of every stored value, report uncached values as NULL get_value(result, warn=FALSE) # get a list of every docstring get_doc(result)
plot(result)
A plot of the connections between functions in the pipeline. The node labels
are the rmonad ids. The red arrow indicates a 'nest' relationship. Green nodes
are OK, orange nodes raised a warning, red nodes (none appear in this graph)
represent errors.
If the pipeline fails, the last valid result is saved, for example
m <- system.file('extdata', 'gff', '15.gff3', package='rmonad') %v>% gff_as_graph
To diagnose the problem, we might first read the error messages:
get_error(m)
Here we see the 4th node raised an error: we having missing values in the 'type' column.
To get a big picture of the pipeline:
mtabulate(m)
The the last three nodes are in the failing state (see column 'OK'). We can see all issues:
missues(m)
We see two warnings in node 3. To see what code is wrapped in node 3:
get_code(m, 3)
So there were issues in reading the file. We can get the value of this node's parent:
get_value(m, get_parents(m, 3)[[1]])[[1]]
If we wanted to explore the raw data, we could jump into a Bash shell, navigate to the file above, and take it apart with UNIX tools (the lines of raw data were not stored in the Rmonad object).
We might next check out the failing columns:
get_code(m, 4:6)
Now that we know exactly where the problem is, we can start working through the data that triggered the problem:
get_value(m, get_parents(m, 5)[[1]]) get_value(m, 5)[[1]]
From this it looks like there is raw sequence data contained in the GFF file. According to the GFF specification, it is legal to have the full sequence data appended to the tabular GFF data. So we probably should either add a test that will raise an appropriate error message or extend our function to handle this case.
rmonad saved the input to the failing function. This makes debugging much
simpler. We did not have to disassemble our pipeline and rerun it step-by-step.
Overall, in rmonad, the output of a pipeline is not just the effluent of the
last node, but the collection of all the nodes along the way. The pipeline
itself becomes data that can be computed upon.
# save data, this is used in testing and examples gff <- list(good_result=result, bad_result0 = bad_result0, bad_result1 = bad_result1) save(gff, file="gff.rda", version=2)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.