[BioC] Few errors with edgeR

[Mark Robinson]

Hi Jung,

FYI, posting guide:
http://www.bioconductor.org/help/mailing-list/posting-guide/

Some comments below.

On 18.04.2012, at 15:05, Javerjung Sandhu wrote:

> Dear List,
> I am getting few errors while running edgeR and would like to know few things too as i am novice to edgeR and R. My errors and questions are as follows:
> 1. I am getting an error "could not find function plotBCV". As i am loading package edgeR and limma, I shouldn't be getting this error. So why i am getting this error?

plotBCV() is new and is not present in edgeR 2.4.x ... maybe you are reading the User's guide from edgeR 2.6.x but not using the corresponding software?  Update your R to 2.15 and then re-install edgeR.

> 2. plot smear also doesn't seem to work.
> THIS LAST POINT IS MOST IMPORTANT TO ME.

Your example below throws an error before your call to plotSmear() ...

>> # To generate a plot of the log-fold change against the log-cpm for each tag
>> DETags <- rownames(DGE_List_Object)[as.logical(de)]
> Error in as.logical(de) :
>  cannot coerce type 'closure' to vector of type 'logical'
>> plotSmear(Exact_Test_Results, de.tags = DETags)
> Error in match(de.tags, rownames(object$table)) :
>  object 'DETags' not found
>> abline(h = c(-2,2), col = "dodgerblue")

You don't show how 'de' was defined, and this is a possible source of your problem; I suggest you check this first.


> 3. For the "topTags" function in edgeR, I want to use my own set of p-values ranging from 0.001 upto 0.01 i.e. i want to do multiple testing. Secondly i want to use atleast two adjust methods namely "bonferroni" and "fdr" for my tests. As the p-value i will be giving to the program i also want to get the list of all those genes who pass that threshhold/ cut off value. Because when i choose n = 6449 (total no. of genes in my dataset). It displays everything but actually i want to get only significantly differentially expressed genes. If anybody of you can help me using the "topTags" function properly, i would be really thankful.

I'm afraid I don't understand what you want to do here.  You want to submit your own P-values?

>> TopTags_Object <- topTags(Exact_Test_Results, n = 6449, adj.p.val = P_values,
> +    adjust.method = "bonferroni", sort.by="p.value")
> Error in topTags(Exact_Test_Results, n = 6449, adj.p.val = P_values, adjust.method = "bonferroni",  :
>  unused argument(s) (adj.p.val = P_values)

Your call to togTags() gives an argument 'adj.p.val'; as the error message states, there is no such argument.  Please read the documentation on what this function does -- ?topTags

Hope that helps,
Mark


> Any help would be really appreciated. Thanks in advance.
> The output from the R console is pasted below:
> R version 2.14.1 (2011-12-22)
> Copyright (C) 2011 The R Foundation for Statistical Computing
> ISBN 3-900051-07-0
> Platform: x86_64-pc-mingw32/x64 (64-bit)
> R is free software and comes with ABSOLUTELY NO WARRANTY.
> You are welcome to redistribute it under certain conditions.
> Type 'license()' or 'licence()' for distribution details.
>  Natural language support but running in an English locale
> R is a collaborative project with many contributors.
> Type 'contributors()' for more information and
> 'citation()' on how to cite R or R packages in publications.
> Type 'demo()' for some demos, 'help()' for on-line help, or
> 'help.start()' for an HTML browser interface to help.
> Type 'q()' to quit R.
> [Previously saved workspace restored]
>> #!/usr/bin/env Rscript
>> #source("http://www.bioconductor.org/biocLite.R")
>> #biocLite("edgeR")
>> library(edgeR)
> Loading required package: limma
> Warning messages:
> 1: package ŒedgeR‚ was built under R version 2.14.2
> 2: package Œlimma‚ was built under R version 2.14.2
>> library(limma)
>> options(width=80)
>> 
>> # Reading file
>> Input_file <- read.delim("wt_ko_library.txt",row.names = "Symbol")
>> 
>> # Creating group factor
>> Group <- factor(c(1,1,1,1,2,2,2,2))
>> #Group <- c(rep("wt",4), rep("ko",4))
>> #lib.size <- c(6449,6449,6449,6449,6449,6449,6449,6449)
>> 
>> # Creating DGEList object
>> DGE_List_Object <- DGEList(counts = Input_file, group = Group)
> Calculating library sizes from column totals.
>> # Filtering lowly expressed tags
>> #DGE_List_Object <- DGE_List_Object[rowSums(DGE_List_Object$counts) >= 5,]
>> #DGE_List_Object$samples$lib.size <- colSums(DGE_List_Object$counts)
>> 
>> # Calculating normalization factors
>> DGE_List_Object <- calcNormFactors(DGE_List_Object)
>> 
>> # Creating MDS plot
>> plotMDS(DGE_List_Object, xlim = c(-2,1), labels=colnames(DGE_List_Object))
>> 
>> # Estimating common dispersion
>> DGE_List_Object <- estimateCommonDisp(DGE_List_Object)
>> 
>> # Estimating tagwise dispersion
>> DGE_List_Object <- estimateTagwiseDisp(DGE_List_Object)
>> 
>> # Plotting tagwise dispersion estimates against log-concentration(i.e. tag abundance)
>> plot(log2(DGE_List_Object$conc$conc.common), DGE_List_Object$tagwise.dispersion, panel.first = grid())
>> abline(h = DGE_List_Object$common.dispersion, col = "dodgerblue", lwd = 3)
>> 
>> # Plot BCV
>> plotBCV(DGE_List_Object, xlab = "Abundance (log2 counts per million)", ylab = "Biological coefficient of variation")
> Error: could not find function "plotBCV"
>> abline(h = sqrt(DGE_List_Object$common.dispersion), col = "firebrick", lwd = 3)
>> 
>> # Doing exactTest test
>> Exact_Test_Results <- exactTest(DGE_List_Object)
>> #names(Exact_Test_Results)
>> 
>> # To generate a plot of the log-fold change against the log-cpm for each tag
>> DETags <- rownames(DGE_List_Object)[as.logical(de)]
> Error in as.logical(de) :
>  cannot coerce type 'closure' to vector of type 'logical'
>> plotSmear(Exact_Test_Results, de.tags = DETags)
> Error in match(de.tags, rownames(object$table)) :
>  object 'DETags' not found
>> abline(h = c(-2,2), col = "dodgerblue")
>> 
>> # Creating TopTags Object
>> P_values <- c(0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01)
>> p.adjust(p = P_values, n = length(P_values))
> [1] 0.010 0.018 0.024 0.028 0.030 0.030 0.030 0.030 0.030 0.030
>> TopTags_Object <- topTags(Exact_Test_Results, n = 6449, adj.p.val = P_values,
> +    adjust.method = "bonferroni", sort.by="p.value")
> Error in topTags(Exact_Test_Results, n = 6449, adj.p.val = P_values, adjust.method = "bonferroni",  :
>  unused argument(s) (adj.p.val = P_values)
>> 
>> # Showing the cpm for the genes that edgeR has identified as the most differentially expressed
>> #DEGenes <- rownames(topTags(Exact_Test_Results))
>> #cpm(DGE_List_Object)[DEGenes,]
>> 
>> # DATA TO BE WRITTEN OR WRITING
>> M_of_DGE_List_Object <- as.matrix(DGE_List_Object)
>> write.table(M_of_DGE_List_Object, file = "DGE_List_Object.txt", sep = "\t")
>> write.table(Exact_Test_Results$table, file = "Table_of_Exact_Test_Results.txt", sep = "\t")
>> M_of_ETR_comparison <- as.matrix(Exact_Test_Results$comparison)
>> write.table(M_of_ETR_comparison, file = "Comparison_Component_of_Exact_Test_Results.txt", sep = "\t")
>> write.table(Exact_Test_Results$genes, file = "Genes_of_Exact_Test_Results.txt", sep = "\t")
>> DF_of_TopTags_Object <- as.data.frame(TopTags_Object, row.names = NULL, optional = FALSE)
> Error in as.data.frame(TopTags_Object, row.names = NULL, optional = FALSE) :
>  object 'TopTags_Object' not found
>> write.table(DF_of_TopTags_Object, file = "TopTags_Object.txt", sep = "\t")
> Error in is.data.frame(x) : object 'DF_of_TopTags_Object' not found
>> 
> Thanks,
> Jung
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor

----------
Prof. Dr. Mark Robinson
Bioinformatics
Institute of Molecular Life Sciences
University of Zurich
Winterthurerstrasse 190
8057 Zurich
Switzerland

v: +41 44 635 4848
f: +41 44 635 6898
e: mark.robinson at imls.uzh.ch
o: Y11-J-16
w: http://tiny.cc/mrobin