Registration and Call for Abstracts Open for Bioc2024

Menu

Newsletters

2016

2015

2014

 

Bioconductor Newsletter

posted by Valerie Obenchain, April 2014

Contents

Introduction

Welcome to the first Bioconductor newsletter! In this space we will provide project updates from the Seattle team and recent contributions from core developers. Topics will include infrastructure reorganizations, new packages and features, conferences, and other areas of interest.

The newsletter will be on a quarterly schedule with the intent of summarizing major developments of the past 3 months. Comments and suggestions are welcome, please send to Valerie at vobencha@fhcrc.org.

Infrastructure

GenomicAlignments package

Last fall Hervé did some housecleaning on the infrastructure packages. All GAlignments* classes and related functions were moved from GenomicRanges and Rsamtools to the new GenomicAlignments package. Classes and functions from GenomicRanges that built on functionality in Rsamtools were also moved. The reorg substantially cleaned up our dependency tree and will be reflected in the upcoming release, BioC 2.14. (ASCII art courtesy of Hervé).

                  IRanges
                  ^     ^
                  |     |
                  |    XVector
                  |     ^
                  |     |
       GenomicRanges   Biostrings
                  ^     ^
                  |     |
                 Rsamtools
                     ^
                     |
                 BSgenome
                     ^
                     |
             GenomicAlignments
                     ^
                     |
                 rtracklayer
                     ^
                     |
               GenomicFeatures

Junction tools in GenomicAlignments

New tools for extracting junctions from genomic alignments have been added to the GenomicAlignments package. The junctions() function, previously known as gaps(), extracts the portion of a read corresponding to the N operation in the CIGAR.

New junction-related operations:

  • summarize the number of reads that overlap each junction
  • retrieve junction motifs based on dinucleotides at junction boundaries (BSgenome package required)
  • import junction files generated by TopHat and STAR aligners
  • determine if the RNA-seq protocol was stranded

GenomeInfoDb package

Converting between chromosome names from different institutions presents a challenge when working with ranges and annotations. The new GenomeInfoDb package is a step towards resolving these inconsistencies by providing a mapping interface between styles (e.g., UCSC, NCBI, Ensembl, …) and organisms.

Low-level functions related to seqnames and SeqnamesStyle were moved from GenomicRanges and AnnotationDbi into GenomeInfoDb. Some functions were deprecated, some renamed. See the NEWS file for a complete listing of activities.

The package is intended to serve as the base for methods querying or transforming seqlevels. Future work will address the problem of auto-detecting incompatible styles between objects in operations that require matching seqlevels.

Sonali, Marc, Hervé and Martin have all contributed to this work.

Currently 9 organisms are supported:

library(GenomeInfoDb)
names(genomeStyles())

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

There are methods for style discovery,

seqlevelsStyle("chr3")

## [1] "UCSC"

and for converting between styles. Here we map from UCSC to NCBI:

mapSeqlevels(c("chrII", "chrIII", "chrM"), "NCBI")

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

Networks / Annotations

AnnotationHub package

The AnnotationHub resource continues to grow steadily. As of March 2014 the hub supports 204 genomes and 98 species from 7 data providers (Ensembl, NCBI, UCSC, dbSNP, HAEMCODE, Encode, Inparanoid 8). The hub is a repository of annotation data parsed from original file format to R objects for convient computing. For example, BED, GTF, and track files from UCSC are available as GRanges. VCF files from dbSNP have been parsed into VCF objects and FASTA files are available as FaFile objects.

Marc, Dan and Paul have led this effort.

RefNet package

A new package to watch for in the Spring release is RefNet. This package offers access to molecular interaction data in computable form and is the first of its kind in Bioconductor. The package offers access to 2 transcription factor networks and Recon 2 (consensus metabolic reconstruction compiled from 5 resources). Queries are made using the PSICQUIC or RefNet.db packages. (Note: RefNet and RefNet.db will be available in BioC 2.14.)

These data can be used to curate or study experiment-specific pathways, explore the function and relationships of differentially expressed genes, and offer insight to a variety of explorations at the molecular level.

RefNet is extensible and we encourage contributions. Networks, known interactions, infection pathways, disease-host interactions, or pathways from manuscripts would be valuable additions. Guidance for contributing data are included in the package.

Credit for RefNet, RefNet.db and PSICQUIC goes to Paul with contributions from Mentored Projects mentees.

Web Site-ish

biocViews terms

The biocViews terms used to categorize packages have been been refreshed. New terms were added, old terms removed and packages with old or invalid terms were updated. The updated web-interface allows searching by biocViews term or subject. A check box option displays terms not associated with any package. Thanks to Sonali and Dan for the work on this.

BiocCheck package

BiocCheck is a new package that encapsulates Bioconductor package guidelines and best practices. It was designed to provide guidance to new package authors and is run by the SinglePackageBuilder during the package review process.

BiocCheck can also be used from the command line

R CMD BiocCheck package

and should be run after R CMD check successfully completes. A full-featured vignette explains the suggestions made by the package. Dan is the architect behind this one.

Programming bits

Hits class

One of my favorite objects is the Hits class. It is a simple, straightforward class that can be used for complicated indexing and combining. It stores a set of “hits” between two vector-like objects and is most familiar as the output of a findOverlaps operation.

Here we use Hits to identify and remove the “shared” regions in a set of ranges. The shared region in the original set is positions 5 to 10. Range 3 is completely within this shared range.

library(GenomicRanges)
gr1 <- GRanges("A", IRanges(c(1, 5, 5), c(10, 15, 10)))
gr1

## GRanges with 3 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [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 
##    <integer>   <integer> 
##  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 
##    <integer>   <integer> 
##  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
##          <Rle> <IRanges>  <Rle>
##   [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 to 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
##          <Rle> <IRanges>  <Rle>
##   [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
##          <Rle> <IRanges>  <Rle>
##   [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
##          <Rle> <IRanges>  <Rle>
##   [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
##          <Rle> <IRanges>  <Rle>
##   [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 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 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.

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.

BioC 2014 is in Boston this year (July 30 - August 1). Watch the web site 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!