Hi everyone,

I am having a small problem with Diffbind, after successfully running through the diffbind algorithm and getting a GRanges object with about 1000 differential binding sties for my contrast every time I call one of the plotting functions:

dba.plotMA

dba.plotVenn

I get the following error:

Error in pv.DBAplotMA(DBA, contrast = contrast, method = method, bMA = !bXY,  : 
  Specified contrast number is greater than number of contrasts

It does not matter how I run through the algorithm it seems I am always left with this error at the end. I'm sure the answer is simple, but I am stumped.

Any help would be much appreciated.

Ben

Hi Ben-

Could you include the output of sessionInfo()?

Also some more of your script. Specifically, the contents of the DBA object before you call dba.plotMA(), and the exact call to dba.plotMA(). If your DBA object is called "myDBA", I just need to see the output of:

> myDBA

Cheers-

Rory

Rory, thanks so much for getting back to me. Please forgive my ignorance on this, I am new to using these tools.

Certainly, I use the following calls to generate the dba objects:

ATAC <- dba(sampleSheet="NeurosphereATAC_TSSintersect_samplesheet_mod1.txt", peakFormat = "narrow")
ATAC_count<-dba.count(ATAC)
ATAC_contrast<-dba.contrast(ATAC_count, categories = DBA_TISSUE, minMembers = 2)
ATAC_analyze<-dba.analyze(ATAC_contrast)
ATAC_report<-dba.report(ATAC_analyze, contrast = 3)

When I print "ATAC_analyze" I get this:

6 Samples, 18179 sites in matrix:
          ID  Tissue Replicate Caller Intervals FRiP
1 50000cells     NSC         1 counts     18179 0.06
2   625cells     NSC         2 counts     18179 0.03
3    Muscle1  Muscle         1 counts     18179 0.10
4    Muscle2  Muscle         2 counts     18179 0.08
5     HelpT1 HelperT         1 counts     18179 0.25
6     HelpT2 HelperT         2 counts     18179 0.10

6 Contrasts:
   Group1 Members1   Group2 Members2 DB.edgeR
1     NSC        2   Muscle        2    11580
2     NSC        2  HelperT        2     8960
3     NSC        2     !NSC        4     2265
4  Muscle        2  HelperT        2    10013
5  Muscle        2  !Muscle        4     1580
6 HelperT        2 !HelperT        4     6272

I want to look at the third contrast so I called the third contrast for dba.report. When I print ATAC_report I get this:

GRanges with 2265 ranges and 6 metadata columns:
        seqnames                 ranges strand   |             Conc         Conc_NSC         Conc_.NSC              Fold              p.value                  FDR
           <Rle>              <IRanges>  <Rle>   |        <numeric>        <numeric>         <numeric>         <numeric>            <numeric>            <numeric>
  11245     chr4   [22415513, 22416108]      *   | 7.55136037763646 9.10122857073246  2.75742577311184  6.34380279762062 1.48869725165757e-17  2.7063027337883e-13
   5589    chr15   [44538034, 44538282]      *   | 6.45245996804235 8.00958843826776   1.3277550985802  6.68183333968755 6.62357308376236e-17  6.0204967544858e-13
  10437     chr3   [54719809, 54720072]      *   | 5.99001680676936 7.55623111284682 0.299744477065765  7.25648663578105 7.31286120757036e-16 4.43135012974738e-12
  13795     chr6   [82351610, 82352424]      *   | 6.90291058153769 8.45032465044132  2.20526967837527  6.24505497206605 4.83078008298791e-15 2.19546877821593e-11
  14979     chr7   [70589167, 70589657]      *   | 5.75408953481454 7.30864698940696 0.755564427600979  6.55308256180598 6.47125730053609e-15 2.35281972932891e-11
    ...      ...                    ...    ... ...              ...              ...               ...               ...                  ...                  ...
    214     chr1 [ 57435100,  57435266]      *   | 3.97841979098216  5.0616971238251  2.79588924967147  2.26580787415363   0.0124395641100217   0.0999508661791857
   1092    chr10 [  8605868,   8606347]      *   | 6.89697266153719 8.11047070405685  5.34271563340228  2.76775507065458   0.0124408020917137   0.0999508661791857
  13968     chr6 [ 98978292,  98979030]      *   | 5.57193225775836 6.46860964206819  4.75872023161371  1.70988941045448   0.0124439359707428   0.0999508661791857
    946     chr1 [182260440, 182260606]      *   | 6.58441639964778 4.86484450245853  7.01543016530217 -2.15058566284363   0.0124478112673787   0.0999508661791857
  16565     chr8 [124103011, 124103330]      *   |  7.4107899639179 5.85136405293126  7.82260337134972 -1.97123931841846   0.0124584917830042   0.0999924600985577
  ---
  seqlengths:
    chr4 chr15  chr3  chr6  chr7  chr9  chr1 chr10  chr8 chr11 chr18  chr5 chr14  chr2 chr13  chrX chr17 chr12 chr16 chr19
      NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA

This is the correct contrast, but then when I go to call the plotting functions such as:

dba.plotMA(ATAC_report)

That is when I get the error.

Here is the sessionInfo() call:

R version 3.1.1 (2014-07-10)
Platform: x86_64-apple-darwin13.1.0 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] proto_0.3-10            directlabels_2013.6.15  quadprog_1.5-5          lattice_0.20-29         ggbiplot_0.55           scales_0.2.4            plyr_1.8.1              devtools_1.6           
 [9] cummeRbund_2.6.1        Gviz_1.8.4              rtracklayer_1.24.2      fastcluster_1.1.13      reshape2_1.4            ggplot2_1.0.0           RSQLite_0.11.4          DBI_0.3.1              
[17] RColorBrewer_1.0-5      fields_7.1              maps_2.3-9              spam_1.0-1              pheatmap_0.7.7          DiffBind_1.10.2         GenomicAlignments_1.0.6 BSgenome_1.32.0        
[25] Rsamtools_1.16.1        Biostrings_2.32.1       XVector_0.4.0           limma_3.20.9            GenomicRanges_1.16.4    GenomeInfoDb_1.0.2      IRanges_1.22.10         BiocGenerics_0.10.0    
[33] BiocInstaller_1.14.2   

loaded via a namespace (and not attached):
 [1] acepack_1.3-3.3          amap_0.8-12              AnnotationDbi_1.26.1     base64enc_0.1-2          BatchJobs_1.4            BBmisc_1.7               Biobase_2.24.0           BiocParallel_0.6.1      
 [9] biomaRt_2.20.0           biovizBase_1.12.3        bitops_1.0-6             brew_1.0-6               caTools_1.17.1           checkmate_1.4            cluster_1.15.3           codetools_0.2-9         
[17] colorspace_1.2-4         dichromat_2.0-0          digest_0.6.4             edgeR_3.6.8              fail_1.2                 foreach_1.4.2            foreign_0.8-61           Formula_1.1-2           
[25] gdata_2.13.3             GenomicFeatures_1.16.3   gplots_2.14.2            gtable_0.1.2             gtools_3.4.1             Hmisc_3.14-5             iterators_1.0.7          KernSmooth_2.23-13      
[33] labeling_0.3             latticeExtra_0.6-26      MASS_7.3-35              matrixStats_0.10.0       munsell_0.4.2            nnet_7.3-8               R.methodsS3_1.6.1        Rcpp_0.11.3             
[41] RCurl_1.95-4.3           rpart_4.1-8              sendmailR_1.2-1          splines_3.1.1            stats4_3.1.1             stringr_0.6.2            survival_2.37-7          tools_3.1.1             
[49] VariantAnnotation_1.10.5 XML_3.98-1.1             zlibbioc_1.10.0    

Thanks so much for your help!

Hi Ben, 

You should pass the DBA object, not the report, to dba.plotMA (and all the other plotting functions as well):

>  dba.plotMA(ATAC_analyze, contrast=3)

Cheers-

Rory

Ah! I see, thanks so much for your help!