Bioconductor Newsletter

## [1] "Arabidopsis_thaliana"     "Caenorhabditis_elegans"  
## [3] "Cyanidioschyzon_merolae"  "Drosophila_melanogaster" 
## [5] "Homo_sapiens"             "Oryza_sativa"            
## [7] "Populus_trichocarpa"      "Saccharomyces_cerevisiae"
## [9] "Zea_mays"

## [1] "UCSC"

##  chrII chrIII   chrM 
##   "II"  "III"   "MT"

## GRanges with 3 ranges and 0 metadata columns:
##       seqnames    ranges strand
##             
##   [1]        A   [1, 10]      *
##   [2]        A   [5, 15]      *
##   [3]        A   [5, 10]      *
##   ---
##   seqlengths:
##     A
##    NA




First we disjoin the ranges and perform `findOverlaps` which
will give us a `Hits` object:

gr2 <- disjoin(gr1)
ov <- findOverlaps(gr1, gr2)
ov

## Hits of length 5
## queryLength: 3
## subjectLength: 3
##   queryHits subjectHits 
##        
##  1         1           1 
##  2         1           2 
##  3         2           2 
##  4         2           3 
##  5         3           2




Convenience functions `countQueryHits' and `countSubjectHits'
summarize the number of times each record in the `query` or
`subject` was hit. 

countSubjectHits(ov)

## [1] 1 3 1




The regions that are not shared between the ranges will be
counted only once (self-hit).

keep_idx <- which(countSubjectHits(ov) == 1L)




Accessors `queryHits` and `subjectHits` extract the corresponding
columns from the `Hits` object:

subjectHits(ov)

## [1] 1 2 2 3 2




The extracted `subjectHits` can be used in a matching operation
to subset the `Hits` object containing just the "unshared" regions:

ov <- ov[subjectHits(ov) %in% keep_idx]
ov

## Hits of length 2
## queryLength: 3
## subjectLength: 3
##   queryHits subjectHits 
##        
##  1         1           1 
##  2         2           3




Again `subjectHits` are extracted but this time from the subset
`Hits` object. We use this index to subset the disjoined original
ranges leaving only the "unshared" regions:

ans_regions <- gr2[subjectHits(ov)]
ans_regions

## GRanges with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##             
##   [1]        A  [ 1,  4]      *
##   [2]        A  [11, 15]      *
##   ---
##   seqlengths:
##     A
##    NA




To finish this off we make a `GRangesList` of the unshared regions.
The list is length 3 corresponding the the 3 ranges in `gr1':

ans_eltlens <- countQueryHits(ov)
ans_skeleton <- PartitioningByEnd(cumsum(ans_eltlens))
relist(ans_regions, ans_skeleton)

## GRangesList of length 3:
## [[1]] 
## GRanges with 1 range and 0 metadata columns:
##       seqnames    ranges strand
##             
##   [1]        A    [1, 4]      *
## 
## [[2]] 
## GRanges with 1 range and 0 metadata columns:
##       seqnames   ranges strand
##   [1]        A [11, 15]      *
## 
## [[3]] 
## GRanges with 0 ranges and 0 metadata columns:
##      seqnames ranges strand
## 
## ---
## seqlengths:
##   A
##  NA




Thanks to Michael and Hervé for the work on the `Hits` class.


## Self-matching

Recently I had a question about splitting on groups/factors. I
thought the solution Hervé came up with was quite elegant and worth 
sharing.

The problem boils down to converting the 'group1' vector to 'group2':

group1 <- c(3, 1, 2, 2)
group2 <- c(1, 2, 3, 3)




The motivating use case was to split a GRanges based on the 
groups in `group1` but not the implied order.

When we split this GRanges

x <- GRanges(c(rep("chr2", 2), rep("chr19", 2)), IRanges(1:4, width=1))
x

## GRanges with 4 ranges and 0 metadata columns:
##       seqnames    ranges strand
##             
##   [1]     chr2    [1, 1]      *
##   [2]     chr2    [2, 2]      *
##   [3]    chr19    [3, 3]      *
##   [4]    chr19    [4, 4]      *
##   ---
##   seqlengths:
##    chr19  chr2
##       NA    NA




by `group1` the object is re-ordered to match the priority
in `group1`.

split(x, group1)

## GRangesList of length 3:
## $1 
## GRanges with 1 range and 0 metadata columns:
##       seqnames    ranges strand
##             
##   [1]     chr2    [2, 2]      *
## 
## $2 
## GRanges with 2 ranges and 0 metadata columns:
##       seqnames ranges strand
##   [1]    chr19 [3, 3]      *
##   [2]    chr19 [4, 4]      *
## 
## $3 
## GRanges with 1 range and 0 metadata columns:
##       seqnames ranges strand
##   [1]     chr2 [1, 1]      *
## 
## ---
## seqlengths:
##  chr19  chr2
##     NA    NA




This is expected behavior but wasn't exactly the result I was
after. I wanted to respect the grouping in `group1` but not the order.
Essentially I wanted to turn `group1` in to `group2` so
I would get this result:

split(x, group2)

## GRangesList of length 3:
## $1 
## GRanges with 1 range and 0 metadata columns:
##       seqnames    ranges strand
##             
##   [1]     chr2    [1, 1]      *
## 
## $2 
## GRanges with 1 range and 0 metadata columns:
##       seqnames ranges strand
##   [1]     chr2 [2, 2]      *
## 
## $3 
## GRanges with 2 ranges and 0 metadata columns:
##       seqnames ranges strand
##   [1]    chr19 [3, 3]      *
##   [2]    chr19 [4, 4]      *
## 
## ---
## seqlengths:
##  chr19  chr2
##     NA    NA




Solution:

One straight-forward solution is to re-order the GRangesList object 
once created. 

grl <- split(x, group1)
grl[as.character(unique(group1))] 




However the data were large so I was hoping to solve the ordering 
problem before creating the list. The second solution was to
transform `group1` with a "self-match". 

match(group1, group1)

## [1] 1 2 3 3




A "self match" on any vector-like object 'x' is an easy and
very efficient way to assign unique ids to unique values in 'x'.
Also the mapping from unique values in 'x' and ids is guaranteed
to be strictly increasing, which is what I needed for my use
case.

xx <- c("Z", "B", "C", "C", "A", "B")
ids <- match(xx, xx)
ids

## [1] 1 2 3 3 5 2

ids[!duplicated(xx)]  # always strictly increasing

## [1] 1 2 3 5




Hervé's 2 line solution to my splitting problem:

grl2 <- split(x, match(group1, group1))
names(grl2) <- unique(group1)




More detail on the `Hits` class can be found on the man page:

    ?Hits

The ["HOWTOs" vignette](http://www.bioconductor.org/packages/devel/bioc/vignettes/GenomicRanges/inst/doc/GenomicRangesHOWTOs.pdf)
in `GenomicRanges` demonstrates other
practical workflows that make use of the extensive `GenomicRanges`
infrastructure.

    library(GenomicRanges)
    browseVignettes("GenomicRanges")
 

## Courses / Conferences

In January Martin, Sean and Hervé headed south to Ribeirão Preto, Brazil
to teach at the 
[Summer Bioinformatics course](http://bioinformatics.fmrp.usp.br/curso2014/) 
hosted by the Ribeirão Preto Medical School. The course was organized by 
Benilton Carvalho and Houtan Noushmehr and focused on the analysis and 
comprehension of high throughput sequence data. 

In Seattle we offer several short courses throughout the year. The last
was in February given by Martin and Sonali. This was a 2 day introductory 
course covering the basics of working with large data, class structure, 
integrated workflows and the analysis of RNAseq differential expresssion data.

Upcoming events and past course materials can be found on the 
[web site](http://www.bioconductor.org/help/course-materials/).


## Looking ahead

The spring release is just weeks away. We plan to branch on April 11th
and release on the 14th. The full schedule is
[posted online](http://www.bioconductor.org/developers/release-schedule/).

BioC 2014 is in Boston this year (July 30 - August 1). Watch the 
[web site](http://www.bioconductor.org/help/course-materials/2014/BioC2014/)
for updates.

The July newsletter will (probably) include copy number packages on the
rise, the `GenomicFileViews` package, and a summary of tally/coverage 
functions. If there are topics you'd like to see covered let us know!