Overall workflow

1. Experimental design
   • Keep it simple!
   • Replication!
   • Avoid or track batch effects
2. Wet-lab preparation
3. High-throughput sequencing
   • Output: FASTQ files of reads and their quality scores
4. Alignment
   • Many different aligners, some specialized for different purposes
   • Output: BAM files of aligned reads
5. Summary
   • e.g., count of reads overlapping regions of interest (e.g., genes)
6. Statistical analysis
7. Comprehension

Notes

• de novo Assembly outside the scope of Bioconductor

Two simple shiny apps

• BAMSpector – display gene models and underlying support across BAM (aligned read) files

  app <- system.file(package="useR2015", "BAMSpector")
  shiny::runApp(app)

• MAPlotExplorer – summarize two-group differential expression, including drill-down of individual genes. Based on CSAMA 2015 lab by Andrzej Oles.

  app <- system.file(package="useR2015", "MAPlotExplorer")
  shiny::runApp(app)

How Bioconductor helps

Annotation

• Gene (e.g., exons-within-transcipts) models
• Identifier mapping (e.g., ‘HNRPC’ to Entrez identifier used in the gene model)

Standard (large) file input & manipulation, e.g., BAM files of aligned reads

Statistical analysis of differential expression

Gene models – TxDb, GRanges, and GRangesList

library("TxDb.Hsapiens.UCSC.hg19.knownGene")
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
txdb

## TxDb object:
## # Db type: TxDb
## # Supporting package: GenomicFeatures
## # Data source: UCSC
## # Genome: hg19
## # Organism: Homo sapiens
## # UCSC Table: knownGene
## # Resource URL: http://genome.ucsc.edu/
## # Type of Gene ID: Entrez Gene ID
## # Full dataset: yes
## # miRBase build ID: GRCh37
## # transcript_nrow: 82960
## # exon_nrow: 289969
## # cds_nrow: 237533
## # Db created by: GenomicFeatures package from Bioconductor
## # Creation time: 2015-03-19 13:55:51 -0700 (Thu, 19 Mar 2015)
## # GenomicFeatures version at creation time: 1.19.32
## # RSQLite version at creation time: 1.0.0
## # DBSCHEMAVERSION: 1.1

methods(class=class(txdb))

##  [1] $                      $<-                    annotatedDataFrameFrom as.list               
##  [5] asBED                  asGFF                  assayData              assayData<-           
##  [9] cds                    cdsBy                  cdsByOverlaps          coerce                
## [13] columns                combine                contents               dbconn                
## [17] dbfile                 dbInfo                 dbmeta                 dbschema              
## [21] disjointExons          distance               exons                  exonsBy               
## [25] exonsByOverlaps        ExpressionSet          extractUpstreamSeqs    featureNames          
## [29] featureNames<-         fiveUTRsByTranscript   genes                  initialize            
## [33] intronsByTranscript    isActiveSeq            isActiveSeq<-          isNA                  
## [37] keys                   keytypes               mapIds                 mappedkeys            
## [41] mapToTranscripts       metadata               microRNAs              nhit                  
## [45] organism               promoters              revmap                 sample                
## [49] sampleNames            sampleNames<-          saveDb                 select                
## [53] seqinfo                seqinfo<-              seqlevels0             show                  
## [57] species                storageMode            storageMode<-          threeUTRsByTranscript 
## [61] transcripts            transcriptsBy          transcriptsByOverlaps  tRNAs                 
## [65] updateObject          
## see '?methods' for accessing help and source code

exons(txdb)

## GRanges object with 289969 ranges and 1 metadata column:
##            seqnames               ranges strand   |   exon_id
##               <Rle>            <IRanges>  <Rle>   | <integer>
##        [1]     chr1       [11874, 12227]      +   |         1
##        [2]     chr1       [12595, 12721]      +   |         2
##        [3]     chr1       [12613, 12721]      +   |         3
##        [4]     chr1       [12646, 12697]      +   |         4
##        [5]     chr1       [13221, 14409]      +   |         5
##        ...      ...                  ...    ... ...       ...
##   [289965]     chrY [27607404, 27607432]      -   |    277746
##   [289966]     chrY [27635919, 27635954]      -   |    277747
##   [289967]     chrY [59358329, 59359508]      -   |    277748
##   [289968]     chrY [59360007, 59360115]      -   |    277749
##   [289969]     chrY [59360501, 59360854]      -   |    277750
##   -------
##   seqinfo: 93 sequences (1 circular) from hg19 genome

exonsBy(txdb, "tx")

## GRangesList object of length 82960:
## $1 
## GRanges object with 3 ranges and 3 metadata columns:
##       seqnames         ranges strand |   exon_id   exon_name exon_rank
##          <Rle>      <IRanges>  <Rle> | <integer> <character> <integer>
##   [1]     chr1 [11874, 12227]      + |         1        <NA>         1
##   [2]     chr1 [12613, 12721]      + |         3        <NA>         2
##   [3]     chr1 [13221, 14409]      + |         5        <NA>         3
## 
## $2 
## GRanges object with 3 ranges and 3 metadata columns:
##       seqnames         ranges strand | exon_id exon_name exon_rank
##   [1]     chr1 [11874, 12227]      + |       1      <NA>         1
##   [2]     chr1 [12595, 12721]      + |       2      <NA>         2
##   [3]     chr1 [13403, 14409]      + |       6      <NA>         3
## 
## $3 
## GRanges object with 3 ranges and 3 metadata columns:
##       seqnames         ranges strand | exon_id exon_name exon_rank
##   [1]     chr1 [11874, 12227]      + |       1      <NA>         1
##   [2]     chr1 [12646, 12697]      + |       4      <NA>         2
##   [3]     chr1 [13221, 14409]      + |       5      <NA>         3
## 
## ...
## <82957 more elements>
## -------
## seqinfo: 93 sequences (1 circular) from hg19 genome

Identifier mapping – OrgDb

library(org.Hs.eg.db)
org.Hs.eg.db

## OrgDb object:
## | DBSCHEMAVERSION: 2.1
## | Db type: OrgDb
## | Supporting package: AnnotationDbi
## | DBSCHEMA: HUMAN_DB
## | ORGANISM: Homo sapiens
## | SPECIES: Human
## | EGSOURCEDATE: 2015-Mar17
## | EGSOURCENAME: Entrez Gene
## | EGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
## | CENTRALID: EG
## | TAXID: 9606
## | GOSOURCENAME: Gene Ontology
## | GOSOURCEURL: ftp://ftp.geneontology.org/pub/go/godatabase/archive/latest-lite/
## | GOSOURCEDATE: 20150314
## | GOEGSOURCEDATE: 2015-Mar17
## | GOEGSOURCENAME: Entrez Gene
## | GOEGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
## | KEGGSOURCENAME: KEGG GENOME
## | KEGGSOURCEURL: ftp://ftp.genome.jp/pub/kegg/genomes
## | KEGGSOURCEDATE: 2011-Mar15
## | GPSOURCENAME: UCSC Genome Bioinformatics (Homo sapiens)
## | GPSOURCEURL: ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19
## | GPSOURCEDATE: 2010-Mar22
## | ENSOURCEDATE: 2015-Mar13
## | ENSOURCENAME: Ensembl
## | ENSOURCEURL: ftp://ftp.ensembl.org/pub/current_fasta
## | UPSOURCENAME: Uniprot
## | UPSOURCEURL: http://www.UniProt.org/
## | UPSOURCEDATE: Tue Mar 17 18:48:15 2015

## 
## Please see: help('select') for usage information

OrgDb objects

• Curated resources, underlying sqlite data base, like TxDb
• make your own: AnnotationForge (but see the AnnotationHub, below!)
• ‘select’ interface: keytypes(), columns(), keys(), select(), mapIds()

select()

• Vector of keys, desired columns

• Specification of key type

  select(org.Hs.eg.db, c("BRCA1", "PTEN"), c("ENTREZID", "GENENAME"), "SYMBOL")

  ##   SYMBOL ENTREZID                       GENENAME
  ## 1  BRCA1      672   breast cancer 1, early onset
  ## 2   PTEN     5728 phosphatase and tensin homolog

  keytypes(org.Hs.eg.db)

  ##  [1] "ENTREZID"     "PFAM"         "IPI"          "PROSITE"      "ACCNUM"       "ALIAS"       
  ##  [7] "ENZYME"       "MAP"          "PATH"         "PMID"         "REFSEQ"       "SYMBOL"      
  ## [13] "UNIGENE"      "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS" "GENENAME"     "UNIPROT"     
  ## [19] "GO"           "EVIDENCE"     "ONTOLOGY"     "GOALL"        "EVIDENCEALL"  "ONTOLOGYALL" 
  ## [25] "OMIM"         "UCSCKG"

  columns(org.Hs.eg.db)

  ##  [1] "ENTREZID"     "PFAM"         "IPI"          "PROSITE"      "ACCNUM"       "ALIAS"       
  ##  [7] "CHR"          "CHRLOC"       "CHRLOCEND"    "ENZYME"       "MAP"          "PATH"        
  ## [13] "PMID"         "REFSEQ"       "SYMBOL"       "UNIGENE"      "ENSEMBL"      "ENSEMBLPROT" 
  ## [19] "ENSEMBLTRANS" "GENENAME"     "UNIPROT"      "GO"           "EVIDENCE"     "ONTOLOGY"    
  ## [25] "GOALL"        "EVIDENCEALL"  "ONTOLOGYALL"  "OMIM"         "UCSCKG"

Related functionality

• mapIds() – special case for mapping from 1 identifier to another

• OrganismDb objects: combined org.*, TxDb.*, and other annotation resources for easy access

  library(Homo.sapiens)
  select(Homo.sapiens, c("BRCA1", "PTEN"), 
         c("TXNAME", "TXCHROM", "TXSTART", "TXEND"), 
         "SYMBOL")

  ##    SYMBOL     TXNAME TXCHROM  TXSTART    TXEND
  ## 1   BRCA1 uc010whl.2   chr17 41196312 41276132
  ## 2   BRCA1 uc002icp.4   chr17 41196312 41277340
  ## 3   BRCA1 uc010whm.2   chr17 41196312 41277340
  ## 4   BRCA1 uc002icu.3   chr17 41196312 41277468
  ## 5   BRCA1 uc010cyx.3   chr17 41196312 41277468
  ## 6   BRCA1 uc002icq.3   chr17 41196312 41277500
  ## 7   BRCA1 uc002ict.3   chr17 41196312 41277500
  ## 8   BRCA1 uc010whn.2   chr17 41196312 41277500
  ## 9   BRCA1 uc010who.3   chr17 41196312 41277500
  ## 10  BRCA1 uc010whp.2   chr17 41196312 41322420
  ## 11  BRCA1 uc010whq.1   chr17 41215350 41256973
  ## 12  BRCA1 uc002idc.1   chr17 41215350 41277468
  ## 13  BRCA1 uc010whr.1   chr17 41215350 41277468
  ## 14  BRCA1 uc002idd.3   chr17 41243117 41276132
  ## 15  BRCA1 uc002ide.1   chr17 41243452 41256973
  ## 16  BRCA1 uc010cyy.1   chr17 41243452 41277340
  ## 17  BRCA1 uc010whs.1   chr17 41243452 41277468
  ## 18  BRCA1 uc010cyz.2   chr17 41243452 41277500
  ## 19  BRCA1 uc010cza.2   chr17 41243452 41277500
  ## 20  BRCA1 uc010wht.1   chr17 41243452 41277500
  ## 21   PTEN uc001kfb.3   chr10 89623195 89728532
  ## 22   PTEN uc021pvw.1   chr10 89623195 89728532

Other annotation resources – biomaRt, AnnotationHub

## NEEDS INTERNET ACCESS !!
library(biomaRt)
head(listMarts(), 3)                      ## list marts
head(listDatasets(useMart("ensembl")), 3) ## mart datasets
ensembl <-                                ## fully specified mart
    useMart("ensembl", dataset = "hsapiens_gene_ensembl")

head(listFilters(ensembl), 3)             ## filters
myFilter <- "chromosome_name"
substr(filterOptions(myFilter, ensembl), 1, 50) ## return values
myValues <- c("21", "22")
head(listAttributes(ensembl), 3)          ## attributes
myAttributes <- c("ensembl_gene_id","chromosome_name")

## assemble and query the mart
res <- getBM(attributes =  myAttributes, filters =  myFilter,
             values =  myValues, mart = ensembl)

library(AnnotationHub)
hub <- AnnotationHub()
hub
query(hub, c("Ensembl", "80", "gtf"))
## ensgtf = display(hub)                   # visual choice
hub["AH47107"]
gtf <- hub[["AH47107"]]
gtf
txdb <- GenomicFeatures::makeTxDbFromGRanges(gtf)

library(AnnotationHub)
hub <- AnnotationHub()
query(hub, "OrgDb")

BAM files of aligned reads – GenomicAlignments

Recall: overall workflow

1. Experimental design
2. Wet-lab preparation
3. High-throughput sequencing
4. Alignment
   • Whole genome, vs. transcriptome
5. Summary
6. Statistical analysis
7. Comprehension

BAM files of aligned reads

• Header

  @HD     VN:1.0  SO:coordinate
  @SQ     SN:chr1 LN:249250621
  @SQ     SN:chr10        LN:135534747
  @SQ     SN:chr11        LN:135006516
  ...
  @SQ     SN:chrY LN:59373566
  @PG     ID:TopHat       VN:2.0.8b       CL:/home/hpages/tophat-2.0.8b.Linux_x86_64/tophat --mate-inner-dist 150 --solexa-quals --max-multihits 5 --no-discordant --no-mixed --coverage-search --microexon-search --library-type fr-unstranded --num-threads 2 --output-dir tophat2_out/ERR127306 /home/hpages/bowtie2-2.1.0/indexes/hg19 fastq/ERR127306_1.fastq fastq/ERR127306_2.fastq

• Alignments

  • ID, flag, alignment and mate

    ERR127306.7941162       403     chr14   19653689        3       72M             =       19652348        -1413  ...
    ERR127306.22648137      145     chr14   19653692        1       72M             =       19650044        -3720  ...

  • Sequence and quality

    ... GAATTGATCAGTCTCATCTGAGAGTAACTTTGTACCCATCACTGATTCCTTCTGAGACTGCCTCCACTTCCC        *'%%%%%#&&%''#'&%%%)&&%%$%%'%%'&*****$))$)'')'%)))&)%%%%$'%%%%&"))'')%))
    ... TTGATCAGTCTCATCTGAGAGTAACTTTGTACCCATCACTGATTCCTTCTGAGACTGCCTCCACTTCCCCAG        '**)****)*'*&*********('&)****&***(**')))())%)))&)))*')&***********)****

  • Tags

    ... AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:72 YT:Z:UU NH:i:2  CC:Z:chr22      CP:i:16189276   HI:i:0
    ... AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:72 YT:Z:UU NH:i:3  CC:Z:=  CP:i:19921600   HI:i:0

• Typically, sorted (by position) and indexed (‘.bai’ files)

GenomicAlignments

• Use an example BAM file (fl could be the path to your own BAM file)

  ## example BAM data
  library(RNAseqData.HNRNPC.bam.chr14)
  ## one BAM file
  fl <- RNAseqData.HNRNPC.bam.chr14_BAMFILES[1]
  ## Let R know that this is a BAM file, not just a character vector
  library(Rsamtools)
  bfl <- BamFile(fl)

• Input the data into R

  aln <- readGAlignments(bfl)
  aln

  ## GAlignments object with 800484 alignments and 0 metadata columns:
  ##            seqnames strand       cigar    qwidth     start       end     width     njunc
  ##               <Rle>  <Rle> <character> <integer> <integer> <integer> <integer> <integer>
  ##        [1]    chr14      +         72M        72  19069583  19069654        72         0
  ##        [2]    chr14      +         72M        72  19363738  19363809        72         0
  ##        [3]    chr14      -         72M        72  19363755  19363826        72         0
  ##        [4]    chr14      +         72M        72  19369799  19369870        72         0
  ##        [5]    chr14      -         72M        72  19369828  19369899        72         0
  ##        ...      ...    ...         ...       ...       ...       ...       ...       ...
  ##   [800480]    chr14      -         72M        72 106989780 106989851        72         0
  ##   [800481]    chr14      +         72M        72 106994763 106994834        72         0
  ##   [800482]    chr14      -         72M        72 106994819 106994890        72         0
  ##   [800483]    chr14      +         72M        72 107003080 107003151        72         0
  ##   [800484]    chr14      -         72M        72 107003171 107003242        72         0
  ##   -------
  ##   seqinfo: 93 sequences from an unspecified genome

  • readGAlignmentPairs() / readGAlignmentsList() if paired-end data
  • Lots of things to do, including all the GRanges / GRangesList operations

  methods(class=class(aln))

  ##   [1] !=                     [                      [<-                    %in%                  
  ##   [5] <                      <=                     ==                     >                     
  ##   [9] >=                     aggregate              anyNA                  append                
  ##  [13] as.character           as.complex             as.data.frame          as.env                
  ##  [17] as.integer             as.list                as.logical             as.numeric            
  ##  [21] as.raw                 c                      cigar                  coerce                
  ##  [25] compare                countOverlaps          coverage               duplicated            
  ##  [29] elementMetadata        elementMetadata<-      end                    eval                  
  ##  [33] export                 extractROWS            findCompatibleOverlaps findOverlaps          
  ##  [37] findSpliceOverlaps     granges                grglist                head                  
  ##  [41] high2low               junctions              length                 mapCoords             
  ##  [45] mapFromAlignments      mapToAlignments        match                  mcols                 
  ##  [49] mcols<-                metadata               metadata<-             mstack                
  ##  [53] names                  names<-                narrow                 njunc                 
  ##  [57] NROW                   overlapsAny            parallelSlotNames      pintersect            
  ##  [61] pmapCoords             pmapFromAlignments     pmapToAlignments       qnarrow               
  ##  [65] qwidth                 ranges                 rank                   relist                
  ##  [69] relistToClass          rename                 rep                    rep.int               
  ##  [73] replaceROWS            rev                    rglist                 rname                 
  ##  [77] rname<-                ROWNAMES               seqinfo                seqinfo<-             
  ##  [81] seqlevelsInUse         seqnames               seqnames<-             shiftApply            
  ##  [85] show                   showAsCell             sort                   split                 
  ##  [89] split<-                start                  strand                 strand<-              
  ##  [93] subset                 subsetByOverlaps       summarizeOverlaps      table                 
  ##  [97] tail                   tapply                 unique                 update                
  ## [101] updateObject           values                 values<-               width                 
  ## [105] window                 window<-               with                   xtfrm                 
  ## see '?methods' for accessing help and source code

• Caveat emptor: BAM files are large. Normally you will restrict the input to particular genomic ranges, or iterate through the BAM file. Key Bioconductor functions (e.g., GenomicAlignments::summarizeOverlaps() do this data management step for you. See next section!

Other formats and packages

Large data – BiocParallel, GenomicFiles

Challenges & solutions

Starting point

• Matrix of counts of reads overlapping each region of interest
• Counts provide statistical information – larger counts indicate greater confidence that the read was observed. Standardized measures (e.g., RPKM) ignore this information and therefore lose statistical power.

Normalization

• Differences in read counts per sample for purely technical reasons
• Simple scaling by total read count inadequate – induces correlations with low-count reads
• General solution: scale by a more robust measure of size, e.g., log geometric mean, quantile, …

Error model

• Poisson ‘shot’ noise of reads sampled from a genome. E.g., longer genes receive more aligned reads compared to shorter genes with identical expression.
• Additional biological variation due to differences between genes and individuals
• Common modeling assumptions: negative binomial variance
  • Dispersion parameter

Limited sample size

• A handful of samples in each treatment
• Many 1000’s of statistical tests
• Challenge – limited statistical power
• Solution – borrow information
  • Estimate variance as weighted average of per gene variance, and average variance of all genes
  • Per-gene variances are estimated precisely, though with some loss of accuracy
  • Example of moderated test statistic

Multiple testing

• Need to adjust for multiple comparisons
• Reducing number of tests enhances statistical power

• Filter genes to exclude from testing using a priori criteria

  • Not biologically interesting
  • Not statistically interesting under the null, e.g., insufficient counts across samples

Work flow

library("airway")
data(airway)
airway

## class: SummarizedExperiment 
## dim: 64102 8 
## exptData(1): ''
## assays(1): counts
## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
## rowRanges metadata column names(0):
## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
## colData names(9): SampleName cell ... Sample BioSample

## main components of SummarizedExperiment
head(assay(airway))

##                 SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520
## ENSG00000000003        679        448        873        408       1138       1047        770
## ENSG00000000005          0          0          0          0          0          0          0
## ENSG00000000419        467        515        621        365        587        799        417
## ENSG00000000457        260        211        263        164        245        331        233
## ENSG00000000460         60         55         40         35         78         63         76
## ENSG00000000938          0          0          2          0          1          0          0
##                 SRR1039521
## ENSG00000000003        572
## ENSG00000000005          0
## ENSG00000000419        508
## ENSG00000000457        229
## ENSG00000000460         60
## ENSG00000000938          0

colData(airway)

## DataFrame with 8 rows and 9 columns
##            SampleName     cell      dex    albut        Run avgLength Experiment    Sample
##              <factor> <factor> <factor> <factor>   <factor> <integer>   <factor>  <factor>
## SRR1039508 GSM1275862   N61311    untrt    untrt SRR1039508       126  SRX384345 SRS508568
## SRR1039509 GSM1275863   N61311      trt    untrt SRR1039509       126  SRX384346 SRS508567
## SRR1039512 GSM1275866  N052611    untrt    untrt SRR1039512       126  SRX384349 SRS508571
## SRR1039513 GSM1275867  N052611      trt    untrt SRR1039513        87  SRX384350 SRS508572
## SRR1039516 GSM1275870  N080611    untrt    untrt SRR1039516       120  SRX384353 SRS508575
## SRR1039517 GSM1275871  N080611      trt    untrt SRR1039517       126  SRX384354 SRS508576
## SRR1039520 GSM1275874  N061011    untrt    untrt SRR1039520       101  SRX384357 SRS508579
## SRR1039521 GSM1275875  N061011      trt    untrt SRR1039521        98  SRX384358 SRS508580
##               BioSample
##                <factor>
## SRR1039508 SAMN02422669
## SRR1039509 SAMN02422675
## SRR1039512 SAMN02422678
## SRR1039513 SAMN02422670
## SRR1039516 SAMN02422682
## SRR1039517 SAMN02422673
## SRR1039520 SAMN02422683
## SRR1039521 SAMN02422677

rowRanges(airway)

## GRangesList object of length 64102:
## $ENSG00000000003 
## GRanges object with 17 ranges and 2 metadata columns:
##        seqnames               ranges strand   |   exon_id       exon_name
##           <Rle>            <IRanges>  <Rle>   | <integer>     <character>
##    [1]        X [99883667, 99884983]      -   |    667145 ENSE00001459322
##    [2]        X [99885756, 99885863]      -   |    667146 ENSE00000868868
##    [3]        X [99887482, 99887565]      -   |    667147 ENSE00000401072
##    [4]        X [99887538, 99887565]      -   |    667148 ENSE00001849132
##    [5]        X [99888402, 99888536]      -   |    667149 ENSE00003554016
##    ...      ...                  ...    ... ...       ...             ...
##   [13]        X [99890555, 99890743]      -   |    667156 ENSE00003512331
##   [14]        X [99891188, 99891686]      -   |    667158 ENSE00001886883
##   [15]        X [99891605, 99891803]      -   |    667159 ENSE00001855382
##   [16]        X [99891790, 99892101]      -   |    667160 ENSE00001863395
##   [17]        X [99894942, 99894988]      -   |    667161 ENSE00001828996
## 
## ...
## <64101 more elements>
## -------
## seqinfo: 722 sequences (1 circular) from an unspecified genome

## e.g., coordinated subset to include dex 'trt'  samples
airway[, airway$dex == "trt"]

## class: SummarizedExperiment 
## dim: 64102 4 
## exptData(1): ''
## assays(1): counts
## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
## rowRanges metadata column names(0):
## colnames(4): SRR1039509 SRR1039513 SRR1039517 SRR1039521
## colData names(9): SampleName cell ... Sample BioSample

## e.g., keep only rows with non-zero counts
airway <- airway[rowSums(assay(airway)) != 0, ]

Key references

• Irizarry R, et al. (2015) Biomedical Data Science. Course Notes, EdX PH525.1x.
• Huber W, et al. (2015) Orchestrating high-throughput genomic analysis with Bioconductor. Nature Methods 12:115-121; doi:10.1038/nmeth.3252 (full-text free with registration).
• Lawrence M, Huber W, Pag&egraves;s H, Aboyoun P, Carlson M, et al. (2013) Software for Computing and Annotating Genomic Ranges. PLoS Comput Biol 9(8): e1003118. doi: 10.1371/journal.pcbi.1003118

sessionInfo()

sessionInfo()

## R version 3.2.1 Patched (2015-06-19 r68553)
## Platform: x86_64-unknown-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.2 LTS
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
##  [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
## [10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] shiny_0.12.0                            ggplot2_1.0.1                          
##  [3] airway_0.102.0                          RNAseqData.HNRNPC.bam.chr14_0.6.0      
##  [5] Homo.sapiens_1.1.2                      TxDb.Hsapiens.UCSC.hg19.knownGene_3.1.2
##  [7] org.Hs.eg.db_3.1.2                      GO.db_3.1.2                            
##  [9] RSQLite_1.0.0                           DBI_0.3.1                              
## [11] OrganismDbi_1.10.0                      GenomicFeatures_1.20.1                 
## [13] AnnotationDbi_1.30.1                    Biobase_2.28.0                         
## [15] GenomicFiles_1.4.0                      BiocParallel_1.2.2                     
## [17] rtracklayer_1.28.4                      GenomicAlignments_1.4.1                
## [19] Rsamtools_1.20.4                        DESeq2_1.8.1                           
## [21] RcppArmadillo_0.5.200.1.0               Rcpp_0.11.6                            
## [23] GenomicRanges_1.20.5                    GenomeInfoDb_1.4.0                     
## [25] Biostrings_2.36.1                       XVector_0.8.0                          
## [27] IRanges_2.2.4                           S4Vectors_0.6.0                        
## [29] BiocGenerics_0.14.0                     AnnotationHub_2.0.2                    
## [31] BiocStyle_1.6.0                         BiocInstaller_1.18.3                   
## 
## loaded via a namespace (and not attached):
##  [1] httr_0.6.1                   splines_3.2.1                Formula_1.2-1               
##  [4] interactiveDisplayBase_1.6.0 latticeExtra_0.6-26          RBGL_1.44.0                 
##  [7] yaml_2.1.13                  lattice_0.20-31              digest_0.6.8                
## [10] RColorBrewer_1.1-2           colorspace_1.2-6             htmltools_0.2.6             
## [13] httpuv_1.3.2                 plyr_1.8.2                   XML_3.98-1.2                
## [16] biomaRt_2.24.0               genefilter_1.50.0            zlibbioc_1.14.0             
## [19] xtable_1.7-4                 snow_0.3-13                  scales_0.2.4                
## [22] annotate_1.46.0              nnet_7.3-9                   proto_0.3-10                
## [25] survival_2.38-2              magrittr_1.5                 mime_0.3                    
## [28] evaluate_0.7                 MASS_7.3-41                  foreign_0.8-63              
## [31] graph_1.46.0                 tools_3.2.1                  formatR_1.2                 
## [34] stringr_1.0.0                munsell_0.4.2                locfit_1.5-9.1              
## [37] cluster_2.0.2                lambda.r_1.1.7               futile.logger_1.4.1         
## [40] grid_3.2.1                   RCurl_1.95-4.6               labeling_0.3                
## [43] bitops_1.0-6                 rmarkdown_0.6.1              codetools_0.2-11            
## [46] gtable_0.1.2                 reshape2_1.4.1               R6_2.0.1                    
## [49] gridExtra_0.9.1              knitr_1.10.5                 Hmisc_3.16-0                
## [52] futile.options_1.0.0         stringi_0.4-1                geneplotter_1.46.0          
## [55] rpart_4.1-9                  acepack_1.3-3.3