[med-svn] [r-bioc-cummerbund] 01/06: Imported Upstream version 2.8.0
Andreas Tille
tille at debian.org
Wed Oct 22 11:33:00 UTC 2014
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-bioc-cummerbund.
commit 185e6e78ef1f78d7af90870f4c58893f609c1fd0
Author: Andreas Tille <tille at debian.org>
Date: Wed Oct 22 13:24:26 2014 +0200
Imported Upstream version 2.8.0
---
DESCRIPTION | 23 +--
NAMESPACE | 8 +-
R/AllClasses.R | 3 +-
R/AllGenerics.R | 9 +-
R/database-setup.R | 128 +++++++++++-
R/methods-CuffData.R | 291 ++++++++++++++++++++++-----
R/methods-CuffFeature.R | 48 ++++-
R/methods-CuffFeatureSet.R | 332 ++++++++++++++++++++++++++++---
R/methods-CuffGene.R | 85 +++++++-
R/methods-CuffSet.R | 59 ++++--
R/tools.R | 77 +++++--
build/vignette.rds | Bin 302 -> 303 bytes
data/sampleData.rda | Bin 0 -> 158 bytes
inst/NEWS | 31 +++
inst/doc/cummeRbund-example-workflow.pdf | Bin 198250 -> 198287 bytes
inst/doc/cummeRbund-manual.R | 130 +++++++-----
inst/doc/cummeRbund-manual.Rnw | 59 +++++-
inst/doc/cummeRbund-manual.pdf | Bin 1546156 -> 1563970 bytes
inst/reports/styles/report_style.css | 124 ++++++++++++
man/CuffData-class.Rd | 2 -
man/CuffFeature-class.Rd | 1 +
man/CuffGene-class.Rd | 3 +
man/CuffSet-class.Rd | 7 +-
man/JSdist.Rd | 5 +-
man/QCplots.Rd | 5 +-
man/csCluster.Rd | 5 +-
man/csClusterPlot.Rd | 5 +-
man/csHeatmap.Rd | 13 ++
man/csScatter.Rd | 2 +-
man/diffData.Rd | 7 +
man/dimensionality.Rd | 7 +-
man/dispersionPlot.Rd | 2 +
man/exploratory.Rd | 63 ++++++
man/expressionPlot.Rd | 5 +-
man/featureNames.Rd | 1 -
man/fpkmMatrix.Rd | 1 +
man/getSig.Rd | 8 +-
man/readCufflinks.Rd | 6 +-
man/replicates.Rd | 1 +
vignettes/cummeRbund-manual.Rnw | 59 +++++-
40 files changed, 1390 insertions(+), 225 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 437a5f1..7647173 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,17 +1,15 @@
Package: cummeRbund
Title: Analysis, exploration, manipulation, and visualization of
Cufflinks high-throughput sequencing data.
-Version: 2.6.1
-Date: 2012-05-04
+Version: 2.8.0
+Date: 2013-04-22
Author: L. Goff, C. Trapnell, D. Kelley
-Description: Allows for persistent storage, access, exploration, and
- manipulation of Cufflinks high-throughput sequencing data. In
- addition, provides numerous plotting functions for commonly
- used visualizations.
-Imports: methods, plyr, BiocGenerics
+Description: Allows for persistent storage, access, exploration, and manipulation of Cufflinks high-throughput sequencing data. In addition, provides numerous plotting functions for commonly used visualizations.
+Imports: methods, plyr, BiocGenerics, Biobase
Depends: R (>= 2.7.0), BiocGenerics (>= 0.3.2), RSQLite, ggplot2,
reshape2, fastcluster, rtracklayer, Gviz
-Suggests: cluster, plyr
+Suggests: cluster, plyr, NMFN, stringr, GenomicFeatures, GenomicRanges,
+ rjson
Maintainer: Loyal A. Goff <lgoff at csail.mit.edu>
License: Artistic-2.0
Collate: AllGenerics.R AllClasses.R database-setup.R methods-CuffSet.R
@@ -19,8 +17,7 @@ Collate: AllGenerics.R AllClasses.R database-setup.R methods-CuffSet.R
methods-CuffFeatureSet.R methods-CuffGene.R
methods-CuffFeature.R tools.R
LazyLoad: yes
-biocViews: Sequencing, RNASeq, RNAseqData, GeneExpression,
- DifferentialExpression, Infrastructure, DataImport,
- DataRepresentation, Visualization, Clustering,
- MultipleComparison, QualityControl
-Packaged: 2014-05-30 06:12:09 UTC; biocbuild
+biocViews: Sequencing, RNASeq, GeneExpression, DifferentialExpression,
+ Infrastructure, DataImport, DataRepresentation, Visualization,
+ Clustering, MultipleComparison, QualityControl
+Packaged: 2014-10-14 02:20:08 UTC; biocbuild
diff --git a/NAMESPACE b/NAMESPACE
index bdc55fd..78c96b2 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,6 +4,7 @@ import(methods)
#importFrom(graphics, plot)
importFrom(plyr, defaults)
+importFrom(Biobase, featureNames)
exportClasses(CuffSet,
CuffData,
@@ -16,6 +17,7 @@ exportClasses(CuffSet,
exportMethods(
runInfo,
+ varModel,
annotation,
genes,
isoforms,
@@ -33,6 +35,7 @@ exportMethods(
csVolcano,
csBoxplot,
csHeatmap,
+ csFoldChangeHeatmap,
csDistHeat,
csCluster,
csDendro,
@@ -71,7 +74,10 @@ exportMethods(
csVolcanoMatrix,
MDSplot,
PCAplot,
- fpkmSCVPlot
+ fpkmSCVPlot,
+ csNMF,
+ conditions,
+ csPie
)
export("readCufflinks",
diff --git a/R/AllClasses.R b/R/AllClasses.R
index d3ed8b4..fa98ddd 100644
--- a/R/AllClasses.R
+++ b/R/AllClasses.R
@@ -35,7 +35,8 @@ setClass("CuffDist",
#CuffSet Class is a 'pointer' container to a group of CuffData elements in a cufflinks dataset
setClass("CuffSet",
representation(DB = "SQLiteConnection",
- runInfo = "data.frame",
+ #runInfo = "data.frame",
+ #varModel = "data.frame",
phenoData = "data.frame",
conditions = "data.frame",
genes = "CuffData",
diff --git a/R/AllGenerics.R b/R/AllGenerics.R
index 1cf13dc..ee294ff 100644
--- a/R/AllGenerics.R
+++ b/R/AllGenerics.R
@@ -13,6 +13,7 @@ setGeneric("samples",function(object) standardGeneric("samples"))
setGeneric("replicates",function(object) standardGeneric("replicates"))
setGeneric("conditions",function(object) standardGeneric("conditions"))
setGeneric("runInfo",function(object) standardGeneric("runInfo"))
+setGeneric("varModel",function(object) standardGeneric("varModel"))
setGeneric("genes",function(object) standardGeneric("genes"))
setGeneric("isoforms",function(object) standardGeneric("isoforms"))
setGeneric("TSS",function(object) standardGeneric("TSS"))
@@ -25,7 +26,7 @@ setGeneric("getGenes",function(object,geneIdList,sampleIdList=NULL) standardGene
setGeneric("getGeneId",function(object,idList) standardGeneric("getGeneId"))
setGeneric("findGene",function(object,query) standardGeneric("findGene"))
setGeneric("getFeatures",function(object,featureIdList,sampleIdList=NULL,...) standardGeneric("getFeatures"))
-setGeneric("getSig", function(object,x,y,alpha=0.05,level="genes") standardGeneric("getSig"))
+setGeneric("getSig", function(object,x,y,alpha=0.05,level="genes",method="BH",useCuffMTC=FALSE) standardGeneric("getSig"))
setGeneric("getSigTable", function(object,alpha=0.05,level='genes') standardGeneric("getSigTable"))
setGeneric("addFeatures",function(object, features, ...) standardGeneric("addFeatures"))
setGeneric("findSimilar",function(object,x,n,...) standardGeneric("findSimilar"))
@@ -38,7 +39,7 @@ setGeneric("sigMatrix",function(object,alpha=0.05,level='genes',orderByDist=F) s
#CuffData generics
setGeneric("features",function(object) standardGeneric("features"))
-setGeneric("featureNames",function(object) standardGeneric("featureNames"))
+#setGeneric("featureNames",function(object) standardGeneric("featureNames"))
setGeneric("fpkm",function(object, features=FALSE,...) standardGeneric("fpkm"))
setGeneric("repFpkm",function(object,features=FALSE,...) standardGeneric("repFpkm"))
setGeneric("count",function(object,...) standardGeneric("count"))
@@ -61,6 +62,7 @@ setGeneric("csVolcanoMatrix",function(object,replicates=FALSE,logMode=TRUE,...)
setGeneric("MDSplot",function(object,replicates=FALSE,logMode=TRUE,...) standardGeneric("MDSplot"))
setGeneric("PCAplot",function(object,x="PC1", y="PC2",replicates=TRUE,pseudocount=1.0,scale=TRUE,...) standardGeneric("PCAplot"))
setGeneric("fpkmSCVPlot",function(object,FPKMLowerBound=1,...) standardGeneric("fpkmSCVPlot"))
+setGeneric("csNMF",function(object,k,...) standardGeneric("csNMF"))
#CuffDist generics
setGeneric("distValues",function(object, x, y,...) standardGeneric("distValues"))
@@ -85,6 +87,7 @@ setGeneric("csBoxplot",function(object, logMode=T, ...) standardGeneric("csBoxpl
#CuffFeatureSet
#################
setGeneric("csHeatmap",function(object,rescaling='none', clustering='none', labCol=T, labRow=T, logMode=T, pseudocount=1.0, border=FALSE, heatscale= c(low='darkred',mid='orange',high='white'), heatMidpoint=NULL, ...) standardGeneric("csHeatmap"))
+setGeneric("csFoldChangeHeatmap",function(object, control_condition, replicate_num=NULL, clustering='none', labCol=T, labRow=T, logMode=F, pseudocount=1.0, border=FALSE, heatscale=c(low='steelblue',mid='white',high='tomato'), heatMidpoint=0,fullnames=T,replicates=FALSE,method='none',heatRange=3, ...) standardGeneric("csFoldChangeHeatmap"))
setGeneric("csDistHeat",function(object, replicates=F, samples.not.genes=T, logMode=T, pseudocount=1.0, heatscale=c(low='lightyellow',mid='orange',high='darkred'), heatMidpoint=NULL, ...) standardGeneric("csDistHeat"))
setGeneric("csCluster",function(object, k, logMode=T, pseudocount=1,...) standardGeneric("csCluster"))
#setGeneric("csClusterPlot",function(clustering, pseudocount=1.0) standardGeneric("csClusterPlot"))
@@ -98,6 +101,8 @@ setGeneric("csDendro",function(object,logMode=T,pseudocount=1,replicates=FALSE,.
##################
setGeneric("genePlot",function(object) standardGeneric("genePlot"))
setGeneric("makeGeneRegionTrack",function(object) standardGeneric("makeGeneRegionTrack"))
+setGeneric("as.GRangesList",function(object,f="transcript") standardGeneric("as.GRangesList"))
+setGeneric("csPie",function(object,...) standardGeneric("csPie"))
##############
diff --git a/R/database-setup.R b/R/database-setup.R
index d680fab..fb93544 100644
--- a/R/database-setup.R
+++ b/R/database-setup.R
@@ -56,7 +56,7 @@ loadRepTable<-function(repTableFile,
fileArgs$file = repTableFile
#Read Run Info file
- full = as.data.frame(do.call(read.table,fileArgs))
+ full = as.data.frame(read.delim(repTableFile))
#print(head(full))
#Fix sample_names
@@ -72,6 +72,37 @@ loadRepTable<-function(repTableFile,
bulk_insert(dbConn,insert_SQL,full)
}
+#ReplicateTable
+loadVarModelTable<-function(VarModelFile,
+ dbConn,
+ path,
+ fileArgs = list(sep=sep, header=header, row.names = row.names, quote=quote, na.string=na.string, ...),
+ sep="\t",
+ na.string = "-",
+ header = TRUE,
+ quote = "",
+ stringsAsFactors=FALSE,
+ row.names=NULL,
+ ...) {
+
+ #Setup and reporting
+ write(paste("Reading Var Model Info ",VarModelFile,sep=""),stderr())
+ fileArgs$file = VarModelFile
+
+ #Read Run Info file
+ full = as.data.frame(read.delim(VarModelFile))
+ #print(head(full))
+
+ #Fix sample_names
+ full$condition<-make.db.names(dbConn,as.character(full$condition),unique=FALSE)
+
+
+ #Load into database (replicates table)
+ write("Writing varModel Table",stderr())
+ insert_SQL<-'INSERT INTO varModel VALUES(:condition, :locus, :compatible_count_mean, :compatible_count_var, :total_count_mean, :total_count_var, :fitted_var)'
+ bulk_insert(dbConn,insert_SQL,full)
+}
+
#Genes
loadGenes<-function(fpkmFile,
diffFile,
@@ -1580,6 +1611,23 @@ CREATE TABLE "isoformReplicateData"(
FOREIGN KEY("sample_name")
REFERENCES "samples"("sample_name")
);
+DROP TABLE IF EXISTS "varModel";
+CREATE TABLE "varModel"(
+ "condition" VARCHAR(45) NOT NULL,
+ "locus" VARCHAR(45) NOT NULL,
+ "compatible_count_mean" FLOAT,
+ "compatible_count_var" FLOAT,
+ "total_count_mean" FLOAT,
+ "total_count_var" FLOAT,
+ "fitted_var" FLOAT,
+ CONSTRAINT "fk_varModel_locus1"
+ FOREIGN KEY("locus")
+ REFERENCES "genes"("locus")
+ CONSTRAINT "fk_varModel_condition1"
+ FOREIGN KEY("condition")
+ REFERENCES "samples"("sample_name")
+);
+
COMMIT;
@@ -1678,6 +1726,8 @@ CREATE INDEX "features.features_strand_index" ON "features"("strand");
CREATE INDEX "features.features_start_end_index" ON "features"("start","end");
CREATE INDEX "features.fk_features_genes1" ON "features"("gene_id");
CREATE INDEX "features.fk_features_isoforms1" ON "features"("isoform_id");
+CREATE INDEX "varModel.varModel_condition1" ON "varModel"("condition");
+CREATE INDEX "varModel.varModel_locus1" ON "varModel"("locus");
'
create.sql <- strsplit(index.text,"\n")[[1]]
@@ -1741,6 +1791,7 @@ readCufflinks<-function(dir = getwd(),
CDSDiff="cds.diff",
promoterFile="promoters.diff",
splicingFile="splicing.diff",
+ varModelFile="var_model.info",
driver = "SQLite",
genome = NULL,
rebuild = FALSE,
@@ -1770,6 +1821,7 @@ readCufflinks<-function(dir = getwd(),
CDSDiff=file.path(dir,CDSDiff)
promoterFile=file.path(dir,promoterFile)
splicingFile=file.path(dir,splicingFile)
+ varModelFile=file.path(dir,varModelFile)
#Check to see whether dbFile exists
@@ -1786,6 +1838,11 @@ readCufflinks<-function(dir = getwd(),
if(file.exists(repTableFile)){
loadRepTable(repTableFile,dbConn)
}
+
+ if(file.exists(varModelFile)){
+ loadVarModelTable(varModelFile,dbConn)
+ }
+
if(!is.null(gtfFile)){
if(!is.null(genome)){
.loadGTF(gtfFile,genome,dbConn)
@@ -1908,6 +1965,73 @@ readCufflinks<-function(dir = getwd(),
dbSendQuery(dbConn,genomeInsertQuery)
}
+.cuff2db <- function(gtfFile, out.file = NULL, verbose = TRUE) {
+
+ #require(rtracklayer)
+ #require(GenomicRanges)
+ #require(GenomicFeatures)
+
+ requiredAttribs <- c("gene_id", "transcript_id", "exon_number")
+
+ if (verbose) message("Importing ", gtfFile)
+ tmp <- import(gtfFile, asRangedData=FALSE)
+
+ #dispose of unspliced unstranded transcripts
+ #tmp <- tmp[ which(strand(tmp) %in% c('+','-')) ]
+
+ # fix the gene IDs
+ #values(tmp)$gene_id <- gsub('CUFF.', '', values(tmp)$gene_id)
+
+ # fix the exon IDs
+ #values(tmp)$transcript_id <- gsub('CUFF.', '', values(tmp)$transcript_id)
+
+ # split the object into transcript and exon pieces
+ #by.type = split(tmp, values(tmp)$type)
+ #browser()
+
+ #make transcripts table
+ tmpT <- split(tmp,
+ values(tmp)$transcript_id)
+ if(verbose) message('Attempting to create the transcripts data.frame')
+ transcripts <- data.frame(
+ tx_id=1:length(tmpT),
+ tx_name=names(tmpT),
+ tx_chrom=as.character(seqnames(unlist(tmpT))[start(tmpT at partitioning)]),
+ tx_strand=as.character(strand(unlist(tmpT))[start(tmpT at partitioning)]),
+ tx_start=sapply(start(ranges(tmpT)), min),
+ tx_end=sapply(end(ranges(tmpT)), max),
+ stringsAsFactors=FALSE
+ )
+
+ #make splicings table
+ tmpS <- split(tmp, values(tmp)$transcript_id)
+ if(verbose) message('Attempting to create the splicings data.frame')
+ splicings <- data.frame(
+ tx_id=rep(1:length(tmpS), elementLengths(tmpS)),
+ exon_rank=as.integer(values(unlist(tmpS))$exon_number),
+ exon_chrom=as.character(seqnames(unlist(tmpS))),
+ exon_strand=as.character(strand(unlist(tmpS))),
+ exon_start=start(unlist(tmpS)),
+ exon_end=end(unlist(tmpS)),
+ stringsAsFactors=FALSE
+ )
+
+ #make genes table
+ if(verbose) message('Attempting to create the genes data.frame')
+ gene_txs <- tapply(values(tmp)$transcript_id, values(tmp)$gene_id, unique)
+ genes <- data.frame(
+ tx_name=unlist(gene_txs),
+ gene_id=rep(names(gene_txs), sapply(gene_txs, length)),
+ stringsAsFactors=FALSE)
+
+ #create the db
+ if (verbose) message("Creating TxDb")
+ tmpdb <- makeTranscriptDb(transcripts, splicings, genes=genes)
+ if (verbose) message("Use saveFeatures() to save the database to a file")
+ return(tmpdb)
+
+}
+
#library(Gviz)
#myGeneId<-'XLOC_000071'
#geneQuery<-paste("SELECT start,end,source AS feature,gene_id as gene,exon_number AS exon,transcript_id as transcript,gene_name as symbol, exon_number as rank, strand FROM features WHERE gene_id ='",myGeneId,"'",sep="")
@@ -1941,4 +2065,4 @@ readCufflinks<-function(dir = getwd(),
#loadIsoforms("isoforms.fpkm_tracking","isoform_exp.diff",dbConn)
#loadTSS("tss_groups.fpkm_tracking","tss_group_exp.diff",dbConn)
#loadCDS("cds.fpkm_tracking","cds_exp.diff",dbConn)
-#date()
\ No newline at end of file
+#date()
diff --git a/R/methods-CuffData.R b/R/methods-CuffData.R
index 6564709..6a69ee6 100644
--- a/R/methods-CuffData.R
+++ b/R/methods-CuffData.R
@@ -72,7 +72,7 @@ setMethod("addFeatures",signature="CuffData",.addFeatures)
#Accessors
###################
.annotation<-function(object){
- featureQuery<-paste("SELECT * FROM ",object at tables$mainTable," x LEFT JOIN ",object at tables$featureTable," xf USING (",object at idField,")",sep="")
+ featureQuery<-paste("SELECT * FROM ",object at tables$mainTable," x LEFT JOIN features xf USING (",object at idField,")",sep="")
dbGetQuery(object at DB, featureQuery)
}
@@ -124,7 +124,7 @@ setMethod("replicates","CuffData",.replicates)
if(!features){
FPKMQuery<-paste("SELECT * FROM ",object at tables$dataTable," WHERE sample_name IN ",sampleString,sep="")
}else{
- FPKMQuery<-paste("SELECT xf.*,xm.*,x.sample_name,x.fpkm,x.conf_hi,x.conf_lo FROM ",object at tables$dataTable," x LEFT JOIN ",object at tables$featureTable," xf ON x.",object at idField,"=xf.",object at idField," LEFT JOIN ",object at tables$mainTable," xm ON x.",object at idField,"=xm.",object at idField," WHERE x.sample_name IN ",sampleString,sep="")
+ FPKMQuery<-paste("SELECT xf.*,xm.*,x.sample_name,x.fpkm,x.conf_hi,x.conf_lo FROM ",object at tables$dataTable," x LEFT JOIN features xf ON x.",object at idField,"=xf.",object at idField," LEFT JOIN ",object at tables$mainTable," xm ON x.",object at idField,"=xm.",object at idField," WHERE x.sample_name IN ",sampleString,sep="")
#print(FPKMQuery)
}
res<-dbGetQuery(object at DB,FPKMQuery)
@@ -157,7 +157,7 @@ setMethod("fpkm","CuffData",.fpkm)
if(!features){
FPKMQuery<-paste("SELECT * FROM ",object at tables$replicateTable," WHERE rep_name IN ",sampleString,sep="")
}else{
- FPKMQuery<-paste("SELECT xf.*,xm.*,x.rep_name,x.raw_frags,x.internal_scaled_frags,x.external_scaled_frags,x.fpkm,x.effective_length,x.status FROM ",object at tables$replicateTable," x LEFT JOIN ",object at tables$featureTable," xf on x.",object at idField,"=xf.",object at idField," LEFT JOIN ",object at tables$mainTable," xm ON x.",object at idField,"=xm.",object at idField," WHERE x.rep_name IN ",sampleString,sep="")
+ FPKMQuery<-paste("SELECT xf.*,xm.*,x.rep_name,x.raw_frags,x.internal_scaled_frags,x.external_scaled_frags,x.fpkm,x.effective_length,x.status FROM ",object at tables$replicateTable," x LEFT JOIN features xf on x.",object at idField,"=xf.",object at idField," LEFT JOIN ",object at tables$mainTable," xm ON x.",object at idField,"=xm.",object at idField," WHERE x.rep_name IN ",sampleString,sep="")
}
#print(FPKMQuery)
res<-dbGetQuery(object at DB,FPKMQuery)
@@ -211,14 +211,19 @@ setMethod("count","CuffData",.count)
}
samp<-samples(object)
- FPKMMatQuery<-paste("select x.",object at idField,", ",sep="")
+ FPKMMatQuery<-paste("select x.",object at idField,", x.gene_short_name, ",sep="")
for (i in samp){
FPKMMatQuery<-paste(FPKMMatQuery,"sum(case when xd.sample_name ='",i,"' then fpkm end) as ",i,",",sep="")
}
FPKMMatQuery<-substr(FPKMMatQuery, 1, nchar(FPKMMatQuery)-1)
FPKMMatQuery<-paste(FPKMMatQuery," from ",object at tables$mainTable," x LEFT JOIN ",object at tables$dataTable," xd on x.",object at idField," = xd.",object at idField," group by x.",object at idField,sep="")
res<-dbGetQuery(object at DB,FPKMMatQuery)
- res<-data.frame(res[,-1],row.names=res[,1])
+ if(fullnames){
+ res<-data.frame(res[,-c(1,2)],row.names=paste(res[,2],res[,1],sep="|"))
+ }else{
+ res<-data.frame(res[,-c(1,2)],row.names=res[,1])
+ }
+
if(!missing(sampleIdList)){
res<-data.frame(res[,sampleIdList],row.names=rownames(res))
colnames(res)<-sampleIdList
@@ -241,14 +246,19 @@ setMethod("fpkmMatrix","CuffData",.fpkmMatrix)
}
samp<-replicates(object)
- FPKMMatQuery<-paste("select x.",object at idField,", ",sep="")
+ FPKMMatQuery<-paste("select x.",object at idField,", x.gene_short_name, ",sep="")
for (i in samp){
FPKMMatQuery<-paste(FPKMMatQuery,"sum(case when xd.rep_name ='",i,"' then fpkm end) as ",i,",",sep="")
}
FPKMMatQuery<-substr(FPKMMatQuery, 1, nchar(FPKMMatQuery)-1)
FPKMMatQuery<-paste(FPKMMatQuery," from ",object at tables$mainTable," x LEFT JOIN ",object at tables$replicateTable," xd on x.",object at idField," = xd.",object at idField," group by x.",object at idField,sep="")
res<-dbGetQuery(object at DB,FPKMMatQuery)
- res<-data.frame(res[,-1],row.names=res[,1])
+
+ if(fullnames){
+ res<-data.frame(res[,-c(1,2)],row.names=paste(res[,2],res[,1],sep="|"))
+ }else{
+ res<-data.frame(res[,-c(1,2)],row.names=res[,1])
+ }
if(!missing(repIdList)){
res<-data.frame(res[,repIdList],row.names=rownames(res))
colnames(res)<-repIdList
@@ -271,14 +281,18 @@ setMethod("repFpkmMatrix","CuffData",.repFpkmMatrix)
}
samp<-samples(object)
- CountMatQuery<-paste("select x.",object at idField,", ",sep="")
+ CountMatQuery<-paste("select x.",object at idField,", x.gene_short_name, ",sep="")
for (i in samp){
CountMatQuery<-paste(CountMatQuery,"sum(case when xd.sample_name ='",i,"' then count end) as ",i,",",sep="")
}
CountMatQuery<-substr(CountMatQuery, 1, nchar(CountMatQuery)-1)
CountMatQuery<-paste(CountMatQuery," from ",object at tables$mainTable," x LEFT JOIN ",object at tables$countTable," xd on x.",object at idField," = xd.",object at idField," group by x.",object at idField,sep="")
res<-dbGetQuery(object at DB,CountMatQuery)
- res<-data.frame(res[,-1],row.names=res[,1])
+ if(fullnames){
+ res<-data.frame(res[,-c(1,2)],row.names=paste(res[,2],res[,1],sep="|"))
+ }else{
+ res<-data.frame(res[,-c(1,2)],row.names=res[,1])
+ }
if(!missing(sampleIdList)){
res<-data.frame(res[,sampleIdList],row.names=rownames(res))
colnames(res)<-sampleIdList
@@ -300,14 +314,18 @@ setMethod("countMatrix","CuffData",.countMatrix)
myLevels<-getRepLevels(object)
}
reps<-replicates(object)
- repCountMatQuery<-paste("select x.",object at idField,", ",sep="")
+ repCountMatQuery<-paste("select x.",object at idField,", x.gene_short_name, ",sep="")
for (i in reps){
repCountMatQuery<-paste(repCountMatQuery,"sum(case when xr.rep_name ='",i,"' then external_scaled_frags end) as ",i,",",sep="")
}
repCountMatQuery<-substr(repCountMatQuery, 1, nchar(repCountMatQuery)-1)
repCountMatQuery<-paste(repCountMatQuery," from ",object at tables$mainTable," x LEFT JOIN ",object at tables$replicateTable," xr on x.",object at idField," = xr.",object at idField," group by x.",object at idField,sep="")
res<-dbGetQuery(object at DB,repCountMatQuery)
- res<-data.frame(res[,-1],row.names=res[,1])
+ if(fullnames){
+ res<-data.frame(res[,-c(1,2)],row.names=paste(res[,2],res[,1],sep="|"))
+ }else{
+ res<-data.frame(res[,-c(1,2)],row.names=res[,1])
+ }
if(!missing(repIdList)){
res<-data.frame(res[,repIdList],row.names=rownames(res))
colnames(res)<-repIdList
@@ -331,7 +349,7 @@ setMethod("repCountMatrix","CuffData",.repCountMatrix)
if(!features){
diffQuery<-paste("SELECT * FROM ",object at tables$expDiffTable,sep="")
}else{
- diffQuery<-paste("SELECT xm.*, xed.*, xf.* FROM ",object at tables$mainTable," xm LEFT JOIN ",object at tables$expDiffTable," xed ON xm.",object at idField,"=xed.",object at idField," LEFT JOIN ",object at tables$featureTable," xf ON xm.",object at idField,"=xf.",object at idField,sep="")
+ diffQuery<-paste("SELECT xm.*, xed.*, xf.* FROM ",object at tables$mainTable," xm LEFT JOIN ",object at tables$expDiffTable," xed ON xm.",object at idField,"=xed.",object at idField," LEFT JOIN features xf ON xm.",object at idField,"=xf.",object at idField,sep="")
}
}else if (missing(x) || missing(y)){
stop("You must supply both x and y or neither.")
@@ -339,7 +357,7 @@ setMethod("repCountMatrix","CuffData",.repCountMatrix)
if(!features){
diffQuery<-paste("SELECT x.",object at idField,", xed.* FROM ",object at tables$mainTable," x LEFT JOIN ",object at tables$expDiffTable," xed on x.",object at idField," = xed.",object at idField," WHERE ((sample_1 = '",x,"' AND sample_2 = '",y,"') OR (sample_1 = '",y,"' AND sample_2 = '",x,"'))",sep="")
}else{
- diffQuery<-paste("SELECT xm.*, xed.*, xf.* FROM ",object at tables$mainTable," xm LEFT JOIN ",object at tables$expDiffTable," xed on xm.",object at idField," = xed.",object at idField," LEFT JOIN ",object at tables$featureTable," xf ON xm.",object at idField,"=xf.",object at idField," WHERE ((sample_1 = '",x,"' AND sample_2 = '",y,"') OR (sample_1 = '",y,"' AND sample_2 = '",x,"'))",sep="")
+ diffQuery<-paste("SELECT xm.*, xed.*, xf.* FROM ",object at tables$mainTable," xm LEFT JOIN ",object at tables$expDiffTable," xed on xm.",object at idField," = xed.",object at idField," LEFT JOIN features xf ON xm.",object at idField,"=xf.",object at idField," WHERE ((sample_1 = '",x,"' AND sample_2 = '",y,"') OR (sample_1 = '",y,"' AND sample_2 = '",x,"'))",sep="")
}
}
dat<-dbGetQuery(object at DB,diffQuery)
@@ -459,6 +477,12 @@ setMethod("getLevels",signature(object="CuffData"),.getLevels)
setMethod("getRepLevels",signature(object="CuffData"),.getRepLevels)
+.getRepConditionLevels<-function(object){
+ levelsQuery<-'SELECT r.sample_name FROM replicates r JOIN samples s ON r.sample_name=s.sample_name ORDER BY s.sample_index ASC'
+ levels<-dbGetQuery(object at DB,levelsQuery)$sample_name
+ levels
+}
+
#Useful SQL commands
#SELECT g.gene_id, g.class_code, g.nearest_ref_id, g.gene_short_name, g.locus, g.length, g.coverage, g.status, gd.sample_name, gd.fpkm, gd.conf_hi, gd.conf_lo FROM genes g LEFT JOIN geneData gd ON g.gene_id = gd.gene_id WHERE (g.gene_id = 'XLOC_000001');
@@ -505,7 +529,7 @@ setMethod("getRepLevels",signature(object="CuffData"),.getRepLevels)
setMethod("csDensity",signature(object="CuffData"),.density)
.scatter<-function(object,x,y,logMode=TRUE,pseudocount=1.0,labels,smooth=FALSE,colorByStatus=FALSE, drawRug=TRUE, ...){
- dat<-fpkmMatrix(object)
+ dat<-fpkmMatrix(object,fullnames=TRUE)
samp<-samples(object)
#check to make sure x and y are in samples
@@ -520,6 +544,16 @@ setMethod("csDensity",signature(object="CuffData"),.density)
}
}
+ #Attach tracking_id and gene_short_name
+ if(!missing(labels)){
+ require(stringr)
+ tracking<-str_split_fixed(rownames(dat),"\\|",2)
+ dat$gene_short_name<-tracking[,1]
+ dat$tracking_id<-tracking[,2]
+
+ labeled.dat<-dat[dat$gene_short_name %in% labels,]
+ }
+
#make plot object
p<-ggplot(dat)
p<- p + aes_string(x=x,y=y)
@@ -544,15 +578,11 @@ setMethod("csDensity",signature(object="CuffData"),.density)
}
#Add highlights from labels
-# if(!missing(labels)){
-# labelIdx<-fData(object)$gene_short_name %in% labels
-# labelfp<-fp[labelIdx,]
-# labelfp$gene_short_name<-fData(object)$gene_short_name[labelIdx]
-# #print(head(labelfp))
-# p <- p + geom_point(data=labelfp,size=1.2,color="red")
-# p <- p + geom_text(data=labelfp,aes(label=gene_short_name),color="red",hjust=0,vjust=0,angle=45,size=2)
-# }
-#
+ if(!missing(labels)){
+ p <- p + geom_point(data=labeled.dat,aes_string(x=x,y=y),size=1.3,color="red")
+ p <- p + geom_text(data=labeled.dat,aes_string(x=x,y=y,label='gene_short_name'),color="red",hjust=0,vjust=0,angle=0,size=4)
+ }
+
#logMode
if(logMode){
p <- p + scale_y_log10() + scale_x_log10()
@@ -565,6 +595,80 @@ setMethod("csDensity",signature(object="CuffData"),.density)
setMethod("csScatter",signature(object="CuffData"), .scatter)
+.scatter2<-function(object,x,y,logMode=TRUE,pseudocount=1.0,labels,smooth=FALSE,alpha=0.05,colorByStatus=FALSE, drawRug=TRUE, ...){
+ samp<-samples(object)
+
+ #check to make sure x and y are in samples
+ if (!all(c(x,y) %in% samp)){
+ stop("One or more values of 'x' or 'y' are not valid sample names!")
+ }
+
+ #Setup query string
+ sampleList<-paste("('",x,"','",y,"')",sep="")
+
+ scatterQuery<-paste("SELECT m.gene_short_name,d.",object at idField,",sum(CASE WHEN d.sample_name='",x,"' THEN d.fpkm END) as x,sum(CASE WHEN d.sample_name='",y,"' THEN d.fpkm END) as y, edd.status, edd.p_value, edd.q_value FROM ",object at tables$mainTable," m LEFT JOIN ",object at tables$dataTable," d ON m.",object at idField,"=d.",object at idField," LEFT JOIN ",object at tables$expDiffTable," edd ON m.",object at idField," = edd.",object at idField," WHERE (edd.sample_1 in ",sampleList," AND (edd.sample_2 [...]
+ #write(scatterQuery,stderr())
+
+ #Retrieve data
+ dat<-dbGetQuery(object at DB,scatterQuery)
+
+ #add pseudocount if necessary
+ if(logMode){
+ for (i in samp){
+ dat$x<-dat$x+pseudocount
+ dat$y<-dat$y+pseudocount
+ }
+ }
+
+ #Flag significant genes
+ dat$significant<-"no"
+ dat$significant[dat$q_value<=alpha]<-"yes"
+ dat$significant<-factor(dat$significant,levels=c("no","yes"))
+
+ #Attach tracking_id and gene_short_name
+ if(!missing(labels)){
+ labeled.dat<-dat[dat$gene_short_name %in% labels,]
+ }
+
+ #make plot object
+ p<-ggplot(dat) + theme_bw()
+ p<- p + aes_string(x='x',y='y')
+
+ #Right now, this does nothing, because 'significant' is not returned from fpkmMatrix object so I don't have this as a feature to draw
+ if(colorByStatus){
+ p<- p + geom_point(aes(color=significant),size=1.2) + scale_color_manual(values=c("black","red"))
+ }else{
+ p<- p + geom_point(size=1.2,alpha=I(1/3))
+ }
+ #Add symmetry line
+ p<- p + geom_abline(intercept=0,slope=1,linetype=2)
+
+ #Add rug
+ if(drawRug){
+ p<- p + geom_rug(size=0.8,alpha=0.01)
+ }
+
+ #add smoother
+ if(smooth){
+ p <- p + stat_smooth(method="lm",fill="blue",alpha=0.2)
+ }
+
+ #Add highlights from labels
+ if(!missing(labels)){
+ p <- p + geom_point(data=labeled.dat,aes_string(x='x',y='y'),size=1.3,color="red")
+ p <- p + geom_text(data=labeled.dat,aes_string(x='x',y='y',label='gene_short_name'),color="red",hjust=0,vjust=0,angle=0,size=4)
+ }
+
+ #logMode
+ if(logMode){
+ p <- p + scale_y_log10() + scale_x_log10()
+ }
+
+ #Add title & Return value
+ p<- p + labs(title=object at tables$mainTable,x=x,y=y)
+ p
+}
+
.scatterMat<-function(object,replicates=FALSE,logMode=TRUE,pseudocount=1.0,hexbin=FALSE,useCounts=FALSE,...){
if(replicates){
if(useCounts){
@@ -579,8 +683,9 @@ setMethod("csScatter",signature(object="CuffData"), .scatter)
dat<-fpkmMatrix(object)
}
}
-
- if(hexbin){
+
+ if(logMode)
+ {
dat<-dat+pseudocount
}
@@ -627,12 +732,13 @@ setMethod("csScatterMatrix",signature(object="CuffData"),.scatterMat)
s2<-unique(dat$sample_2)
p<-ggplot(dat)
- if(showSignificant){
- p<- p + geom_point(aes(x=log2_fold_change,y=-log10(p_value),color=significant),size=0.8)
- }else{
+ if(showSignificant==FALSE){
p<- p + geom_point(aes(x=log2_fold_change,y=-log10(p_value)),size=1.2)
+ }else{
+ p<- p + geom_point(aes(x=log2_fold_change,y=-log10(p_value),color=significant),size=1.2)
}
#Add title and return
+ p<- p + theme_bw()
p<- p + labs(title=paste(object at tables$mainTable,": ",s2,"/",s1,sep=""))
p<- p + scale_colour_manual(values = c("black","red"))
@@ -642,6 +748,7 @@ setMethod("csScatterMatrix",signature(object="CuffData"),.scatterMat)
p <- p + xlab(bquote(paste(log[2],"(fold change)",sep=""))) +
ylab(bquote(paste(-log[10],"(p value)",sep="")))
p
+
}
setMethod("csVolcano",signature(object="CuffData"), .volcano)
@@ -667,7 +774,7 @@ setMethod("csVolcano",signature(object="CuffData"), .volcano)
mapping <- defaults(mapping, aes_string(x = "log2_fold_change", y = "-log10(p_value)", color="significant"))
class(mapping) <- "uneval"
- p <-ggplot(dat) + geom_point(mapping,na.rm=TRUE,size=0.8) + scale_colour_manual(values = c("black","red")) + geom_text(aes(x=0,y=15,label=label),data=filler) + facet_grid(sample_1~sample_2)
+ p <-ggplot(dat) + geom_point(mapping,na.rm=TRUE,size=0.8) + scale_colour_manual(values = c("black","red")) + facet_grid(sample_1~sample_2)
p<- p + geom_vline(aes(x=0),linetype=2)
@@ -679,7 +786,7 @@ setMethod("csVolcano",signature(object="CuffData"), .volcano)
setMethod("csVolcanoMatrix",signature(object="CuffData"),.volcanoMatrix)
-.distheat<-function(object, replicates=F, samples.not.genes=T, logMode=T, pseudocount=1.0, heatscale=c(low='lightyellow',mid='orange',high='darkred'), heatMidpoint=NULL, ...) {
+.distheat<-function(object, replicates=F, samples.not.genes=T, logMode=T, pseudocount=1.0, heatscale=c(low='lightyellow',mid='orange',high='darkred'), heatMidpoint=NULL, sigDigits=3, ...) {
# get expression from a sample or gene perspective
if(replicates){
obj.fpkm<-repFpkmMatrix(object,fullnames=T)
@@ -706,7 +813,7 @@ setMethod("csVolcanoMatrix",signature(object="CuffData"),.volcanoMatrix)
# make data frame
dist.df = melt(as.matrix(obj.dists),varnames=c("X1","X2"))
-
+ dist.df$value<-as.numeric(format(dist.df$value,digits=sigDigits))
# initialize
g = ggplot(dist.df, aes(x=X1, y=X2, fill=value))
@@ -731,7 +838,7 @@ setMethod("csVolcanoMatrix",signature(object="CuffData"),.volcanoMatrix)
g = g + scale_fill_gradient2(low=heatscale[1], mid=heatscale[2], high=heatscale[3], midpoint=heatMidpoint, name="JS Distance")
}
if(samples.not.genes){
- g <- g + geom_text(aes(label=format(value,digits=3)))
+ g <- g + geom_text(aes(label=value))
}
# return
g
@@ -819,6 +926,7 @@ setMethod("dispersionPlot",signature(object="CuffData"),.dispersionPlot)
.MDSplot<-function(object,replicates=FALSE,logMode=TRUE,pseudocount=1.0){
if(replicates){
dat<-repFpkmMatrix(object)
+ #repData<-sapply(replicates(object),function(x){strsplit(x,"_")[[1]][1]}) This is to color by condition and not replicate...
}else{
dat<-fpkmMatrix(object)
}
@@ -838,7 +946,7 @@ setMethod("dispersionPlot",signature(object="CuffData"),.dispersionPlot)
setMethod("MDSplot",signature(object="CuffData"),.MDSplot)
#Not sure if I want to include this or not..
-.PCAplot<-function(object,x="PC1", y="PC2",replicates=FALSE,pseudocount=1.0,scale=TRUE,...){
+.PCAplot<-function(object,x="PC1", y="PC2",replicates=FALSE,pseudocount=1.0,scale=TRUE,showPoints=TRUE,...){
if(replicates){
fpkms<-repFpkmMatrix(object)
}else{
@@ -849,7 +957,10 @@ setMethod("MDSplot",signature(object="CuffData"),.MDSplot)
dat <- data.frame(obsnames=row.names(PC$x), PC$x)
#dat$shoutout<-""
#dat$shoutout[matchpt(PC$rotation,PC$x)$index]<-rownames(pca$x[matchpt(pca$rotation,pca$x)$index,])
- plot <- ggplot(dat, aes_string(x=x, y=y)) + geom_point(alpha=.4, size=0.8, aes(label=obsnames))
+ plot <- ggplot(dat, aes_string(x=x, y=y))
+ if(showPoints){
+ plot<- plot + geom_point(alpha=.4, size=0.8, aes(label=obsnames))
+ }
plot <- plot + geom_hline(aes(0), size=.2) + geom_vline(aes(0), size=.2) #+ geom_text(aes(label=shoutout),size=2,color="red")
datapc <- data.frame(varnames=rownames(PC$rotation), PC$rotation)
mult <- min(
@@ -862,8 +973,8 @@ setMethod("MDSplot",signature(object="CuffData"),.MDSplot)
)
plot <- plot +
#coord_equal() +
- geom_text(data=datapc, aes(x=v1, y=v2, label=varnames), size = 3, vjust=1, color="red")
- plot <- plot + geom_segment(data=datapc, aes(x=0, y=0, xend=v1, yend=v2), arrow=arrow(length=unit(0.2,"cm")), alpha=0.75, color="red") + theme_bw()
+ geom_text(data=datapc, aes(x=v1, y=v2, label=varnames,color=varnames), vjust=1)
+ plot <- plot + geom_segment(data=datapc, aes(x=0, y=0, xend=v1, yend=v2,color=varnames), arrow=arrow(length=unit(0.2,"cm")), alpha=0.75) + theme_bw()
plot
}
@@ -893,15 +1004,27 @@ setMethod('PCAplot',signature(object="CuffData"),.PCAplot)
p
}
-.fpkmSCVPlot<-function(object,FPKMLowerBound=1){
+.fpkmSCVPlot<-function(object,FPKMLowerBound=1,showPool=FALSE){
+ if(!showPool){
+ showPool=any(table(.getRepConditionLevels(object))<2) # Counts number of replicates per condition and requires a minimum of 2. If any n<2, values are determined from pooling across samples.
+ if(showPool){
+ warning("At least one of your conditions does not have enough replicates to estimate variance. Estimating variance across all conditions instead.")
+ }
+ }
dat<-repFpkm(object)
colnames(dat)[1]<-"tracking_id"
dat<-dat[,c('tracking_id','sample_name','fpkm')]
dat<-dat[dat$fpkm>0,]
+ if(showPool){
+ subsetCols<-c('tracking_id')
+ }else{
+ subsetCols<-c('tracking_id','sample_name')
+ }
+
#Option 3 (tapply on log10(replicateFPKM) values)
- dat.means<-tapply(dat$fpkm,dat[,c('tracking_id','sample_name')],function(x){mean(x,na.rm=T)})
- dat.sd<-tapply(dat$fpkm,dat[,c('tracking_id','sample_name')],function(x){sd(x,na.rm=T)})
+ dat.means<-tapply(dat$fpkm,dat[,subsetCols],function(x){mean(x,na.rm=T)})
+ dat.sd<-tapply(dat$fpkm,dat[,subsetCols],function(x){sd(x,na.rm=T)})
#write("Calculating replicate fpkm mean...",stderr())
dat.means<-melt(dat.means)
#write("Calculating replicate fpkm stdev...",stderr())
@@ -912,16 +1035,22 @@ setMethod('PCAplot',signature(object="CuffData"),.PCAplot)
colnames(dat)[colnames(dat)=="dat.sd$stdev"]<-'stdev'
dat<-dat[!is.na(dat$stdev) & !is.na(dat$fpkm),]
dat<-dat[dat$fpkm>0 & dat$stdev>0,]
- dat$sample_name<-factor(dat$sample_name,levels=samples(object))
+ if(!showPool){
+ dat$sample_name<-factor(dat$sample_name,levels=samples(object))
+ }
p <-ggplot(dat,aes(x=fpkm,y=(stdev/fpkm)^2),na.rm=T)
#p <-ggplot(dat,aes(x=log10(fpkm+1),y=log10(stdev)),na.rm=T)
- p <- p + #geom_point(aes(color=sample_name),size=1,na.rm=T) +
- stat_smooth(aes(color=sample_name,fill=sample_name),na.rm=T,method='auto',fullrange=T) +
- scale_x_log10() +
- scale_y_continuous(name=bquote(CV^2)) +
- xlab(bquote(paste(log[10],"FPKM",sep=" "))) +
- theme_bw() + xlim(c(log10(FPKMLowerBound),max(log10(dat$fpkm)))) + labs(title=object at type)
+ if(showPool){
+ p <- p + stat_smooth(na.rm=T,method='auto',fullrange=T)
+ }else{
+ p <- p + #geom_point(aes(color=sample_name),size=1,na.rm=T) +
+ stat_smooth(aes(color=sample_name,fill=sample_name),na.rm=T,method='auto',fullrange=T)
+ }
+ p <- p + scale_x_log10() +
+ scale_y_continuous(name=bquote(CV^2)) +
+ xlab(bquote(paste(log[10],"FPKM",sep=" "))) +
+ theme_bw() + xlim(c(log10(FPKMLowerBound),max(log10(dat$fpkm)))) + labs(title=object at type)
p
}
@@ -964,6 +1093,31 @@ setMethod("fpkmSCVPlot",signature(object="CuffData"),.fpkmSCVPlot)
setMethod("csSpecificity",signature(object="CuffData"),.specificity)
+######################
+# Exploratory Analysis
+######################
+
+.nmf<-function(object,k,logMode=T,pseudocount=1,maxiter=1000,replicates=FALSE,fullnames=FALSE){
+ require(NMFN)
+ if(missing(k)) stop("Please provide a rank value for factorization (arg=k)")
+
+ if(replicates){
+ m=repFpkmMatrix(object,fullnames=fullnames)
+ }else{
+ m=fpkmMatrix(object,fullnames=fullnames)
+ }
+
+ if(logMode)
+ {
+ m = log10(m+pseudocount)
+ }
+
+ myNMF<-nnmf(m,k=k,maxiter=maxiter)
+ return (myNMF)
+}
+
+setMethod("csNMF",signature(object="CuffData"),.nmf)
+
#############
# GSEA helper methods
#############
@@ -989,7 +1143,50 @@ setMethod("csSpecificity",signature(object="CuffData"),.specificity)
#res
}
-setMethod("makeRnk",signature(object="CuffData"),.makeRnk)
+.makeRnkTest<-function(object,x,y,filename,...){
+ samp<-samples(object)
+ #check to make sure x and y are in samples
+ if (!all(c(x,y) %in% samp)){
+ stop("One or more values of 'x' or 'y' are not valid sample names!")
+ }
+ query<-paste("SELECT gd.gene_id,g.gene_short_name,gd.test_stat FROM genes g LEFT JOIN geneExpDiffData gd ON g.gene_id=gd.gene_id WHERE ((gd.sample_1='",x,"' AND gd.sample_2='",y,"') OR (gd.sample_2='",x,"' AND gd.sample_1='",y,"')) ORDER BY gd.test_stat DESC",sep="")
+ #print(query)
+ res<-dbGetQuery(object at DB,query)
+ #Remove gene_id field
+ res<-res[,-1]
+ #Remove rows with "NA" for gene_short_name
+ res<-res[!is.na(res$gene_short_name),]
+ #Write to file
+ write.table(res,file=filename,sep="\t",quote=F,...,row.names=F,col.names=F)
+}
+
+.makeRnkTestStat<-function(object,x,y,filename,...){
+ samp<-samples(object)
+ #check to make sure x and y are in samples
+ if (!all(c(x,y) %in% samp)){
+ stop("One or more values of 'x' or 'y' are not valid sample names!")
+ }
+ query<-paste("SELECT gd.gene_id,g.gene_short_name,gd.sample_1, gd.sample_2,gd.test_stat FROM genes g LEFT JOIN geneExpDiffData gd ON g.gene_id=gd.gene_id WHERE ((gd.sample_1 ='",x,"' AND gd.sample_2='",y,"') OR (gd.sample_2='",x,"' AND gd.sample_1='",y,"')) GROUP BY gd.gene_id",sep="")
+ res<-dbGetQuery(object at DB,query)
+ if(unique(res$sample_2)==x){
+ res2<-res
+ res2$test_stat<--(res$test_stat)
+ res2$sample_1<-res$sample_2
+ res2$sample_2<-res$sample_1
+ res<-res2
+ }
+ #Order by test_stat
+ res<-res[order(res$test_stat,decreasing=F),]
+ #Remove gene_id field
+ res<-res[,c('gene_short_name','test_stat')]
+ #Remove rows with "NA" for gene_short_name
+ res<-res[!is.na(res$gene_short_name),]
+ #Write to file
+ write.table(res,file=filename,sep="\t",quote=F,...,row.names=F,col.names=F)
+ #head(res)
+}
+
+setMethod("makeRnk",signature(object="CuffData"),.makeRnkTestStat)
#############
diff --git a/R/methods-CuffFeature.R b/R/methods-CuffFeature.R
index 7929203..e88b830 100644
--- a/R/methods-CuffFeature.R
+++ b/R/methods-CuffFeature.R
@@ -92,6 +92,39 @@ setMethod("fpkm",signature="CuffFeature",.fpkm)
setMethod("fpkmMatrix",signature(object="CuffFeature"),.fpkmMatrix)
+#TODO: Replicate FPKM and Replicate FPKM Matrix for cuffFeature objects
+
+.repFpkmMatrix<-function(object,fullnames=FALSE,repIdList){
+ #Sample subsetting
+ if(!missing(repIdList)){
+ if (!all(repIdList %in% replicates(object))){
+ stop("Replicate does not exist!")
+ }else{
+ myReps<-repIdList
+ }
+ }else{
+ myReps<-replicates(object)
+ }
+ if(fullnames){
+ res<-repFpkm(object,features=TRUE)
+ res$tracking_id<-paste(res$gene_short_name,res[,1],sep="|")
+ }else{
+ res<-repFpkm(object)
+ colnames(res)[1]<-"tracking_id"
+ }
+ selectedRows<-c('tracking_id','rep_name','fpkm')
+ res<-res[,selectedRows]
+ res<-melt(res)
+ res<-dcast(res,tracking_id~rep_name)
+ res<-data.frame(res[,-1],row.names=res[,1])
+ if(!missing(repIdList)){
+ res<-res[,myReps]
+ }
+ res
+}
+
+setMethod("repFpkmMatrix",signature(object="CuffFeature"),.repFpkmMatrix)
+
.samples<-function(object){
res<-fpkm(object)$sample_name
res
@@ -99,6 +132,13 @@ setMethod("fpkmMatrix",signature(object="CuffFeature"),.fpkmMatrix)
setMethod("samples","CuffFeature",.samples)
+.replicates<-function(object){
+ res<-repFpkm(object)$rep_name
+ res
+}
+
+setMethod("replicates","CuffFeature",.replicates)
+
#setMethod("diff","CuffFeature",function(object){
# return(object at diff)
@@ -166,7 +206,7 @@ setMethod("count",signature(object="CuffFeature"),.count)
p<-ggplot(dat,aes(x=sample_name,y=fpkm,fill=sample_name))
- p <- p + geom_bar()
+ p <- p + geom_bar(stat="identity")
if(replicates){
p <- p + geom_point(aes(x=sample_name,y=fpkm),size=3,shape=18,colour="black",data=repDat)
@@ -215,7 +255,7 @@ setMethod("count",signature(object="CuffFeature"),.count)
setMethod("expressionBarplot",signature(object="CuffFeature"),.barplot)
-.expressionPlot<-function(object,logMode=FALSE,pseudocount=1.0, drawSummary=FALSE, sumFun=mean_cl_boot, showErrorbars=TRUE,showStatus=TRUE,replicates=FALSE,...){
+.expressionPlot<-function(object,logMode=FALSE,pseudocount=1.0, drawSummary=FALSE, sumFun=mean_cl_boot, showErrorbars=TRUE,showStatus=TRUE,replicates=FALSE,facet=TRUE,...){
#Coloring scheme for quant flags
quant_types<-c("OK","FAIL","LOWDATA","HIDATA","TOOSHORT")
quant_types<-factor(quant_types,levels=quant_types)
@@ -285,7 +325,9 @@ setMethod("expressionBarplot",signature(object="CuffFeature"),.barplot)
#Recolor quant flags
#for some reason this doesn't work (ggplot2 problem)
#p<- p+ scale_colour_manual(name='quant_status',values=quant_colors)
- p<-p+facet_wrap('tracking_id')
+ if(facet){
+ p<-p+facet_wrap('tracking_id')
+ }
p
}
diff --git a/R/methods-CuffFeatureSet.R b/R/methods-CuffFeatureSet.R
index d8dbc17..c83223c 100644
--- a/R/methods-CuffFeatureSet.R
+++ b/R/methods-CuffFeatureSet.R
@@ -251,6 +251,20 @@ setMethod("repCountMatrix",signature(object="CuffFeatureSet"),.repCountMatrix)
setMethod("diffData",signature(object="CuffFeatureSet"),.diffData)
+.diffTable<-function(object,logCutoffValue=99999){
+ measureVars<-c('status','value_1','value_2','log2_fold_change','test_stat','p_value','q_value','significant')
+ all.diff<-diffData(object,features=TRUE)
+ all.diff$log2_fold_change[all.diff$log2_fold_change>=logCutoffValue]<-Inf
+ all.diff$log2_fold_change[all.diff$log2_fold_change<=-logCutoffValue]<--Inf
+ all.diff.melt<-melt(all.diff,measure.vars=measureVars)
+ #all.diff.melt<-all.diff.melt[!grepl("^value_",all.diff.melt$variable),]
+ all.diff.cast<-dcast(all.diff.melt,formula=...~sample_2+sample_1+variable)
+ all.diff.cast
+}
+
+setMethod("diffTable",signature(object="CuffFeatureSet"),.diffTable)
+
+
.count<-function(object,features=FALSE){
if (features){
return (merge(object at annotation,object at count,by=1))
@@ -310,6 +324,22 @@ setMethod("annotation","CuffFeatureSet",function(object){
m = log10(m+pseudocount)
}
+
+ if(is.function(rescaling))
+ {
+ m=rescaling(m)
+ } else {
+ if(rescaling=='column'){
+ m=scale(m, center=T)
+ m[is.nan(m)] = 0
+ }
+ if(rescaling=='row'){
+ m=t(scale(t(m),center=T))
+ m[is.nan(m)] = 0
+ }
+ }
+
+
## I have supplied the default cluster and euclidean distance (JSdist) - and chose to cluster after scaling
## if you want a different distance/cluster method-- or to cluster and then scale
## then you can supply a custom function
@@ -327,20 +357,7 @@ setMethod("annotation","CuffFeatureSet",function(object){
## this is just reshaping into a ggplot format matrix and making a ggplot layer
- if(is.function(rescaling))
- {
- m=rescaling(m)
- } else {
- if(rescaling=='column'){
- m=scale(m, center=T)
- m[is.nan(m)] = 0
- }
- if(rescaling=='row'){
- m=t(scale(t(m),center=T))
- m[is.nan(m)] = 0
- }
- }
-
+
rows=dim(m)[1]
cols=dim(m)[2]
@@ -430,6 +447,135 @@ setMethod("annotation","CuffFeatureSet",function(object){
setMethod("csHeatmap",signature("CuffFeatureSet"),.ggheat)
+.fcheatmap<-function(object, control_condition, replicate_num=NULL, clustering='none', labCol=T, labRow=T, logMode=F, pseudocount=1.0,
+ border=FALSE, heatscale=c(low='steelblue',mid='white',high='tomato'), heatMidpoint=0,fullnames=T,replicates=FALSE,method='none',heatRange=3, ...) {
+ ## the function can be be viewed as a two step process
+ ## 1. using the rehape package and other funcs the data is clustered, scaled, and reshaped
+ ## using simple options or by a user supplied function
+ ## 2. with the now resahped data the plot, the chosen labels and plot style are built
+
+ if(replicates){
+ if (is.null(replicate_num)){
+ print ("Error: if replicates == TRUE, you must specify both a control condition and a replicate number")
+ return()
+ }
+ m=repFpkmMatrix(object,fullnames=fullnames)
+ selected_rep <- paste(control_condition,replicate_num,sep="_")
+
+ #remove genes with no expression in any condition
+ m=m[!apply(m,1,sum)==0,]
+ m=log2(m/m[,selected_rep])
+ m=m[,names(m) != selected_rep]
+ }else{
+ m=fpkmMatrix(object,fullnames=fullnames)
+ #remove genes with no expression in any condition
+ m=m[!apply(m,1,sum)==0,]
+ m=log2(m/m[,control_condition])
+ m=m[,names(m) != control_condition]
+ }
+
+ m_vals <- unlist(as.list(m))[is.finite(unlist(as.list(m)))]
+ m_max <- max(m_vals)
+ m_min <- max(m_vals)
+ range_lim <- max(c(abs(m_max), abs(m_min)))
+ range_lim <- min(c(heatRange, range_lim))
+ m_max <- range_lim
+ m_min <- -range_lim
+ m[m < m_min] <- m_min
+ m[m > m_max] <- m_max
+ m[is.na(m)] <- 0
+
+ ## I have supplied the default cluster and euclidean distance (JSdist) - and chose to cluster after scaling
+ ## if you want a different distance/cluster method-- or to cluster and then scale
+ ## then you can supply a custom function
+
+ if(!is.function(method)){
+ method = dist
+ }
+
+ if(clustering=='row')
+ m=m[hclust(method(m))$order, ]
+ if(clustering=='column')
+ m=m[,hclust(method(t(m)))$order]
+ if(clustering=='both')
+ m=m[hclust(method(m))$order ,hclust(method(t(m)))$order]
+
+ rows=dim(m)[1]
+ cols=dim(m)[2]
+
+ melt.m=cbind(rowInd=rep(1:rows, times=cols), colInd=rep(1:cols, each=rows), melt(m))
+
+ g=ggplot(data=melt.m)
+
+ ## add the heat tiles with or without a white border for clarity
+
+ if(border==TRUE)
+ g2=g+geom_rect(aes(xmin=colInd-1,xmax=colInd,ymin=rowInd-1,ymax=rowInd, fill=value),colour='grey')
+ if(border==FALSE)
+ g2=g+geom_rect(aes(xmin=colInd-1,xmax=colInd,ymin=rowInd-1,ymax=rowInd, fill=value))
+
+ ## add axis labels either supplied or from the colnames rownames of the matrix
+
+ if(labCol==T)
+ {
+ g2=g2+scale_x_continuous(breaks=(1:cols)-0.5, labels=colnames(m))
+ }
+ if(labCol==F)
+ {
+ g2=g2+scale_x_continuous(breaks=(1:cols)-0.5, labels=rep('',cols))
+ }
+
+
+ if(labRow==T)
+ {
+ g2=g2+scale_y_continuous(breaks=(1:rows)-0.5, labels=rownames(m))
+ }
+ if(labRow==F)
+ {
+ g2=g2+scale_y_continuous(breaks=(1:rows)-0.5, labels=rep('',rows))
+ }
+
+ # Get rid of the ticks, they get way too dense with lots of rows
+ g2 <- g2 + theme(axis.ticks = element_blank())
+
+ ## get rid of grey panel background and gridlines
+
+ g2=g2+theme(panel.grid.minor=element_line(colour=NA), panel.grid.major=element_line(colour=NA),
+ panel.background=element_rect(fill=NA, colour=NA))
+
+ ##adjust x-axis labels
+ g2=g2+theme(axis.text.x=element_text(angle=-90, hjust=0))
+
+ #write(paste(c("Length of heatscale is :", length(heatscale))), stderr())
+ if(replicates){
+ legendTitle <- bquote(paste(log[2], frac("FPKM",.(selected_rep),sep="")))
+ }else{
+ legendTitle <- bquote(paste(log[2], frac("FPKM",.(control_condition),sep="")))
+ }
+ #legendTitle <- paste(expression(plain(log)[10])," FPKM + ",pseudocount,sep="")
+
+ if (length(heatscale) == 2){
+ g2 <- g2 + scale_fill_gradient(low=heatscale[1], high=heatscale[2], name=legendTitle)
+ } else if (length(heatscale) == 3) {
+ if (is.null(heatMidpoint))
+ {
+ heatMidpoint = (max(m) + min(m)) / 2.0
+ #write(heatMidpoint, stderr())
+ }
+
+ g2 <- g2 + scale_fill_gradient2(low=heatscale[1], mid=heatscale[2], high=heatscale[3], midpoint=heatMidpoint, name=legendTitle)
+ }
+
+ #g2<-g2+scale_x_discrete("",breaks=tracking_ids,labels=gene_short_names)
+
+
+ ## finally add the fill colour ramp of your choice (default is blue to red)-- and return
+ return (g2)
+
+}
+
+setMethod("csFoldChangeHeatmap",signature("CuffFeatureSet"),.fcheatmap)
+
# Distance Heatmaps
.distheat<-function(object, replicates=F, samples.not.genes=T, logMode=T, pseudocount=1.0, heatscale=c(low='lightyellow',mid='orange',high='darkred'), heatMidpoint=NULL, ...) {
# get expression from a sample or gene perspective
@@ -493,7 +639,8 @@ setMethod("csDistHeat", signature("CuffFeatureSet"), .distheat)
#Scatterplot
.scatter<-function(object,x,y,logMode=TRUE,pseudocount=0.0,labels, smooth=FALSE,colorByStatus=FALSE,...){
- dat<-fpkmMatrix(object)
+ dat<-fpkmMatrix(object,fullnames=T)
+
samp<-samples(object)
#check to make sure x and y are in samples
@@ -508,6 +655,16 @@ setMethod("csDistHeat", signature("CuffFeatureSet"), .distheat)
}
}
+ #Attach tracking_id and gene_short_name
+ if(!missing(labels)){
+ require(stringr)
+ tracking<-str_split_fixed(rownames(dat),"\\|",2)
+ dat$gene_short_name<-tracking[,1]
+ dat$tracking_id<-tracking[,2]
+
+ labeled.dat<-dat[dat$gene_short_name %in% labels,]
+ }
+
#make plot object
p<-ggplot(dat)
p<- p + aes_string(x=x,y=y)
@@ -519,15 +676,11 @@ setMethod("csDistHeat", signature("CuffFeatureSet"), .distheat)
}
#Add highlights from labels
-# if(!missing(labels)){
-# labelIdx<-fData(object)$gene_short_name %in% labels
-# labelfp<-fp[labelIdx,]
-# labelfp$gene_short_name<-fData(object)$gene_short_name[labelIdx]
-# #print(head(labelfp))
-# p <- p + geom_point(data=labelfp,size=1.2,color="red")
-# p <- p + geom_text(data=labelfp,aes(label=gene_short_name),color="red",hjust=0,vjust=0,angle=45,size=2)
-# }
-#
+ if(!missing(labels)){
+ p <- p + geom_point(data=labeled.dat,aes_string(x=x,y=y),size=1.3,color="red")
+ p <- p + geom_text(data=labeled.dat,aes_string(x=x,y=y,label='gene_short_name'),color="red",hjust=0,vjust=0,angle=0,size=4)
+ }
+
#logMode
if(logMode){
p <- p + scale_y_log10() + scale_x_log10()
@@ -547,6 +700,8 @@ setMethod("csDistHeat", signature("CuffFeatureSet"), .distheat)
p
}
+#.scatter(sigGenes,'P7_lincBrn1b_KO_Brain','P7_WT_Brain',labels=c('Arc'))
+
setMethod("csScatter",signature(object="CuffFeatureSet"), .scatter)
#Volcano plot
@@ -777,14 +932,18 @@ setMethod("expressionPlot",signature(object="CuffFeatureSet"),.expressionPlot)
# c
#}
-.cluster<-function(object, k, logMode=T, pseudocount=1,...){
+.cluster<-function(object, k, logMode=T, method='none', pseudocount=1,...){
require(cluster)
m<-as.data.frame(fpkmMatrix(object))
m<-m[rowSums(m)>0,]
if(logMode){
m<-log10(m+pseudocount)
}
- n<-JSdist(makeprobs(t(m)))
+
+ if(!is.function(method)){
+ method = function(mat){JSdist(makeprobs(t(m)))}
+ }
+ n<-method(m)
clusters<-pam(n,k, ...)
#clsuters<-pamk(n,krange=2:20)
class(clusters)<-"list"
@@ -794,33 +953,59 @@ setMethod("expressionPlot",signature(object="CuffFeatureSet"),.expressionPlot)
setMethod("csCluster",signature(object="CuffFeatureSet"),.cluster)
-.ratioCluster<-function(object,k,ratioTo=NULL,...){
+.ratioCluster<-function(object,k,ratioTo=NULL,pseudocount=0.0001,...){
require(cluster)
- m<-as.data.frame(fpkmMatrix(object))
+ m<-as.data.frame(fpkmMatrix(object)+pseudocount)
#TODO: ensure that ratioTo is in colnames(fpkmMatrix(object))
m.ratio<-m/m[[ratioTo]]
m.log.ratio<-log2(m.ratio)
n<-dist(m.log.ratio)
clusters<-pam(n,k,...)
class(clusters)<-"list"
- clusters$ratio<-m
+ clusters$fpkm<-m
clusters
}
-csClusterPlot<-function(clustering,pseudocount=1.0,drawSummary=TRUE, sumFun=mean_cl_boot){
+csClusterPlot<-function(clustering,pseudocount=1.0,logMode=FALSE,drawSummary=TRUE, sumFun=mean_cl_boot){
m<-clustering$fpkm+pseudocount
m$ids<-rownames(clustering$fpkm)
m$cluster<-factor(clustering$clustering)
m.melt<-melt(m,id.vars=c("ids","cluster"))
c<-ggplot(m.melt)
- c<-c+geom_line(aes(x=variable,y=value,color=cluster,group=ids)) + facet_wrap('cluster',scales='free')+scale_y_log10()
+ c<-c+geom_line(aes(x=variable,y=value,color=cluster,group=ids)) + theme_bw() + facet_wrap('cluster',scales='free_y')
if(drawSummary){
c <- c + stat_summary(aes(x=variable,y=value,group=1),fun.data=sumFun,color="black",fill="black",alpha=0.2,size=1.1,geom="smooth")
}
- c<-c + scale_color_hue(l=50,h.start=200)
+ if(logMode){
+ c<-c + scale_y_log10()
+ }
+ c<-c + scale_color_hue(l=50,h.start=200) + theme(axis.text.x=element_text(angle=-90,hjust=0))
c
}
+.findK<-function(object, k.range=c(2:20), logMode=T, pseudocount=1,...){
+ require(cluster)
+ m<-as.data.frame(fpkmMatrix(object))
+ m<-m[rowSums(m)>0,]
+ if(logMode){
+ m<-log10(m+pseudocount)
+ }
+ n<-JSdist(makeprobs(t(m)))
+ myWidths<-c()
+ for (k in k.range){
+ print(k)
+ myWidths<-c(myWidths,pam(n,k,...)$silinfo$avg.width)
+ }
+ plot(k.range,myWidths)
+}
+
+
+
+
+######################
+#
+######################
+
.dendro<-function(object,logMode=T,pseudocount=1,replicates=FALSE){
if(replicates){
fpkmMat<-repFpkmMatrix(object)
@@ -869,6 +1054,30 @@ setMethod("csDendro",signature(object="CuffFeatureSet"),.dendro)
# c
#}
+######################
+# Exploratory Analysis
+######################
+.nmf<-function(object,k,logMode=T,pseudocount=1,maxiter=1000,replicates=FALSE,fullnames=FALSE){
+ require(NMFN)
+ if(missing(k)) stop("Please provide a rank value for factorization (arg=k)")
+
+ if(replicates){
+ m=repFpkmMatrix(object,fullnames=fullnames)
+ }else{
+ m=fpkmMatrix(object,fullnames=fullnames)
+ }
+
+ if(logMode)
+ {
+ m = log10(m+pseudocount)
+ }
+
+ myNMF<-nnmf(m,k=k,maxiter=maxiter)
+ return (myNMF)
+}
+
+setMethod("csNMF",signature(object="CuffFeatureSet"),.nmf)
+
.density<-function(object, logMode = TRUE, pseudocount=0.0, labels, features=FALSE, replicates=FALSE,...){
if(replicates){
dat<-repFpkm(object,features=features)
@@ -896,6 +1105,65 @@ setMethod("csDendro",signature(object="CuffFeatureSet"),.dendro)
setMethod("csDensity",signature(object="CuffFeatureSet"),.density)
+#################
+# Dimensionality Reduction (use sparingly for Gene Sets)
+#################
+
+.MDSplot<-function(object,replicates=FALSE,logMode=TRUE,pseudocount=1.0){
+ if(replicates){
+ dat<-repFpkmMatrix(object)
+ #repData<-sapply(replicates(object),function(x){strsplit(x,"_")[[1]][1]}) This is to color by condition and not replicate...
+ }else{
+ dat<-fpkmMatrix(object)
+ }
+
+ if(logMode){
+ dat<-log10(dat+pseudocount)
+ }
+
+ d<-JSdist(makeprobs(dat))
+ fit <- cmdscale(d,eig=TRUE, k=2)
+ res<-data.frame(names=rownames(fit$points),M1=fit$points[,1],M2=fit$points[,2])
+ p <- ggplot(res)
+ p <- p + geom_point(aes(x=M1,y=M2,color=names)) + geom_text(aes(x=M1,y=M2,label=names,color=names)) + theme_bw()
+ p
+}
+
+setMethod("MDSplot",signature(object="CuffFeatureSet"),.MDSplot)
+
+.PCAplot<-function(object,x="PC1", y="PC2",replicates=FALSE,pseudocount=1.0,scale=TRUE,showPoints=TRUE,...){
+ if(replicates){
+ fpkms<-repFpkmMatrix(object)
+ }else{
+ fpkms<-fpkmMatrix(object)
+ }
+ fpkms<-log10(fpkms+pseudocount)
+ PC<-prcomp(fpkms,scale=scale,...)
+ dat <- data.frame(obsnames=row.names(PC$x), PC$x)
+ #dat$shoutout<-""
+ #dat$shoutout[matchpt(PC$rotation,PC$x)$index]<-rownames(pca$x[matchpt(pca$rotation,pca$x)$index,])
+ plot <- ggplot(dat, aes_string(x=x, y=y))
+ if(showPoints){
+ plot<- plot + geom_point(alpha=.4, size=0.8, aes(label=obsnames))
+ }
+ plot <- plot + geom_hline(aes(0), size=.2) + geom_vline(aes(0), size=.2) #+ geom_text(aes(label=shoutout),size=2,color="red")
+ datapc <- data.frame(varnames=rownames(PC$rotation), PC$rotation)
+ mult <- min(
+ (max(dat[,y]) - min(dat[,y])/(max(datapc[,y])-min(datapc[,y]))),
+ (max(dat[,x]) - min(dat[,x])/(max(datapc[,x])-min(datapc[,x])))
+ )
+ datapc <- transform(datapc,
+ v1 = .7 * mult * (get(x)),
+ v2 = .7 * mult * (get(y))
+ )
+ plot <- plot +
+ #coord_equal() +
+ geom_text(data=datapc, aes(x=v1, y=v2, label=varnames,color=varnames), vjust=1)
+ plot <- plot + geom_segment(data=datapc, aes(x=0, y=0, xend=v1, yend=v2,color=varnames), arrow=arrow(length=unit(0.2,"cm")), alpha=0.75) + theme_bw()
+ plot
+}
+
+setMethod('PCAplot',signature(object="CuffFeatureSet"),.PCAplot)
#################
#Misc #
diff --git a/R/methods-CuffGene.R b/R/methods-CuffGene.R
index 5fc828a..9f76284 100644
--- a/R/methods-CuffGene.R
+++ b/R/methods-CuffGene.R
@@ -131,22 +131,99 @@ setMethod("makeGeneRegionTrack",signature(object="CuffGene"),.makeGeneRegionTrac
setMethod("genePlot",signature(object="CuffGene"),.plot)
+
+
+#################
+#Feature Plotting
+#################
+.ideogram<-function(object){
+ myStart<-min(object at features$start)
+ myEnd<-max(object at features$end)
+ mychr<-unique(object at features$seqnames)
+ p<-plotIdeogram(genome=object at genome,subchr=mychr,zoom.region=c(myStart,myEnd))
+ p
+
+}
+
+.plot2<-function(object,...){
+ #Ideogram
+ ideoTrack<-.ideogram(object)
+
+ #Expression levels
+ expressionTrack<-expressionPlot(isoforms(object),facet=T,...)+theme_bw() + theme(legend.position='none')
+ hasAxis(expressionTrack)<-TRUE
+
+ #Transcript Models
+ #modelTrack<-autoplot(.asGRangesList(object),aes(fill=transcript,group=transcript),gap.geom="arrow") + theme_bw() + scale_fill_hue(l=50,h.start=200) + scale_color_hue(l=50,h.start=200)
+ modelTrack<-ggplot(as.GRangesList(object),)
+
+ hasAxis(modelTrack)<-TRUE
+
+ #Plot it all...
+ tracks(ideoTrack,modelTrack,expressionTrack,heights=c(1,3,3),fixed=c(TRUE,TRUE,FALSE),main=unique(object at annotation$gene_short_name))
+
+}
+
+.pie<-function(object,level="isoforms",pseudocount=0.0001,...){
+ dat<-fpkm(slot(object,level))
+ colnames(dat)[1]<-'tracking_id'
+ #dat$fpkm<-dat$fpkm+pseudocount
+ #print(dat)
+ p<-ggplot(dat,aes(x="",y=fpkm,fill=tracking_id))
+
+ p<- p + geom_bar(stat="identity",position="fill",line="black")
+
+ p<- p + coord_polar(theta='y')
+
+ p<-p + scale_fill_hue(l=50,h.start=200) + scale_color_hue(l=50,h.start=200) + theme_bw()
+
+ p<-p + facet_wrap('sample_name') + theme(axis.text.x = element_blank(),aspect.ratio=1)
+
+ p
+}
+
+setMethod("csPie",signature(object="CuffGene"),.pie)
+
#################
#Coersion methods
#################
#As GRanges
-.asGRanges<-function(object){
+.as.GRanges<-function(from){
#featCols<-c('seqnames','start','end','source','gene_id','exon_number','isoform_id','isoform_id','exon_number','strand')
- feats<-object at features
+ feats<-from at features
#newColnames<-c('seqnames','start','end','feature','gene','exon','transcript','symbol','rank','strand')
+ #ExpressionValues (transcript)
+ fpkm<-cbind(isoform_id=rownames(fpkmMatrix(isoforms(from))),fpkmMatrix(isoforms(from)))
+ feats<-merge(feats,fpkm)
colnames(feats)[colnames(feats)=='isoform_id']<-'transcript'
colnames(feats)[colnames(feats)=='gene_id']<-'gene'
colnames(feats)[colnames(feats)=='exon_number']<-'exon'
colnames(feats)[colnames(feats)=='source']<-'feature'
feats$symbol<-feats$transcript
- corCols<-c('seqnames','start','end','strand')
+ corCols<-c('seqnames','start','end','strand','width')
myGR<-GRanges(Rle(feats$seqnames),ranges=IRanges(feats$start,end=feats$end),strand=Rle(feats$strand),elementMetadata=feats[,!colnames(feats) %in% corCols])
+ colnames(elementMetadata(myGR))<-colnames(feats[,!colnames(feats) %in% corCols])
myGR
}
-#As GRangesList
\ No newline at end of file
+setAs("CuffGene","GRanges",.as.GRanges)
+
+#As GRangesList
+.as.GRangesList<-function(object,f="transcript"){
+ gr<-as(object,"GRanges")
+ grl<-split(gr,f)
+ grl
+}
+setMethod("as.GRangesList",signature(object="CuffGene"),.as.GRangesList)
+
+
+#######################
+# ggbio integration
+#######################
+#setMethod("ggplot", "CuffGene", function(data, ...){
+# df <- mold(as.GRangesList(data))
+# g <- ggplot(df, ...)
+# g$.data <- as.GRangesList(data)
+# g <- ggbio(g)
+# g
+# })
\ No newline at end of file
diff --git a/R/methods-CuffSet.R b/R/methods-CuffSet.R
index 07447e0..6a997bb 100644
--- a/R/methods-CuffSet.R
+++ b/R/methods-CuffSet.R
@@ -11,7 +11,7 @@
setMethod("initialize","CuffSet",
function(.Object,
DB,
- runInfo=data.frame(),
+ #runInfo=data.frame(),
phenoData=data.frame(),
conditions=data.frame(),
genes,
@@ -24,7 +24,7 @@ setMethod("initialize","CuffSet",
...){
.Object<-callNextMethod(.Object,
DB = DB,
- runInfo=runInfo,
+ #runInfo=runInfo,
phenoData=phenoData,
conditions = conditions,
genes = genes,
@@ -110,6 +110,13 @@ setMethod("DB","CuffSet",function(object){
setMethod("runInfo","CuffSet",.runInfo)
+.varModel<-function(object){
+ varModelQuery<-"SELECT * from varModel"
+ dbGetQuery(object at DB,varModelQuery)
+}
+
+setMethod("varModel","CuffSet",.varModel)
+
#setMethod("phenoData","CuffSet",function(object){
# return(object at phenoData)
# })
@@ -324,8 +331,10 @@ setMethod("getGene",signature(object="CuffSet"),.getGene)
}
sampleString<-substr(sampleString,1,nchar(sampleString)-1)
sampleString<-paste(sampleString,")",sep="")
+
- #ID Search String (SQL)
+ #idQuery<-paste('SELECT DISTINCT x.gene_id from genes x WHERE (x.gene_id IN ',origIdString,' OR x.gene_short_name IN ',origIdString,')',sep="")
+
idString<-'('
for (i in geneIdList){
idString<-paste(idString,"'",i,"',",sep="")
@@ -333,12 +342,14 @@ setMethod("getGene",signature(object="CuffSet"),.getGene)
idString<-substr(idString,1,nchar(idString)-1)
idString<-paste(idString,")",sep="")
- whereStringGene<-paste('WHERE (x.gene_id IN ',idString,' OR x.gene_short_name IN ',idString,')',sep="")
- whereStringGeneFPKM<-paste('WHERE (x.gene_id IN ',idString,' OR x.gene_short_name IN ',idString,')',sep="")
- whereStringGeneDiff<-paste('WHERE (x.gene_id IN ',idString,' OR x.gene_short_name IN ',idString,')',sep="")
- whereString<-paste('WHERE (x.gene_id IN ',idString,' OR g.gene_short_name IN ',idString,')',sep="")
- whereStringFPKM<-paste('WHERE (x.gene_id IN ',idString,' OR g.gene_short_name IN ',idString,')',sep="")
- whereStringDiff<-paste('WHERE (x.gene_id IN ',idString,' OR g.gene_short_name IN ',idString,')',sep="")
+ #write(idString,stderr())
+
+ whereStringGene<-paste('WHERE x.gene_id IN ',idString,sep="")
+ whereStringGeneFPKM<-paste('WHERE x.gene_id IN ',idString,sep="")
+ whereStringGeneDiff<-paste('WHERE x.gene_id IN ',idString,sep="")
+ whereString<-paste('WHERE x.gene_id IN ',idString,sep="")
+ whereStringFPKM<-paste('WHERE x.gene_id IN ',idString,sep="")
+ whereStringDiff<-paste('WHERE x.gene_id IN ',idString,sep="")
if(!is.null(sampleIdList)){
whereStringGene<-whereStringGene
@@ -350,9 +361,6 @@ setMethod("getGene",signature(object="CuffSet"),.getGene)
}
- #dbQueries
- idQuery<-paste("SELECT DISTINCT gene_id from genes x ",whereStringGene,sep="")
-
#geneAnnotationQuery<-paste("SELECT * from genes x LEFT JOIN geneFeatures xf ON x.gene_id=xf.gene_id ", whereStringGene,sep="")
geneAnnotationQuery<-paste("SELECT * from genes x LEFT JOIN geneFeatures xf USING (gene_id) ", whereStringGene,sep="")
geneFPKMQuery<-paste("SELECT y.* from genes x JOIN geneData y ON x.gene_id=y.gene_id ", whereStringGeneFPKM,sep="")
@@ -474,7 +482,7 @@ setMethod("getGene",signature(object="CuffSet"),.getGene)
res<-new("CuffGeneSet",
#TODO: Fix ids so that it only displays those genes in CuffGeneSet
- ids=as.character(dbGetQuery(object at DB,idQuery)),
+ ids=geneIdList,
annotation=genes.annot,
fpkm=genes.fpkm,
diff=genes.diff,
@@ -707,7 +715,7 @@ setMethod("getFeatures",signature(object="CuffSet"),.getFeatures)
#
#}
-.getSig<-function(object,x,y,alpha=0.05,level='genes'){
+.getSig<-function(object,x,y,alpha=0.05,level='genes',method="BH",useCuffMTC=FALSE){
mySamp<-samples(slot(object,level))
if(level %in% c('promoters','splicing','relCDS')){
@@ -723,12 +731,14 @@ setMethod("getFeatures",signature(object="CuffSet"),.getFeatures)
stop("One or more values of 'x' or 'y' are not valid sample names!")
}
}
-
queryString<-paste("(",paste(mySamp,collapse="','",sep=""),")",sep="'")
- sql<-paste("SELECT ",slot(object,level)@idField,",p_value from ", diffTable," WHERE sample_1 IN ",queryString," AND sample_2 IN ",queryString, " AND STATUS='OK'",sep="")
+ sql<-paste("SELECT ",slot(object,level)@idField,",p_value,q_value from ", diffTable," WHERE sample_1 IN ",queryString," AND sample_2 IN ",queryString, " AND STATUS='OK'",sep="")
#print(sql)
sig<-dbGetQuery(object at DB,sql)
- sig$q_value<-p.adjust(sig$p_value,method="BH")
+ if(!missing(x) && !missing(y) && !useCuffMTC){
+ sig$q_value<-p.adjust(sig$p_value,method=method)
+ #print(sig[order(sig$q_value,decreasing=T),])
+ }
sig<-sig[sig$q_value<=alpha,]
sigGenes<-unique(sig[[slot(object,level)@idField]])
sigGenes
@@ -768,6 +778,10 @@ setMethod("getSig",signature(object="CuffSet"),.getSig)
setMethod("getSigTable",signature(object="CuffSet"),.getSigTable)
+####################
+#QC Visualizations on entire dataset
+#####################
+
.sigMatrix<-function(object,alpha=0.05,level='genes',orderByDist=FALSE){
if(level %in% c('promoters','splicing','relCDS')){
diffTable<-slot(object,level)@table
@@ -809,6 +823,17 @@ setMethod("getSigTable",signature(object="CuffSet"),.getSigTable)
setMethod("sigMatrix",signature(object="CuffSet"),.sigMatrix)
+#dispersionPlot on cuffSet objects actually draws from the varModel table and is the preferred way to visualize the dispersion and model fitting from cuffdiff 2.1 or greater
+.dispersionPlot<-function(object){
+ dat<-varModel(object)
+ p<-ggplot(dat)
+ p<-p + geom_point(aes(x=compatible_count_mean,y=compatible_count_var,color=condition),alpha=0.3,size=0.8) + geom_line(aes(x=compatible_count_mean,y=fitted_var),lwd=0.5,color="black") + facet_wrap('condition') +scale_y_log10() + scale_x_log10() + theme_bw() + guides(color=FALSE)
+ p
+}
+
+setMethod("dispersionPlot",signature(object="CuffSet"),.dispersionPlot)
+
+
#Find similar genes
.findSimilar<-function(object,x,n,distThresh,returnGeneSet=TRUE,...){
#x can be either a gene_id, gene_short_name or a vector of FPKM values (fake gene expression profile)
diff --git a/R/tools.R b/R/tools.R
index 958fd6b..a6e41e3 100644
--- a/R/tools.R
+++ b/R/tools.R
@@ -3,7 +3,7 @@
# Author: lgoff
###############################################################################
-JSdist<-function(mat){
+JSdist<-function(mat,...){
res<-matrix(0,ncol=dim(mat)[2],nrow=dim(mat)[2])
# col_js <- matrix(0,ncol=dim(mat)[2],nrow=1)
@@ -23,7 +23,9 @@ JSdist<-function(mat){
res[j,i] = sqrt(JSdiv)
}
}
- as.dist(res)
+ res<-as.dist(res,...)
+ attr(res,"method")<-"JSdist"
+ res
}
JSdistVec<-function(p,q){
@@ -32,6 +34,12 @@ JSdistVec<-function(p,q){
JSdist
}
+JSdivVec<-function(p,q){
+ JSdiv<-shannon.entropy((p+q)/2)-(shannon.entropy(p)+shannon.entropy(q))*0.5
+ #JSdist<-sqrt(JSdiv)
+ JSdiv
+}
+
JSdistFromP<-function(mat,q){
#row_js<-apply(mat,MARGIN=1,shannon.entropy)
res<-apply(mat,MARGIN=1,function(p) {
@@ -51,7 +59,7 @@ shannon.entropy <- function(p) {
if (min(p) < 0 || sum(p) <=0)
return(Inf)
p.norm<-p[p>0]/sum(p)
- -sum( log10(p.norm)*p.norm)
+ -sum( log2(p.norm)*p.norm)
}
makeprobs<-function(a){
@@ -95,24 +103,51 @@ makeprobs<-function(a){
p
}
-#.volcanoMatrix <- function(data){
-# part1<-data[,c('gene_id','sample_1','sample_2','value_1','value_2','test_stat','p_value','q_value')]
-# part2<-data.frame(gene_id=part1$gene_id,sample_1=part1$sample_2,sample_2=part1$sample_1,value_1=part1$value_2,value_2=part1$value_1,test_stat=-part1$test_stat,p_value=part1$p_value,q_value=part1$q_value)
-# data<-rbind(part1,part2)
-# myLevels<-union(data$sample_1,data$sample_2)
-# data$sample_1<-factor(data$sample_1,levels=myLevels)
-# data$sample_2<-factor(data$sample_2,levels=myLevels)
-# data$log2_fold_change<-log2(data$value_2/data$value_1)
-# filler<-data.frame(sample_1=factor(myLevels,levels=myLevels),sample_2=factor(myLevels,levels=myLevels),label=as.character(myLevels))
-# filler$label<-as.character(filler$label)
-# mapping <- defaults(mapping, aes_string(x = "log2_fold_change", y = "-log10(p_value)"))
-# class(mapping) <- "uneval"
-#
-# p <-ggplot(data) + geom_point(mapping,na.rm=TRUE,size=0.8) + geom_text(aes(x=0,y=15,label=label),data=filler) + facet_grid(sample_1~sample_2)
-#
-# p
-#
-#}
+.hclustToJSON<-function(d,...){
+ require(rjson)
+ d$call<-NULL
+ res<-toJSON(d,...)
+ res
+}
+
+
+.dfToJSONarray <- function(dtf){
+ clnms <- colnames(dtf)
+
+ name.value <- function(i){
+ quote <- '';
+ if(class(dtf[, i])!='numeric'){
+ quote <- '"';
+ }
+
+ paste('"', i, '" : ', quote, dtf[,i], quote, sep='')
+ }
+
+ objs <- apply(sapply(clnms, name.value), 1, function(x){paste(x, collapse=', ')})
+ objs <- paste('{', objs, '}')
+
+ res <- paste('[', paste(objs, collapse=', '), ']')
+
+ return(res)
+}
+
+.specificity<-function(mat,logMode=T,pseudocount=1,relative=FALSE,...){
+ if(logMode){
+ mat<-log10(mat+pseudocount)
+ }
+ mat<-t(makeprobs(t(mat)))
+ d<-diag(ncol(mat))
+ res<-apply(d,MARGIN=1,function(q){
+ JSdistFromP(mat,q)
+ })
+ colnames(res)<-paste(colnames(mat),"_spec",sep="")
+
+ if(relative){
+ res<-res/max(res)
+ }
+ 1-res
+}
+
#multiplot <- function(..., plotlist=NULL, cols) {
# require(grid)
diff --git a/build/vignette.rds b/build/vignette.rds
index 74f8238..0da12ed 100644
Binary files a/build/vignette.rds and b/build/vignette.rds differ
diff --git a/data/sampleData.rda b/data/sampleData.rda
new file mode 100644
index 0000000..fdac3bd
Binary files /dev/null and b/data/sampleData.rda differ
diff --git a/inst/NEWS b/inst/NEWS
index cd777ad..1377da1 100644
--- a/inst/NEWS
+++ b/inst/NEWS
@@ -1,3 +1,34 @@
+v2.7.1
+ Bugfixes:
+ - Fixed 'fullnames' argument to cuffData::*Matrix() methods so that it does what it's supposed to do.
+ - Added 'showPool' argument to fpkmSCVPlot. When TRUE, empirical mean and standard deviation are determined across all conditions as opposed to cross-replicate. This is set to TRUE anytime you have n<2 replicates per condition.
+ - Added stat="identity" to expressionBarplot to comply with ggplot 0.9.3 enforcement.
+ - 'labels' argument to csScatter is now working as it's supposed to. You can pass a vector of 'gene_short_name' identifiers to labels and these will be specifically called out in red text on scatterplot.
+ - Added repFpkmMatrix() and replicates() methods to CuffFeature objects.
+ - Removed unnecessary Joins to optimize retrieval speed for several key queries.
+ - Fixed bug in csVolcano matrix that forced ylimits to be c(0,15)
+
+ New Features:
+ - Added csNMF() method for CuffData and CuffFeatureSet objects to perform non-negative matrix factorization. As of now, it's merely a wrapper around the default settings for NMFN::nnmf(), but hope to expand in the future.
+ * Does not adjust sparsity of matrices after output, must be done by user as needed.
+ - Added csPie() method for CuffGene objects. Allows for visualization of relative isoform, CDS, and promoter usage proportions as a pie chart by condition (or optionally as stacked bar charts by adding + coord_cartesian() ).
+ - Added 'method' argument to csCluster and csHeatmap to allow custom distance functions for clustering. Default = "none" = JSdist(). You can now provide a function that returns a 'dist' object on rows of a matrix.
+ - Added varModel.info tracking for compatibility with cuffdiff >=2.1. Will now find varModel.info file if exists, and incorporate into database.
+ - dispersionPlot() method added for CuffSet object. This now appropriately draws from varModel.info and is the preferred visualization for dispersion of RNA-Seq data with cummeRbund.
+ - Added diffTable() method to CuffData and CuffFeatureSet objects to allow a 'one-table' snapshot of results for all Features (CuffData) or a set of Features (CuffFeatureSet). This table outputs key values including gene name,
+ gene short name, expression estimates and per-comparison fold-change, p-value, q-value, and significance values (yes/no). A convenient 'data-dump' function to merge across several tables.
+ - Added coercion methods for CuffGene objects to create GRanges and GRangeslist objects (more BioC friendly!). Will work on making this possible on CuffFeatureSet and CuffFeature objects as well.
+ - Added pass-through to select p.adjust method for getSig (method argument to getSig)
+ - Added ability to revert to cuffdiff q-values for specific paired-wise interrogations with getSig as opposed to re-calculating new ones (useCuffMTC; default=FALSE)
+
+ Notes:
+ - Removed generic for 'featureNames'. Now appropriately uses featureNames generic from Biobase. As a consequence, Biobase is now a dependency.
+ - Added passthrough to as.dist(...) in JSdist(...)
+ - Added 'logMode' argument to csClusterPlot.
+ - Added 'showPoints' argument to PCAplot to allow disabling of gene values in PCA plot. If false, only sample projections are plotted.
+ - Added 'facet' argument to expressionPlot to disable faceting by feature_id.
+ - shannon.entropy now uses log2 instead of log10 to constrain specificity scores between 0 and 1.
+
v1.99.6
Notes:
- 'annotation' and "annotation<-" generics were moved to BiocGenerics 0.3.2. Now using appropriate generic function, but requiring BiocGenerics >= 0.3.2
diff --git a/inst/doc/cummeRbund-example-workflow.pdf b/inst/doc/cummeRbund-example-workflow.pdf
index d074132..f0e5c9c 100644
Binary files a/inst/doc/cummeRbund-example-workflow.pdf and b/inst/doc/cummeRbund-example-workflow.pdf differ
diff --git a/inst/doc/cummeRbund-manual.R b/inst/doc/cummeRbund-manual.R
index ffb4ed8..8917425 100644
--- a/inst/doc/cummeRbund-manual.R
+++ b/inst/doc/cummeRbund-manual.R
@@ -299,7 +299,16 @@ head(repFpkm(TSS(myGenes)))
###################################################
-### code chunk number 38: geneset_plots_1
+### code chunk number 38: create_geneset_3
+###################################################
+myGeneset.pluri<-getGenes(cuff,myGeneIds,sampleIdList=c("hESC","iPS"))
+
+myGeneset.pluri
+
+
+
+###################################################
+### code chunk number 39: geneset_plots_1
###################################################
h<-csHeatmap(myGenes,cluster='both')
h
@@ -308,7 +317,7 @@ h.rep
###################################################
-### code chunk number 39: geneset_plots_heatmap
+### code chunk number 40: geneset_plots_heatmap
###################################################
h<-csHeatmap(myGenes,cluster='both')
h
@@ -318,20 +327,20 @@ print(h)
###################################################
-### code chunk number 40: geneset_plots_heatmap_rep
+### code chunk number 41: geneset_plots_heatmap_rep
###################################################
print(h.rep)
###################################################
-### code chunk number 41: geneset_plots_1.5
+### code chunk number 42: geneset_plots_1.5
###################################################
b<-expressionBarplot(myGenes)
b
###################################################
-### code chunk number 42: geneset_plots_barplot
+### code chunk number 43: geneset_plots_barplot
###################################################
b<-expressionBarplot(myGenes)
b
@@ -339,14 +348,14 @@ print(b)
###################################################
-### code chunk number 43: geneset_plots_2
+### code chunk number 44: geneset_plots_2
###################################################
s<-csScatter(myGenes,"Fibroblasts","hESC",smooth=T)
s
###################################################
-### code chunk number 44: geneset_plots_scatter
+### code chunk number 45: geneset_plots_scatter
###################################################
s<-csScatter(myGenes,"Fibroblasts","hESC",smooth=T)
s
@@ -354,14 +363,14 @@ print(s)
###################################################
-### code chunk number 45: geneset_plots_3
+### code chunk number 46: geneset_plots_3
###################################################
v<-csVolcano(myGenes,"Fibroblasts","hESC")
v
###################################################
-### code chunk number 46: geneset_plots_volcano
+### code chunk number 47: geneset_plots_volcano
###################################################
v<-csVolcano(myGenes,"Fibroblasts","hESC")
v
@@ -369,7 +378,7 @@ print(v)
###################################################
-### code chunk number 47: geneset_plots_4
+### code chunk number 48: geneset_plots_4
###################################################
ih<-csHeatmap(isoforms(myGenes),cluster='both',labRow=F)
ih
@@ -378,7 +387,7 @@ th
###################################################
-### code chunk number 48: geneset_plots_isoform_heatmap
+### code chunk number 49: geneset_plots_isoform_heatmap
###################################################
ih<-csHeatmap(isoforms(myGenes),cluster='both',labRow=F)
ih
@@ -388,26 +397,26 @@ print(ih)
###################################################
-### code chunk number 49: geneset_plots_TSS_heatmap
+### code chunk number 50: geneset_plots_TSS_heatmap
###################################################
print(th)
###################################################
-### code chunk number 50: geneset_plots_5
+### code chunk number 51: geneset_plots_5
###################################################
den<-csDendro(myGenes)
###################################################
-### code chunk number 51: geneset_plots_dendro
+### code chunk number 52: geneset_plots_dendro
###################################################
den<-csDendro(myGenes)
plot(den)
###################################################
-### code chunk number 52: gene_level_1
+### code chunk number 53: gene_level_1
###################################################
myGeneId<-"PINK1"
myGene<-getGene(cuff,myGeneId)
@@ -417,7 +426,7 @@ head(fpkm(isoforms(myGene)))
###################################################
-### code chunk number 53: gene_plots_1
+### code chunk number 54: gene_plots_1
###################################################
gl<-expressionPlot(myGene)
gl
@@ -433,7 +442,7 @@ gl.cds.rep
###################################################
-### code chunk number 54: gene_plots_line
+### code chunk number 55: gene_plots_line
###################################################
gl<-expressionPlot(myGene)
gl
@@ -450,28 +459,28 @@ gl.cds.rep
###################################################
-### code chunk number 55: gene_plots_replicate_line
+### code chunk number 56: gene_plots_replicate_line
###################################################
print(gl.rep)
###################################################
-### code chunk number 56: gene_plots_iso_replicate_line
+### code chunk number 57: gene_plots_iso_replicate_line
###################################################
print(gl.iso.rep)
###################################################
-### code chunk number 57: gene_plots_cds_replicate_line
+### code chunk number 58: gene_plots_cds_replicate_line
###################################################
print(gl.cds.rep)
###################################################
-### code chunk number 58: gene_plots_2
+### code chunk number 59: gene_plots_2
###################################################
gb<-expressionBarplot(myGene)
gb
@@ -480,7 +489,7 @@ gb.rep
###################################################
-### code chunk number 59: gene_plots_bar
+### code chunk number 60: gene_plots_bar
###################################################
gb<-expressionBarplot(myGene)
gb
@@ -490,20 +499,20 @@ print(gb)
###################################################
-### code chunk number 60: gene_plots_bar_rep
+### code chunk number 61: gene_plots_bar_rep
###################################################
print(gb.rep)
###################################################
-### code chunk number 61: gene_plots_3
+### code chunk number 62: gene_plots_3
###################################################
igb<-expressionBarplot(isoforms(myGene),replicates=T)
igb
###################################################
-### code chunk number 62: gene_plots_bar_isoforms
+### code chunk number 63: gene_plots_bar_isoforms
###################################################
igb<-expressionBarplot(isoforms(myGene),replicates=T)
igb
@@ -511,22 +520,37 @@ print(igb)
###################################################
-### code chunk number 63: features_1
+### code chunk number 64: gene_plots_4
+###################################################
+gp<-csPie(myGene,level="isoforms")
+gp
+
+
+###################################################
+### code chunk number 65: gene_plots_pie
+###################################################
+gp<-csPie(myGene,level="isoforms")
+gp
+print(gp)
+
+
+###################################################
+### code chunk number 66: features_1
###################################################
head(features(myGene))
###################################################
-### code chunk number 64: features_2
+### code chunk number 67: features_2
###################################################
genetrack<-makeGeneRegionTrack(myGene)
plotTracks(genetrack)
###################################################
-### code chunk number 65: features_3
+### code chunk number 68: features_3
###################################################
-trackList<-list()
+trackList<-list()
myStart<-min(features(myGene)$start)
myEnd<-max(features(myGene)$end)
myChr<-unique(features(myGene)$seqnames)
@@ -563,14 +587,14 @@ plotTracks(trackList,from=myStart-2000,to=myEnd+2000)
###################################################
-### code chunk number 66: sig_mat_1
+### code chunk number 69: sig_mat_1
###################################################
mySigMat<-sigMatrix(cuff,level='genes',alpha=0.05)
###################################################
-### code chunk number 67: sig_mat_plot_1
+### code chunk number 70: sig_mat_plot_1
###################################################
mySigMat<-sigMatrix(cuff,level='genes',alpha=0.05)
@@ -578,7 +602,7 @@ print(mySigMat)
###################################################
-### code chunk number 68: get_sig_1
+### code chunk number 71: get_sig_1
###################################################
mySigGeneIds<-getSig(cuff,alpha=0.05,level='genes')
head(mySigGeneIds)
@@ -586,7 +610,7 @@ length(mySigGeneIds)
###################################################
-### code chunk number 69: get_sig_2
+### code chunk number 72: get_sig_2
###################################################
hESC_vs_iPS.sigIsoformIds<-getSig(cuff,x='hESC',y='iPS',alpha=0.05,level='isoforms')
head(hESC_vs_iPS.sigIsoformIds)
@@ -594,7 +618,7 @@ length(hESC_vs_iPS.sigIsoformIds)
###################################################
-### code chunk number 70: get_sig_3
+### code chunk number 73: get_sig_3
###################################################
mySigGenes<-getGenes(cuff,mySigGeneIds)
mySigGenes
@@ -602,21 +626,21 @@ mySigGenes
###################################################
-### code chunk number 71: get_sig_4
+### code chunk number 74: get_sig_4
###################################################
mySigTable<-getSigTable(cuff,alpha=0.01,level='genes')
head(mySigTable,20)
###################################################
-### code chunk number 72: dist_heat_1
+### code chunk number 75: dist_heat_1
###################################################
myDistHeat<-csDistHeat(genes(cuff))
###################################################
-### code chunk number 73: dist_heat_plot_1
+### code chunk number 76: dist_heat_plot_1
###################################################
myDistHeat<-csDistHeat(genes(cuff))
@@ -624,14 +648,14 @@ print(myDistHeat)
###################################################
-### code chunk number 74: dist_heat_2
+### code chunk number 77: dist_heat_2
###################################################
myRepDistHeat<-csDistHeat(genes(cuff),replicates=T)
###################################################
-### code chunk number 75: dist_heat_plot_2
+### code chunk number 78: dist_heat_plot_2
###################################################
myRepDistHeat<-csDistHeat(genes(cuff),replicates=T)
@@ -639,7 +663,7 @@ print(myRepDistHeat)
###################################################
-### code chunk number 76: dim_reduction_1
+### code chunk number 79: dim_reduction_1
###################################################
genes.PCA<-PCAplot(genes(cuff),"PC1","PC2")
genes.MDS<-MDSplot(genes(cuff))
@@ -649,7 +673,7 @@ genes.MDS.rep<-MDSplot(genes(cuff),replicates=T)
###################################################
-### code chunk number 77: gene_PCA
+### code chunk number 80: gene_PCA
###################################################
genes.PCA<-PCAplot(genes(cuff),"PC1","PC2")
genes.MDS<-MDSplot(genes(cuff))
@@ -660,28 +684,28 @@ genes.MDS.rep<-MDSplot(genes(cuff),replicates=T)
###################################################
-### code chunk number 78: gene_MDS
+### code chunk number 81: gene_MDS
###################################################
print(genes.MDS)
###################################################
-### code chunk number 79: gene_PCA_rep
+### code chunk number 82: gene_PCA_rep
###################################################
print(genes.PCA.rep)
###################################################
-### code chunk number 80: gene_MDS_rep
+### code chunk number 83: gene_MDS_rep
###################################################
print(genes.MDS.rep)
###################################################
-### code chunk number 81: geneset_cluster_1
+### code chunk number 84: geneset_cluster_1
###################################################
ic<-csCluster(myGenes,k=4)
head(ic$cluster)
@@ -690,27 +714,27 @@ icp
###################################################
-### code chunk number 82: geneset_plots_cluster
+### code chunk number 85: geneset_plots_cluster
###################################################
print(icp)
###################################################
-### code chunk number 83: specificity_1
+### code chunk number 86: specificity_1
###################################################
myGenes.spec<-csSpecificity(myGenes)
head(myGenes.spec)
###################################################
-### code chunk number 84: similar_1
+### code chunk number 87: similar_1
###################################################
mySimilar<-findSimilar(cuff,"PINK1",n=20)
mySimilar.expression<-expressionPlot(mySimilar,logMode=T,showErrorbars=F)
###################################################
-### code chunk number 85: similar_plots_1
+### code chunk number 88: similar_plots_1
###################################################
mySimilar<-findSimilar(cuff,"PINK1",n=20)
mySimilar.expression<-expressionPlot(mySimilar,logMode=T,showErrorbars=F)
@@ -718,7 +742,7 @@ print(mySimilar.expression)
###################################################
-### code chunk number 86: similar_2
+### code chunk number 89: similar_2
###################################################
myProfile<-c(500,0,400)
mySimilar2<-findSimilar(cuff,myProfile,n=10)
@@ -726,7 +750,7 @@ mySimilar2.expression<-expressionPlot(mySimilar2,logMode=T,showErrorbars=F)
###################################################
-### code chunk number 87: similar_plots_2
+### code chunk number 90: similar_plots_2
###################################################
myProfile<-c(500,0,400)
mySimilar2<-findSimilar(cuff,myProfile,n=10)
@@ -735,13 +759,13 @@ print(mySimilar2.expression)
###################################################
-### code chunk number 88: close_connection
+### code chunk number 91: close_connection
###################################################
end<-sqliteCloseConnection(cuff at DB)
###################################################
-### code chunk number 89: session
+### code chunk number 92: session
###################################################
sessionInfo()
diff --git a/inst/doc/cummeRbund-manual.Rnw b/inst/doc/cummeRbund-manual.Rnw
index 25de050..0b146c8 100644
--- a/inst/doc/cummeRbund-manual.Rnw
+++ b/inst/doc/cummeRbund-manual.Rnw
@@ -6,7 +6,7 @@
%\VignettePackage{cummeRbund}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%\SweaveOpts{prefix.string='graphics/cummeRbund-manual'}
+%\SweaveOpts{prefix.string=graphics/cummeRbund-manual}
\documentclass[10pt]{article}
\usepackage[margin=0.75in]{geometry}
\usepackage{amsmath}
@@ -63,7 +63,7 @@
\title{CummeRbund: Visualization and Exploration of Cufflinks High-throughput Sequencing Data}
\author{Loyal A. Goff, Cole Trapnell, David Kelley}
-\date{September 24, 2012}
+\date{May 7, 2014}
\begin{document}
<<init, echo=FALSE>>=
options(width=65)
@@ -74,20 +74,22 @@ options(width=65)
\clearpage
\section{Requirements}
NOTE: cummeRbund 2.0 was designed in conjunction with the release of cufflinks 2.0. While we attempted to preserve backwards-compatability, it is highly recommended
-that you update your cufflinks version >=2.0 to take full advantage of the improvements in modeling, reporting, and visualization that have been incorporated.
+that you update your cufflinks installation to version $\ge$2.0 to take full
+advantage of the improvements in modeling, reporting, and visualization that have been incorporated.
\begin{itemize}
- \item Cufflinks $\ge$ v2.0.0
+ \item Cufflinks $\ge$ v2.0.1
\item SQLite
- \item R $\ge$ v2.7.0
+ \item R $\ge$ v3.0
\item Packages:
\begin{itemize}
\item \Rpackage{RSQLite}
- \item \Rpackage{ggplot2 v0.9.2}
+ \item \Rpackage{ggplot2 $\ge$ v0.9.3}
\item \Rpackage{reshape2}
\item \Rpackage{plyr}
\item \Rpackage{fastcluster}
\item \Rpackage{rtracklayer}
\item \Rpackage{Gviz}
+ \item \Rpackage{BiocGenerics $\ge$ 0.3.2}
\item Recommended:
\begin{itemize}
\item \Rpackage{Hmisc}
@@ -299,6 +301,9 @@ disp
\end{center}
\end{figure}
+Alternatively a call to \code{dispersionPlot(cuff)} directly will allow you to
+visualize the full model fit.
+
The squared coefficient of variation is a normalized measure of cross-replicate
variability that can be useful for evaluating the quality your RNA-seq data.
Differences in $CV^2$ can result in lower numbers of differentially expressed
@@ -629,6 +634,17 @@ head(repFpkm(TSS(myGenes)))
As of \Rpackage{cummeRbund} $v2.0$ \Rclass{CuffGeneSet} classes can be created from any type of identifier ('gene\_id','isoform\_id','TSS\_group\_id', or 'CDS\_id'). When you pass a list of identifiers that are not gene\_id to \Rmethod{getGenes()}, the function attempts to lookup the parent gene\_id for each feature and returns \emph{all} relevant
information for the given genes and all of their sub-features (not just the sub-features passed to \Rmethod{getGenes()}). If you are interested in just retrieving information for a given set of features, please use the new \Rmethod{getFeatures()} method described later.
+More recent versions of cummeRbund allow for subsetting of conditions as well,
+by passing a vector of condition names to getGenes using the \Rfunarg{sampleIdList}
+argument.
+
+<<create_geneset_3>>=
+myGeneset.pluri<-getGenes(cuff,myGeneIds,sampleIdList=c("hESC","iPS"))
+
+myGeneset.pluri
+
+@
+
\subsection{Geneset level plots}
There are several plotting functions available for gene-set-level visualization:
@@ -878,7 +894,26 @@ igb
<<gene_plots_3>>
print(igb)
@
- \includegraphics[width=0.4\textwidth]{cummeRbund-manual-gene_plots_bar_isoforms}
+
+ \includegraphics[width=0.4\textwidth]{cummeRbund-manual-gene_plots_bar_isoforms}}
+ \end{center}
+\end{figure}
+
+<<gene_plots_4,include=FALSE>>=
+gp<-csPie(myGene,level="isoforms")
+gp
+@
+
+\begin{figure}[htp]
+ \begin{center}
+ \subfloat[Pie charts showing relative proportion of individual isoforms for a
+ single gene across conditions.]{
+
+<<label=gene_plots_pie,fig=TRUE,echo=FALSE,include=FALSE>>=
+<<gene_plots_4>>
+print(gp)
+@
+ \includegraphics[width=0.4\textwidth]{cummeRbund-manual-gene_plots_pie}
}
\end{center}
@@ -887,6 +922,7 @@ print(igb)
\clearpage
\subsubsection{Gene Feature plots}
+
If you included both the genome build and gtfFile in your call to \Rmethod{readCufflinks()} then you will be able to access some of the transcript-structure level features that are now being integrated into cummeRbund. For now, these features are extended only to the single gene, \Rclass{CuffGene} objects.
Feature data are loaded into the \'features\' table of the cuffData.db database. When a \Rclass{CuffGene} object is created using \Rmethod{getGene()}, all relative features are selected from this table and a \'features\' slot is added to the resulting object.
@@ -905,8 +941,11 @@ plotTracks(genetrack)
We can then use some of the additional features from the \Rpackage{Gviz} package
to add additional tracks from an external data source.
+\textit{Note: This feature is now deprecated owing to developmental constraints.
+It is still provided with cummeRbund, but is no longer supported.}
+
<<features_3,fig=TRUE>>=
-trackList<-list()
+trackList<-list()
myStart<-min(features(myGene)$start)
myEnd<-max(features(myGene)$end)
myChr<-unique(features(myGene)$seqnames)
@@ -1114,6 +1153,10 @@ genes.MDS.rep<-MDSplot(genes(cuff),replicates=T)
\end{center}
\end{figure}
+CummeRbund also includes a convenience wrapper around the \Rpackage{NMFN}
+function \Rfunction{nnmf} for non-negative matrix factorization. You can use the
+\Rfunction{csNMF()} method for either CuffData and CuffFeatureSet objects.
+
\clearpage
\subsection{Partitioning}
diff --git a/inst/doc/cummeRbund-manual.pdf b/inst/doc/cummeRbund-manual.pdf
index 0d75cba..4c01705 100644
Binary files a/inst/doc/cummeRbund-manual.pdf and b/inst/doc/cummeRbund-manual.pdf differ
diff --git a/inst/reports/styles/report_style.css b/inst/reports/styles/report_style.css
new file mode 100644
index 0000000..9b66fa2
--- /dev/null
+++ b/inst/reports/styles/report_style.css
@@ -0,0 +1,124 @@
+body, td {
+ font-family: sans-serif;
+ background-color: white;
+ font-size: 12px;
+ margin: 8px;
+}
+
+tt, code, pre {
+ font-family: 'DejaVu Sans Mono', 'Droid Sans Mono', 'Lucida Console', Consolas, Monaco, monospace;
+}
+
+h1 {
+ font-size:2.2em;
+}
+
+h2 {
+ font-size:1.8em;
+}
+
+h3 {
+ font-size:1.4em;
+}
+
+h4 {
+ font-size:1.0em;
+}
+
+h5 {
+ font-size:0.9em;
+}
+
+h6 {
+ font-size:0.8em;
+}
+
+a:visited {
+ color: rgb(50%, 0%, 50%);
+}
+
+pre {
+ margin-top: 0;
+ max-width: 95%;
+ border: 1px solid #ccc;
+ white-space: pre-wrap;
+}
+
+pre code {
+ display: block; padding: 0.5em;
+}
+
+code.r, code.cpp {
+ background-color: #F8F8F8;
+}
+
+table, td, th {
+ border: none;
+}
+
+blockquote {
+ color:#666666;
+ margin:0;
+ padding-left: 1em;
+ border-left: 0.5em #EEE solid;
+}
+
+hr {
+ height: 0px;
+ border-bottom: none;
+ border-top-width: thin;
+ border-top-style: dotted;
+ border-top-color: #999999;
+}
+
+ at media print {
+ * {
+ background: transparent !important;
+ color: black !important;
+ filter:none !important;
+ -ms-filter: none !important;
+ }
+
+ body {
+ font-size:12pt;
+ max-width:100%;
+ }
+
+ a, a:visited {
+ text-decoration: underline;
+ }
+
+ hr {
+ visibility: hidden;
+ page-break-before: always;
+ }
+
+ pre, blockquote {
+ padding-right: 1em;
+ page-break-inside: avoid;
+ }
+
+ tr, img {
+ page-break-inside: avoid;
+ }
+
+ img {
+ max-width: 100% !important;
+ }
+
+ @page :left {
+ margin: 15mm 20mm 15mm 10mm;
+ }
+
+ @page :right {
+ margin: 15mm 10mm 15mm 20mm;
+ }
+
+ p, h2, h3 {
+ orphans: 3; widows: 3;
+ }
+
+ h2, h3 {
+ page-break-after: avoid;
+ }
+}
\ No newline at end of file
diff --git a/man/CuffData-class.Rd b/man/CuffData-class.Rd
index 89364a0..56dff85 100644
--- a/man/CuffData-class.Rd
+++ b/man/CuffData-class.Rd
@@ -5,9 +5,7 @@
\alias{dim,CuffData-method}
\alias{getFeatures,CuffData-method}
\alias{DB,CuffData-method}
-\alias{diffTable}
\alias{makeRnk}
-\alias{diffTable,CuffData-method}
\alias{makeRnk,CuffData-method}
\alias{annotation,CuffData-method}
diff --git a/man/CuffFeature-class.Rd b/man/CuffFeature-class.Rd
index f7d2214..b97311b 100644
--- a/man/CuffFeature-class.Rd
+++ b/man/CuffFeature-class.Rd
@@ -29,6 +29,7 @@ Objects can be created by calls of the form \code{new("CuffFeature", annotation,
\section{Methods}{
\describe{
\item{fpkmMatrix}{\code{signature(object="CuffFeature")}: ...}
+ \item{repFpkmMatrix}{\code{signature(object = "CuffFeature")}: ... }
\item{length}{\code{signature(x = "CuffFeature")}: ... }
}
}
diff --git a/man/CuffGene-class.Rd b/man/CuffGene-class.Rd
index 54011a6..35fdcbb 100644
--- a/man/CuffGene-class.Rd
+++ b/man/CuffGene-class.Rd
@@ -14,6 +14,8 @@
\alias{makeGeneRegionTrack,CuffGene-method}
\alias{genePlot}
\alias{genePlot,CuffGene-method}
+\alias{csPie}
+\alias{csPie,CuffGene-method}
\title{Class "CuffGene"}
\description{
@@ -47,6 +49,7 @@ Class \code{"\linkS4class{CuffFeature}"}, directly.
\item{length}{\code{signature(object="CuffFeature")}: Part of length validation (internal use only)}
\item{makeGeneRegionTrack}{\code{signature(object="CuffFeature")}: Creates a GeneRegionTrack object (see package Gviz) from a CuffGene object.}
\item{genePlot}{\code{signature(object="CuffFeature")}: Internal use only.}
+ \item{csPie}{\code{signature(object="CuffGene")}: Allows for visualization of relative isoform proportion as a pie chart by condition (or optionally as stacked bar charts by adding + coord_cartesian()}
}
}
diff --git a/man/CuffSet-class.Rd b/man/CuffSet-class.Rd
index 59f00a5..5eb24d5 100644
--- a/man/CuffSet-class.Rd
+++ b/man/CuffSet-class.Rd
@@ -11,6 +11,8 @@
\alias{promoters}
\alias{splicing}
\alias{relCDS}
+\alias{conditions}
+\alias{varModel}
\alias{DB,CuffSet-method}
\alias{genes,CuffSet-method}
@@ -20,6 +22,8 @@
\alias{promoters,CuffSet-method}
\alias{splicing,CuffSet-method}
\alias{relCDS,CuffSet-method}
+\alias{conditions,CuffSet-method}
+\alias{varModel,CuffSet-method}
\alias{[,CuffSet-method}
@@ -41,9 +45,9 @@ Available methods are primary accessors to retrieve CuffGeneSet or CuffGene obje
\item{\code{TSS}:}{Object of class \code{"CuffData"} ~~ }
\item{\code{CDS}:}{Object of class \code{"CuffData"} ~~ }
\item{\code{promoters}:}{Object of class \code{"CuffDist"} ~~ }
- \item{\code{runInfo}:}{Object of class \code{"data.frame"} ~~ }
\item{\code{splicing}:}{Object of class \code{"CuffDist"} ~~ }
\item{\code{relCDS}:}{Object of class \code{"CuffDist"} ~~ }
+ \item{\code{conditions}:}{Object of class \code{"data.frame"} ~~ }
}
}
\section{Methods}{
@@ -62,6 +66,7 @@ Available methods are primary accessors to retrieve CuffGeneSet or CuffGene obje
\item{promoters}{\code{signature(object = "CuffSet")}: Access @promoters slot }
\item{splicing}{\code{signature(object = "CuffSet")}: Access @splicing slot }
\item{relCDS}{\code{signature(object = "CuffSet")}: Access @relCDS slot }
+ \item{varModel}{\code{signature(object = "CuffSet")}: Access varModel info }
}
}
diff --git a/man/JSdist.Rd b/man/JSdist.Rd
index 150b9ca..f83a5c5 100644
--- a/man/JSdist.Rd
+++ b/man/JSdist.Rd
@@ -9,13 +9,16 @@ JSdist takes a matrix of expression probabilites (calculated directly or output
of the pairwise Jensen-Shannon distances between columns
}
\usage{
-JSdist(mat)
+JSdist(mat,...)
}
\arguments{
\item{mat}{
A matrix of expression probabilities (e.g. from makeprobs())
}
+ \item{\dots}{
+Passthrough argument to as.dist()
+}
}
\details{
Returns pairwise Jensen-Shannon distance (in the form of a dist object) for a matrix of probabilities (by column)
diff --git a/man/QCplots.Rd b/man/QCplots.Rd
index 749d374..03faa48 100644
--- a/man/QCplots.Rd
+++ b/man/QCplots.Rd
@@ -12,7 +12,7 @@ A collection of ggplot2 visualizations for quality control assessment of cuffdif
}
\usage{
-\S4method{fpkmSCVPlot}{CuffData}(object,FPKMLowerBound=1)
+\S4method{fpkmSCVPlot}{CuffData}(object,FPKMLowerBound=1, showPool = FALSE)
}
\arguments{
@@ -22,6 +22,9 @@ An object of class CuffData.
\item{FPKMLowerBound}{
A lower limit cutoff for FPKM values from which a fit of squared Coefficient of variation (default: 1)
}
+ \item{showPool}{
+Logical argument whether to display variability across all replicates independent of condition (TRUE) or the cross-replicate variability for each condition (FALSE)
+}
}
\details{
None
diff --git a/man/csCluster.Rd b/man/csCluster.Rd
index b5a1383..aa7a42b 100644
--- a/man/csCluster.Rd
+++ b/man/csCluster.Rd
@@ -10,7 +10,7 @@ Returns a ggplot2 plot object with geom_line layer plotting FPKM values over con
This does not return any of the clustering information directly, but if you want it, you can retrieve it from the ggplot object returned.
}
\usage{
-\S4method{csCluster}{CuffFeatureSet}(object,k,logMode=T,pseudocount=1,...)
+\S4method{csCluster}{CuffFeatureSet}(object,k,logMode=T,method = "none",pseudocount=1,...)
}
\arguments{
\item{object}{
@@ -22,6 +22,9 @@ Number of pre-defined clusters to attempt to find.
\item{logMode}{
A logical value whether or not to log-transform the FPKM values prior to clustering.
}
+ \item{method}{
+Distance function to use when computing cluster solution. Default "none" will use the Jensen-Shannon distance (JSdist). Provide a function that returns a dist object on rows.
+}
\item{pseudocount}{
Value added to FPKM to avoid log-transform issues.
}
diff --git a/man/csClusterPlot.Rd b/man/csClusterPlot.Rd
index 983b622..e74d7fd 100644
--- a/man/csClusterPlot.Rd
+++ b/man/csClusterPlot.Rd
@@ -8,7 +8,7 @@ csClusterPlot
Replaces the default plotting behavior of the old csCluster. Takes as an argument the output of csCluster and plots expression profiles of features facet by cluster.
}
\usage{
-csClusterPlot(clustering, pseudocount=1.0,drawSummary=TRUE,sumFun=mean_cl_boot)
+csClusterPlot(clustering, pseudocount=1.0,logMode=FALSE,drawSummary=TRUE,sumFun=mean_cl_boot)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
@@ -18,6 +18,9 @@ The output of csCluster. (Must be the output of csCluster. Only this data format
\item{pseudocount}{
Value added to FPKM to avoid log transformation issues.
}
+ \item{logMode}{
+Logical argument whether to plot FPKM with log axis (Y-axis).
+}
\item{drawSummary}{
Logical value whether or not to draw a summary line for each cluster (by default this is the cluster mean)
}
diff --git a/man/csHeatmap.Rd b/man/csHeatmap.Rd
index 1cd38b6..23afd0d 100644
--- a/man/csHeatmap.Rd
+++ b/man/csHeatmap.Rd
@@ -1,6 +1,8 @@
\name{csHeatmap}
\alias{csHeatmap}
+\alias{csFoldChangeHeatmap}
\alias{csHeatmap,CuffFeatureSet-method}
+\alias{csFoldChangeHeatmap,CuffFeatureSet-method}
\title{
csHeatmap
@@ -11,12 +13,20 @@ Creates a ggplot plot object with a geom_tile layer of FPKM values per feature a
\usage{
\S4method{csHeatmap}{CuffFeatureSet}(object, rescaling='none', clustering='none', labCol=T, labRow=T, logMode=T, pseudocount=1.0,
border=FALSE, heatscale= c(low='lightyellow',mid='orange',high='darkred'), heatMidpoint=NULL, fullnames = T, replicates=FALSE,method='none', ...)
+\S4method{csFoldChangeHeatmap}{CuffFeatureSet}(object, control_condition, replicate_num=NULL, clustering='none', labCol=T, labRow=T, logMode=F, pseudocount=1.0,
+ border=FALSE, heatscale=c(low='steelblue',mid='white',high='tomato'), heatMidpoint=0,fullnames=T,replicates=FALSE,method='none',heatRange=3, ...)
}
\arguments{
\item{object}{
An object of class 'CuffFeatureSet' or 'CuffGeneSet'
}
+ \item{control_condition}{
+A character argument indicating which condition should be used as the denominator for fold change. (e.g. "Day0", "Control", etc)
+}
+ \item{replicate_num}{
+If replicates == TRUE, you must specify both a control condition and a replicate number to use as the denominator.
+}
\item{rescaling}{
Rescaling can either be 'row' or 'column' OR you can pass rescale a function that operates on a matrix to do your own rescaling. Default is 'none'.
}
@@ -53,6 +63,9 @@ A logical value whether or not to plot individual replicates or aggregate condit
\item{method}{
Function to be used for clustering. Default is JS-distance. You can pass your own function to this argument as long as the output is an instance of the 'dist' class and is applied to the rows of the input matrix.
}
+ \item{heatRange}{
+Numerical argument for upper bound on log fold change to be visualized.
+}
\item{\dots}{
Additional arguments to csHeatmap
}
diff --git a/man/csScatter.Rd b/man/csScatter.Rd
index 2a0dc59..b0715d7 100644
--- a/man/csScatter.Rd
+++ b/man/csScatter.Rd
@@ -27,7 +27,7 @@ Sample name for x axis
Sample name for y axis
}
\item{logMode}{
-Logical argument to log2-transform data (default: T )
+Logical argument to render axes on log10 scale (default: T )
}
\item{replicates}{
Logical argument whether or not to draw individual replicate values instead of condition values. (default: T )
diff --git a/man/diffData.Rd b/man/diffData.Rd
index a0fd7f7..4616715 100644
--- a/man/diffData.Rd
+++ b/man/diffData.Rd
@@ -1,6 +1,9 @@
\name{diffData}
\alias{diffData}
\alias{diffData,CuffData-method}
+\alias{diffTable}
+\alias{diffTable,CuffData-method}
+\alias{diffTable,CuffFeatureSet-method}
\title{
Differential comparison data
@@ -10,6 +13,7 @@ An accessor method to retrieve differential expression data from a 'CuffData', '
}
\usage{
\S4method{diffData}{CuffData}(object, x, y, features=FALSE)
+\S4method{diffTable}{CuffData}(object,logCutoffValue=99999)
}
\arguments{
@@ -26,6 +30,9 @@ See 'x'
\item{features}{
A logical value that returns all feature-level data as part of data.frame when true. object must be of class 'CuffData'.
}
+ \item{logCutoffValue}{
+Cutoff value for FC estimates to convert to [-]Inf values. Should never really be needed...
+}
\item{\dots}{
Additional arguments.
diff --git a/man/dimensionality.Rd b/man/dimensionality.Rd
index 331a2dc..a399b3c 100644
--- a/man/dimensionality.Rd
+++ b/man/dimensionality.Rd
@@ -1,8 +1,10 @@
\name{Dimensionality Reduction}
\alias{MDSplot}
\alias{MDSplot,CuffData-method}
+\alias{MDSplot,CuffFeatureSet-method}
\alias{PCAplot}
\alias{PCAplot,CuffData-method}
+\alias{PCAplot,CuffFeatureSet-method}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
Dimensionality reduction utilities
@@ -12,7 +14,7 @@ Dimensionality reduction plots for feature selection and extraction for cummeRbu
}
\usage{
\S4method{MDSplot}{CuffData}(object,replicates=FALSE,logMode=TRUE,pseudocount=1.0)
-\S4method{PCAplot}{CuffData}(object,x="PC1", y="PC2",replicates=FALSE,pseudocount=1.0,scale=TRUE,...)
+\S4method{PCAplot}{CuffData}(object,x="PC1", y="PC2",replicates=FALSE,pseudocount=1.0,scale=TRUE,showPoints = TRUE,...)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
@@ -37,6 +39,9 @@ A logical value to indicate whether or not individual replicate expression estim
\item{scale}{
For PCAplot, a logical value passed directly to prcomp.
}
+ \item{showPoints}{
+For PCAplot, a logical value whether or not to display individual gene values on final PCA plot.
+}
\item{\dots}{
Additional passthrough arguments (may not be fully implemented yet).
}
diff --git a/man/dispersionPlot.Rd b/man/dispersionPlot.Rd
index 014aee2..bb4f237 100644
--- a/man/dispersionPlot.Rd
+++ b/man/dispersionPlot.Rd
@@ -1,6 +1,7 @@
\name{dispersionPlot}
\alias{dispersionPlot}
\alias{dispersionPlot,CuffData-method}
+\alias{dispersionPlot,CuffSet-method}
\title{
Mean count vs dispersion plot
@@ -10,6 +11,7 @@ A scatter plot comparing the mean counts against the estimated dispersion for a
}
\usage{
\S4method{dispersionPlot}{CuffData}(object)
+\S4method{dispersionPlot}{CuffSet}(object)
}
\arguments{
diff --git a/man/exploratory.Rd b/man/exploratory.Rd
new file mode 100644
index 0000000..7df95ce
--- /dev/null
+++ b/man/exploratory.Rd
@@ -0,0 +1,63 @@
+\name{Exploratory Analysis}
+\alias{csNMF}
+\alias{csNMF,CuffData-method}
+\alias{csNMF,CuffFeatureSet-method}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{
+Methods for Exploratory Analysis in cummeRbund
+}
+\description{
+Exploratory analysis methods for cummeRbund RNA-Seq data.
+}
+\usage{
+\S4method{csNMF}{CuffData}(object,k,logMode=T,pseudocount=1,maxiter=1000,replicates=FALSE,fullnames=FALSE)
+\S4method{csNMF}{CuffFeatureSet}(object,k,logMode=T,pseudocount=1,maxiter=1000,replicates=FALSE,fullnames=FALSE)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{object}{
+The output of class CuffData or CuffFeatureSet from which to draw expression estimates. (e.g. genes(cuff) or custom feature set via getGenes() or getFeatures() )
+}
+ \item{k}{
+rank value for factorization
+}
+ \item{logMode}{
+Logical value whether or not to use log-transformed FPKM values. [Default: TRUE]
+}
+ \item{pseudocount}{
+Value added to FPKM to avoid log transformation issues.
+}
+ \item{maxiter}{
+Maximum number of iterations for factorization [Default: 1000]
+}
+ \item{replicates}{
+A logical value to indicate whether or not individual replicate expression estimates will be used.
+}
+ \item{fullnames}{
+Logical passthrough value to fpkmMatrix whether or not to concatenate gene_short_name with tracking_id. [Default: FALSE]
+}
+}
+\details{
+csNMF is a convenience method to invoke the nnmf() method from package:NMFN. This performs non-negative matrix factorization on the provided data and can be useful for many downstream applications.
+
+}
+\value{
+csNMF returns W, H - decomposed matrices of input FPKM values. (See package:NMFN for details)
+}
+\references{
+None.
+}
+\author{
+Loyal A. Goff
+}
+\note{
+None.
+}
+
+
+\examples{
+ data(sampleData)
+ csNMF(sampleGeneSet,4)
+
+}
+
diff --git a/man/expressionPlot.Rd b/man/expressionPlot.Rd
index c1ce315..823ab3f 100644
--- a/man/expressionPlot.Rd
+++ b/man/expressionPlot.Rd
@@ -11,7 +11,7 @@ Expression Plot
A line plot (optionally with confidence intervals) detailing FPKM expression levels across conditions for a given gene(s) or feature(s)
}
\usage{
-\S4method{expressionPlot}{CuffFeature}(object, logMode=FALSE, pseudocount=1.0, drawSummary=FALSE, sumFun=mean_cl_boot, showErrorbars=TRUE, showStatus=TRUE, replicates=FALSE, ...)
+\S4method{expressionPlot}{CuffFeature}(object, logMode=FALSE, pseudocount=1.0, drawSummary=FALSE, sumFun=mean_cl_boot, showErrorbars=TRUE, showStatus=TRUE, replicates=FALSE, facet = TRUE,...)
}
\arguments{
@@ -39,6 +39,9 @@ A logical value whether or not to draw visual queues for quantification status o
\item{replicates}{
A logical value whether or not to plot individual replicates or aggregate condition values.
}
+ \item{facet}{
+A logical value whether or not to facet the plot by feature id (default=TRUE).
+}
\item{\dots}{
Additional arguments
}
diff --git a/man/featureNames.Rd b/man/featureNames.Rd
index 96e52cf..d57ddae 100644
--- a/man/featureNames.Rd
+++ b/man/featureNames.Rd
@@ -1,5 +1,4 @@
\name{featureNames}
-\alias{featureNames}
\alias{featureNames,CuffData-method}
\title{
diff --git a/man/fpkmMatrix.Rd b/man/fpkmMatrix.Rd
index 64a742e..8972ccb 100644
--- a/man/fpkmMatrix.Rd
+++ b/man/fpkmMatrix.Rd
@@ -5,6 +5,7 @@
\alias{repFpkmMatrix}
\alias{repFpkmMatrix,CuffData-method}
\alias{repFpkmMatrix,CuffFeatureSet-method}
+\alias{repFpkmMatrix,CuffFeature-method}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
diff --git a/man/getSig.Rd b/man/getSig.Rd
index dcbb63a..fc9c652 100644
--- a/man/getSig.Rd
+++ b/man/getSig.Rd
@@ -10,7 +10,7 @@ getSig
Returns the identifiers of significant genes in a vector format.
}
\usage{
-\S4method{getSig}{CuffSet}(object,x,y,alpha=0.05,level='genes')
+\S4method{getSig}{CuffSet}(object,x,y,alpha=0.05,level='genes',method="BH", useCuffMTC=FALSE)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
@@ -29,6 +29,12 @@ An alpha value by which to filter multiple-testing corrected q-values to determi
\item{level}{
Feature level to be queried for significance (must be one of c('genes','isoforms','TSS','CDS')
}
+ \item{method}{
+Multiple testing method to be used for correction. (default: "BH")
+}
+ \item{useCuffMTC}{
+Logical vector whether or not to use the multiple-testing corrected q-values from the cuffdiff analysis directly, or calculate new q-values from a subset of tests.
+}
}
\details{
diff --git a/man/readCufflinks.Rd b/man/readCufflinks.Rd
index 46700b0..c6133bb 100644
--- a/man/readCufflinks.Rd
+++ b/man/readCufflinks.Rd
@@ -16,7 +16,8 @@ readCufflinks(dir = getwd(), dbFile = "cuffData.db", gtfFile = NULL,
CDSFPKM = "cds.fpkm_tracking", CDSExpDiff = "cds_exp.diff", CDSCount="cds.count_tracking", CDSRep="cds.read_group_tracking",
CDSDiff = "cds.diff",
promoterFile = "promoters.diff",
- splicingFile = "splicing.diff",
+ splicingFile = "splicing.diff",
+ varModelFile = "var_model.info",
driver = "SQLite",
genome = NULL,
rebuild = FALSE,verbose=FALSE, ...)
@@ -95,6 +96,9 @@ promoters.diff file (distribution tests on promoters)
\item{splicingFile}{
splicing.diff (distribution tests on isoforms)
}
+ \item{varModelFile}{
+varModel.info (emitted in cuffdiff >= v2.1)
+}
\item{driver}{
Driver for backend database. (Currently only "SQLite" is supported).
}
diff --git a/man/replicates.Rd b/man/replicates.Rd
index 65c6b33..b2d9495 100644
--- a/man/replicates.Rd
+++ b/man/replicates.Rd
@@ -3,6 +3,7 @@
\alias{replicates,CuffSet-method}
\alias{replicates,CuffData-method}
\alias{replicates,CuffFeatureSet-method}
+\alias{replicates,CuffFeature-method}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
Get replicate sample list from CuffData object
diff --git a/vignettes/cummeRbund-manual.Rnw b/vignettes/cummeRbund-manual.Rnw
index 25de050..0b146c8 100644
--- a/vignettes/cummeRbund-manual.Rnw
+++ b/vignettes/cummeRbund-manual.Rnw
@@ -6,7 +6,7 @@
%\VignettePackage{cummeRbund}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%\SweaveOpts{prefix.string='graphics/cummeRbund-manual'}
+%\SweaveOpts{prefix.string=graphics/cummeRbund-manual}
\documentclass[10pt]{article}
\usepackage[margin=0.75in]{geometry}
\usepackage{amsmath}
@@ -63,7 +63,7 @@
\title{CummeRbund: Visualization and Exploration of Cufflinks High-throughput Sequencing Data}
\author{Loyal A. Goff, Cole Trapnell, David Kelley}
-\date{September 24, 2012}
+\date{May 7, 2014}
\begin{document}
<<init, echo=FALSE>>=
options(width=65)
@@ -74,20 +74,22 @@ options(width=65)
\clearpage
\section{Requirements}
NOTE: cummeRbund 2.0 was designed in conjunction with the release of cufflinks 2.0. While we attempted to preserve backwards-compatability, it is highly recommended
-that you update your cufflinks version >=2.0 to take full advantage of the improvements in modeling, reporting, and visualization that have been incorporated.
+that you update your cufflinks installation to version $\ge$2.0 to take full
+advantage of the improvements in modeling, reporting, and visualization that have been incorporated.
\begin{itemize}
- \item Cufflinks $\ge$ v2.0.0
+ \item Cufflinks $\ge$ v2.0.1
\item SQLite
- \item R $\ge$ v2.7.0
+ \item R $\ge$ v3.0
\item Packages:
\begin{itemize}
\item \Rpackage{RSQLite}
- \item \Rpackage{ggplot2 v0.9.2}
+ \item \Rpackage{ggplot2 $\ge$ v0.9.3}
\item \Rpackage{reshape2}
\item \Rpackage{plyr}
\item \Rpackage{fastcluster}
\item \Rpackage{rtracklayer}
\item \Rpackage{Gviz}
+ \item \Rpackage{BiocGenerics $\ge$ 0.3.2}
\item Recommended:
\begin{itemize}
\item \Rpackage{Hmisc}
@@ -299,6 +301,9 @@ disp
\end{center}
\end{figure}
+Alternatively a call to \code{dispersionPlot(cuff)} directly will allow you to
+visualize the full model fit.
+
The squared coefficient of variation is a normalized measure of cross-replicate
variability that can be useful for evaluating the quality your RNA-seq data.
Differences in $CV^2$ can result in lower numbers of differentially expressed
@@ -629,6 +634,17 @@ head(repFpkm(TSS(myGenes)))
As of \Rpackage{cummeRbund} $v2.0$ \Rclass{CuffGeneSet} classes can be created from any type of identifier ('gene\_id','isoform\_id','TSS\_group\_id', or 'CDS\_id'). When you pass a list of identifiers that are not gene\_id to \Rmethod{getGenes()}, the function attempts to lookup the parent gene\_id for each feature and returns \emph{all} relevant
information for the given genes and all of their sub-features (not just the sub-features passed to \Rmethod{getGenes()}). If you are interested in just retrieving information for a given set of features, please use the new \Rmethod{getFeatures()} method described later.
+More recent versions of cummeRbund allow for subsetting of conditions as well,
+by passing a vector of condition names to getGenes using the \Rfunarg{sampleIdList}
+argument.
+
+<<create_geneset_3>>=
+myGeneset.pluri<-getGenes(cuff,myGeneIds,sampleIdList=c("hESC","iPS"))
+
+myGeneset.pluri
+
+@
+
\subsection{Geneset level plots}
There are several plotting functions available for gene-set-level visualization:
@@ -878,7 +894,26 @@ igb
<<gene_plots_3>>
print(igb)
@
- \includegraphics[width=0.4\textwidth]{cummeRbund-manual-gene_plots_bar_isoforms}
+
+ \includegraphics[width=0.4\textwidth]{cummeRbund-manual-gene_plots_bar_isoforms}}
+ \end{center}
+\end{figure}
+
+<<gene_plots_4,include=FALSE>>=
+gp<-csPie(myGene,level="isoforms")
+gp
+@
+
+\begin{figure}[htp]
+ \begin{center}
+ \subfloat[Pie charts showing relative proportion of individual isoforms for a
+ single gene across conditions.]{
+
+<<label=gene_plots_pie,fig=TRUE,echo=FALSE,include=FALSE>>=
+<<gene_plots_4>>
+print(gp)
+@
+ \includegraphics[width=0.4\textwidth]{cummeRbund-manual-gene_plots_pie}
}
\end{center}
@@ -887,6 +922,7 @@ print(igb)
\clearpage
\subsubsection{Gene Feature plots}
+
If you included both the genome build and gtfFile in your call to \Rmethod{readCufflinks()} then you will be able to access some of the transcript-structure level features that are now being integrated into cummeRbund. For now, these features are extended only to the single gene, \Rclass{CuffGene} objects.
Feature data are loaded into the \'features\' table of the cuffData.db database. When a \Rclass{CuffGene} object is created using \Rmethod{getGene()}, all relative features are selected from this table and a \'features\' slot is added to the resulting object.
@@ -905,8 +941,11 @@ plotTracks(genetrack)
We can then use some of the additional features from the \Rpackage{Gviz} package
to add additional tracks from an external data source.
+\textit{Note: This feature is now deprecated owing to developmental constraints.
+It is still provided with cummeRbund, but is no longer supported.}
+
<<features_3,fig=TRUE>>=
-trackList<-list()
+trackList<-list()
myStart<-min(features(myGene)$start)
myEnd<-max(features(myGene)$end)
myChr<-unique(features(myGene)$seqnames)
@@ -1114,6 +1153,10 @@ genes.MDS.rep<-MDSplot(genes(cuff),replicates=T)
\end{center}
\end{figure}
+CummeRbund also includes a convenience wrapper around the \Rpackage{NMFN}
+function \Rfunction{nnmf} for non-negative matrix factorization. You can use the
+\Rfunction{csNMF()} method for either CuffData and CuffFeatureSet objects.
+
\clearpage
\subsection{Partitioning}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-bioc-cummerbund.git
More information about the debian-med-commit
mailing list