avs.vse {RTN} | R Documentation

## Variant set enrichment (VSE) analysis.

### Description

The VSE method tests the enrichment of an AVS for a particular trait in a genomic annotation.

### Usage

    avs.vse(object, annotation, maxgap=0, pValueCutoff=0.05,
    pAdjustMethod="bonferroni", boxcox=TRUE, lab="annotation",
    glist=NULL, minSize=100, verbose=TRUE)

### Arguments

`object` | an object. When this function is implemented as the S4 method of class `AVS-class`, this argument is an object of class 'AVS'.
`annotation` | a data frame with genomic annotations listing chromosome coordinates to which a particular property or function has been attributed. It should include the following columns: <CHROM>, <START>, <END> and <ID>. The <ID> column can be any genomic identifier, while values in <CHROM> should be listed in ['chr1', 'chr2', 'chr3' ..., 'chrX']. Both <START> and <END> columns correspond to chromosome positions mapped to the human genome assembly used to build the AVS object.
`maxgap` | a single integer value specifying the max distant (kb) between the AVS and the annotation used to compute the enrichment analysis.
`pValueCutoff` | a single numeric value specifying the cutoff for p-values considered significant.
`pAdjustMethod` | a single character value specifying the p-value adjustment method to be used (see 'p.adjust' for details).
`boxcox` | a single logical value specifying to use Box-Cox procedure to find a transformation of the null that approaches normality (when boxcox=TRUE) or not (when boxcox=FALSE). See `powerTransform` and `bcPower`.
`lab` | a single character value specifying a name for the annotation dataset (this option is overrided if 'glist' is used).
`glist` | an optional list with character vectors mapped to the 'annotation' data via <ID> column. This option can be used to run a batch mode for gene sets and regulons.
`minSize` | if 'glist' is provided, this argument is a single integer or numeric value specifying the minimum number of elements for each gene set in the 'glist'. Gene sets with fewer than this number are removed from the analysis.
`verbose` | a single logical value specifying to display detailed messages (when verbose=TRUE) or not (when verbose=FALSE).

### Value

a data frame in the slot "results", see 'what' options in `avs.get`.

### Author(s)

Mauro Castro

### See Also

`AVS-class`

### Examples

    ## Not run:
    # This example requires the RTNdata package! (currently available under request)
    library(RTNdata.LDHapMapRel27.hg18)
    library(Fletcher2013b)
    library(TxDb.Hsapiens.UCSC.hg18.knownGene)

    ##################################################
    ### Build AVS and random AVSs (mapped to hg18)
    ##################################################

    #--- step 1: load 'risk SNPs' data (e.g. BCa risk SNPs from the GWAS catalog)
    data(bcarisk, package = "RTNdata.LDHapMapRel27.hg18")

    #--- step 2: build an AVS and 1000 matched random AVSs for the input 'risk SNPs'
    bcavs <- avs.preprocess.LDHapMapRel27.hg18(bcarisk, nrand=1000)

    ##################################################
    ### Example of VSE analysis for ERa and FOXA1
    ### cistromes (one genomic annotation each time)
    ##################################################

    #--- step 1: load a precomputed AVS (same 'bcavs' object as above!)
    data(bcavs, package="RTNdata.LDHapMapRel27.hg18")

    #--- step 2: load cistrome data from the Fletcher2013b package
    #NOTE: Fletcher2013b is a large data package, but only two 'bed files'
    #are used to illustrate this analysis (ESR1bdsites and FOXA1bdsites).
    #these bed files provide ERa and FOXA1 binding sites mapped by
    #ChIP-seq experiments
    data(miscellaneous)

    #--- step 3: run the avs.vse pipeline
    bcavs <- avs.vse(bcavs, annotation=ESR1bdsites$bdsites,
                     pValueCutoff=0.001, lab="ERa")
    bcavs <- avs.vse(bcavs, annotation=FOXA1bdsites$bdsites,
                     pValueCutoff=0.001, lab="FOXA1")

    #--- step 4: generate the VSE plots
    avs.plot2(bcavs,"vse",height=2.2)

    ##################################################
    ### Example of VSE analysis for sets of genomic
    ### annotations (e.g. regulons, gene sets, etc.)
    ##################################################

    #--- step 1: load the precomputed AVS (same 'bcavs' object as above!)
    data(bcavs, package="RTNdata.LDHapMapRel27.hg18")

    #--- step 2: load genomic annotation for all genes
    genemap <- as.data.frame(genes(TxDb.Hsapiens.UCSC.hg18.knownGene))
    genemap <- genemap[,c("seqnames","start","end","gene_id")]
    colnames(genemap) <- c("CHROM","START","END","ID")

    #--- step 3: load a TNI object, or any other source of regulons (e.g. gene sets)
    #--- and prepare a gene set list
    #--- (gene ids should be the same as in the 'genemap' object)
    data("rtni1st")
    glist <- tni.get(rtni1st,what="refregulons",idkey="ENTREZ")
    glist <- glist[ c("FOXA1","GATA3","ESR1") ] #reduce the list for demonstration!

    #--- step 4: run the avs.vse pipeline
    bcavs<-avs.vse(bcavs, annotation=genemap, glist=glist, pValueCutoff=0.05)

    #--- step 5: generate the VSE plots
    avs.plot2(bcavs,"vse",height=2.5)

    ### NOTE REGARDING THIS EXAMPLE ####
    #- This example is for demonstration purposes only;
    #- we recommend using the EVSE/eQTL approach when analysing genes/regulons.
    #- Also, the AVS object here is not the same as the one used in the study that
    #- extended the method (doi:10.1038/ng.3458), so the results are not comparable;
    #- (here fewer risk SNPs are considered, and without the eQTL step).
    ####################################

    ## End(Not run)

[Package RTN version 2.11.5 Index]