sydSeq_1.1.7.tar

sydSeq/ 40777      0      0           0 12526447621   5356 5sydSeq/DESCRIPTION100777      0      0        1010 12552232711   7143 0Package: sydSeq
Type: Package
Title: Tools for analysing RNA-Seq Data
Version: 1.1.7
Date: 2015-07-14
Author: Ellis Patrick
Maintainer: Ellis Patrick <ellis.patrick@sydney.edu.au>
Description: sydSeq is a set of tools that can help with the analysis of RNA-Seq data. 
URL: http://www.ellispatrick.com/r-packages/
Depends: R (>= 2.14.0), Biobase, GenomicFeatures, GenomicRanges
Imports: limma, DESeq, edgeR, googleVis
License: GPL (>= 2)
LazyData: no
LazyLoad: yes
Packaged: 2015-07-14 03:39:06 UTC; ellis
sydSeq/man/ 40777      0      0           0 12526447715   6135 5sydSeq/man/exClust.Rd100777      0      0        3007 12405237430  10136 0\name{exClust}
\Rdversion{1}
\alias{exClust}
\title{Estimates constitutive exons using both annotation and data}
\description{
Use clustering to estimate constitutive exons.
}
\usage{
exClust(simpleGRanges, cut = 1.5)
}
\arguments{
\item{simpleGRanges}{An object of type \code{GRanges} which has been simplified and had a counts matrix merged.}
  \item{cut}{The level to cut clustering tree.}
}
  \details{
Method clusters exons and cuts clustering tree at value \code{cut}.
}
\value{
An object of type \code{GRanges} containing the data specific constitutive exons.
}
\author{ Ellis Patrick }
\examples{
\dontrun{
hsapiensEnsembl <-makeTranscriptDbFromBiomart(biomart = "ensembl",dataset = "hsapiens_gene_ensembl")
saveFeatures(hsapiensEnsembl, file="hsapiensEnsembl.sqlite")
simpleEnsembl = simplifyGFeatures(hsapiensEnsembl)

###  Read in Data from bam files
bamDir = '.'
bamFls <- list.files(bamDir, "bam$", full=TRUE)
names(bamFls) <- sub("\\..*", "", basename(bamFls))
counter <- function(fl, gnModel)
{
    aln <- readGappedAlignments(fl)
	seqlevels(aln) <- sub("chr","", seqlevels(aln))
    strand(aln) <- "*" # for strand-blind sample prep protocol
    hits <- countOverlaps(aln, gnModel)
    counts <- countOverlaps(gnModel, aln[hits==1])
    names(counts) <- names(gnModel)
    counts
}
counts <- sapply(bamFls, counter, simpleEnsembl)

simpleEnsembl= mergeCounts(counts,simpleEnsembl)
save(simpleEnsembl,file = 'simpleEnsembl.RData')
exClustAnno = exClust(simpleEnsembl,cut = 2)
}
}
sydSeq/man/geneCounts.Rd100777      0      0        1007 12405237430  10617 0\name{geneCounts}
\Rdversion{1}
\alias{geneCounts}
\title{Gene counts.}
\description{
Get gene counts from GRanges object
}
\usage{
geneCounts(Anno)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
\item{Anno}{An annotation object of type \code{GRanges} which has been simplified and had a counts matrix merged.}
}
\value{
An object of type \code{matrix}
}
\author{ Ellis Patrick }

% Add one or more standard keywords, see file 'KEYWORDS' in the
% R documentation directory.
sydSeq/man/makeUI.Rd100777      0      0        1655 12405237430   7671 0\name{makeUI}
\Rdversion{1}
\alias{makeUI}
\title{Make UI annotation}
\description{
Make a UI annotation from GRanges object.
}
\usage{
makeUI(simpleGRanges, gFeature)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
\item{simpleGRanges}{An object of type \code{GRanges} which has been simplified}
  \item{gFeature}{a genome annotation of type \code{gFeature}}
}
  \details{

Method makes UI definition.
}
\value{
An object of type \code{GRanges} containing UI constitutive exons.
}
\author{ Ellis Patrick }
\examples{
\dontrun{
hsapiensEnsembl <-makeTranscriptDbFromBiomart(biomart = "ensembl",dataset = "hsapiens_gene_ensembl")
saveFeatures(hsapiensEnsembl, file="hsapiensEnsembl.sqlite")
simpleEnsembl <- simplifyGFeatures(hsapiensEnsembl)
annoUI <- makeUI(simpleEnsembl,hsapiensEnsembl)
}}
% Add one or more standard keywords, see file 'KEYWORDS' in the
% R documentation directory.
sydSeq/man/mergeCounts.Rd100777      0      0        1067 12405237430  11006 0\name{mergeCounts}
\Rdversion{1}
\alias{mergeCounts}
\title{Merge Counts}
\description{
Merges a counts matrix with a GRanges object.
}
\usage{
mergeCounts(c, GRanges)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
\item{c}{a matrix of counts}
  \item{cut}{An object of type \code{GRanges} }
}
  \details{

Merges a GRanges object with a counts matrix.
}
\value{
An object of type \code{GRanges}
}
\author{ Ellis Patrick }

% Add one or more standard keywords, see file 'KEYWORDS' in the
% R documentation directory.
sydSeq/man/pMim.Rd100777      0      0        4436 12405237430   7420 0\name{pMim}
\Rdversion{1}
\alias{pMim}
\title{Infer function of miRNA in RNA-Seg and miRNA-Seq data}
\description{
pMim jointly ranks miRNA and pathways with respect to a condition of interest.
}
\usage{
pMim(DataMi, DataG, classes, targets, pathways, Zmi=NULL, Zg=NULL, stat = 'cor', corP = 'Fisher', nInt = 2, randomTarg=FALSE, randomPath=FALSE, verbose == TRUE)}
\arguments{
\item{DataMi}{A matrix of normalised miRNA expression values.}
\item{DataG}{A matrix of normalised gene expression values.}
  \item{classes}{A vector containing two class labels for each sample. The names of which are the sample names.}
\item{targets}{A named list of the genes targeted by the named miRNA.}
\item{pathways}{A named list of the genes in each named pathway.}
\item{Zmi}{A vector of z-scores from a DE test on miRNA.}
\item{Zg}{A vector of z-scores from a DE test on genes.}
\item{stat}{A string describing which gene measure to use, 'cor' or 'de'.}
\item{nInt}{The minimum number of genes in each mir-pathway.}
\item{randomTarg}{Logical of whether gene target information should be randomised}
\item{randomPath}{Logical of whether pathway information should be randomised}
\item{verbose}{Should warnings be printed}
}
  \details{
Method ranks mir-pathways with respect to a condition of interest.
}
\value{
A list containing a formatted results matrix and the unformatted scores for each mir-pathway. Small scores provide evidence that the mir-pathway is interesting.
}
\author{ Ellis Patrick }
\examples{
### Generate Data
nMir = 10
nGene = 1000
nPath = 100
n = 20
DataMi = matrix(rnorm(nMir*n,1000,sqrt(1000)),nMir,n)
rownames(DataMi) = paste('mir',1:nMir,sep = '-')
DataG = matrix(rnorm(nGene*n,1000,sqrt(1000)),nGene,n)
rownames(DataG) = paste('gene',1:nGene,sep = '-')

targets = sample(rownames(DataG),nMir*100,replace = TRUE)
targets = split(targets,rep(rownames(DataMi),rep(100,nMir)))
targets = lapply(targets,unique)

pathways = sample(rownames(DataG),nPath*100,replace = TRUE)
pathways = split(pathways,rep(paste('path',1:nPath,sep = '-'),rep(100,nPath)))
pathways = lapply(pathways,unique)

classes = rep(0:1,c(10,10))
names(classes) = colnames(DataMi) = colnames(DataG) = paste('sample',1:n,sep = '-')

### Run pMim
output = pMim(DataMi, DataG, classes, targets, pathways)
}sydSeq/man/pMimHTML.Rd100777      0      0        2111 12405237430  10071 0\name{pMimHTML}
\Rdversion{1}
\alias{pMimHTML}
\title{Output HTML files for pMim}
\description{
This function outputs HTML files for the pMim function
}
\usage{
pMimHTML(output,dir = 'mirpathways',cutoff = 0.05,filename = 'pMimResults.html',GeneSymbol = NULL,outputHTML=FALSE)}
\arguments{
\item{output}{A list of output that was generated by pMim.}
\item{dir}{The directory you would like the information about the genes in each mirpathway to be written.}
  \item{cut-off}{Mir-pathways must have an intergrative score greater than this to be included in results.}
\item{filename}{The name you would like you html file to be called.}
\item{GeneSymbol}{A named vector that can transform the rownames of the Gene Data matrix to Gene Symbols.}
\item{outputHTML}{If TRUE this outputs a call to gvisTable.}
}
  \details{
This function write the results from a call to pMim to html files.
}
\value{
If outputHTML is TRUE this outputs the results from pMim as a call to gvisTable. Otherwise this simple writes the results to a html file.
}
\author{ Ellis Patrick }
\examples{}

sydSeq/man/simplifyGFeatures.Rd100777      0      0        1477 12405237430  12162 0\name{simplifyGFeatures}
\Rdversion{1}
\alias{simplifyGFeatures}
\title{Simplifies a GFeatures object}
\description{
Simplifies a GFeatures object removing regions that overlap multiple genes and splits exons into unique regions.
}
\usage{
simplifyGFeatures(gFeature)}
%- maybe also 'usage' for other objects documented here.
\arguments{
\item{gFeature}{An object of type \code{gFeature} }
}
\details{
.
}
\value{
An object of type \code{GRanges} }
\author{ Ellis Patrick }
\examples{
\dontrun{
hsapiensEnsembl <-makeTranscriptDbFromBiomart(biomart = "ensembl",dataset = "hsapiens_gene_ensembl")
saveFeatures(hsapiensEnsembl, file="hsapiensEnsembl.sqlite")
simpleEnsembl = simplifyGFeatures(hsapiensEnsembl)
}}
% Add one or more standard keywords, see file 'KEYWORDS' in the
% R documentation directory.
sydSeq/man/TshrinkPlus.Rd100777      0      0        4503 12405237430  10777 0\name{TshrinkPlus}
\Rdversion{1}
\alias{TshrinkPlus}
\title{Perform a moderated t-test using additional information}
\description{
TshrinkPlus is a performs a moderated t-test using additional information for normalised RNA-Seq counts.
}
\usage{
TshrinkPlus(Data, classes, meanPlus = NULL, variancePlus = NULL, verbose = TRUE,nn = 0.7)}
\arguments{
\item{Data}{A matrix of normlised counts with genes as row names and samples as column names.}
  \item{classes}{A vector containing two class labels for each sample. The names of which are the sample sample names.}
\item{meanPlus}{A vector of average gene expression from the additional data set named accordingly.}
\item{variancePlus}{A vector of gene variances from the additional data set named accordingly.}
\item{verbose}{Should warnings be printed}
\item{nn}{A parameter for kernal smoothing}
}
  \details{
Method performs a moderated t-test. The method uses additional information to fit the common variance curve if either \code{variancePlus} is not NULL or if both \code{variancePlus} and \code{meanPlus} are not NULL.
}
\value{
A list containg the pvalues, degrees of freedom and shrinkage parameter for the tests.
}
\author{ Ellis Patrick }
\examples{

set.seed(123)
# Generate true mean and variance
m = exp(rgamma(10000,5,1))
v = rnorm(10000)
v = pmax(m+1/40*m^2*exp(v),m)
# simulate "normalised counts"
Data = matrix(rnorm(10000*8,m,sqrt(v)),10000,8)
colnames(Data) = paste('s',1:8,sep ='')
rownames(Data) = paste('g',1:10000,sep ='')

# simulate additional data

Additional.Data = matrix(rnorm(10000*8,m,sqrt(pmax(m+1/40*m^2*exp(v),m))),10000,8)
rownames(Additional.Data) = paste('g',1:10000,sep ='')

# check mean vs variance
plot(apply(Data,1,mean),apply(Data,1,var),log = 'xy',xlab = 'Mean',ylab = 'Variance')
abline(0,1)

# calculate additional mean and variance vector
meanPlus = apply(Additional.Data,1,mean)
variancePlus = apply(Additional.Data,1,var)

# filter out genes with small means
variancePlus = variancePlus[meanPlus>5]
meanPlus = meanPlus[meanPlus>5]
Data = Data[rowMeans(Data)>10,]

# make class vector
classes = rep(1:2,c(4,4))
names(classes) = colnames(Data)

# test using tshrink and tshrinkplus
tshrinkplus = TshrinkPlus(Data = Data,classes = classes ,meanPlus,variancePlus)
tshrink = TshrinkPlus(Data,classes = classes)
}
sydSeq/NAMESPACE100777      0      0         433 12405237430   6644 0# Default NAMESPACE created by R
# Remove the previous line if you edit this file

# Export all names
exportPattern(".")

# Import all packages listed as Imports or Depends
import(
  limma,
  DESeq,
  edgeR,
  Biobase,
  GenomicFeatures,
  GenomicRanges,
googleVis
)
sydSeq/R/ 40777      0      0           0 12526447664   5566 5sydSeq/R/exClust.R100777      0      0        4776 12405237430   7436 0exClust = function(simpleGRanges,cut= 1.5){
	geneList <- IRanges::split(simpleGRanges,as.character(elementMetadata(simpleGRanges)$gene_id))
	simpleList <- IRanges::lapply(geneList,Clust,cut)

anno = eval(parse(text=(paste('c(',paste(paste('simpleList[[',1:length(simpleList),']]',sep = ''),collapse=','),')',sep = ''))))
anno
}


Clust =function(simpleGRanges, cut) {  
  
  major_group <- function(grouping) {
    group_id = names(table(grouping))[table(grouping) == 
                                        max(table(grouping))][1]
    return(names(grouping)[grouping == group_id])
  }
  counts = (elementMetadata(simpleGRanges)[, grep("counts", 
                                                  IRanges::colnames(elementMetadata(simpleGRanges)))])
  if (is.null(counts)) {
    cat("need counts matrix/n")
    break
  }
  IRanges::colnames(counts) = gsub("counts.", "", IRanges::colnames(counts))
  counts = IRanges::as.data.frame(counts)
  ncount = dim(counts)[2]
  counts = as.matrix(counts, ncol = ncount)
  NROW = nrow(counts)
  rownames(counts) = 1:nrow(counts)
  ir = ranges(simpleGRanges)
  width = width(ir)
  if (dim(counts)[1] > 2 & sum(rowSums(counts) > 5) > 1) {
    x.gs = colSums(counts)
    x.ge = rowSums(counts)
    x.g = sum(counts)
    MU = x.ge %*% t(x.gs)/x.g
    MU[MU < 0] = 0
    Z = (2 * sqrt(counts + 3/8) - 2 * sqrt(MU + 3/8))
    corZ = (Z[, ] %*% t(Z[, ]))/(dim(Z)[2] - 1)
    width = width[rowSums(counts) != 0]
    Z = Z[rowSums(counts) != 0, ]
    corZ = corZ[rowSums(counts) != 0, rowSums(counts) != 
                  0]
    counts = counts[rowSums(counts) != 0, ]
    probe_corr = corZ
    probe_corr[is.na(probe_corr)] = 0
    probe_corr_dist <- as.dist(1 - probe_corr)
    probe_tree <- hclust(probe_corr_dist, method = "complete")
    cut_trees = cutree(probe_tree, h = cut)
    countsW = counts/matrix(width, nrow(counts), ncol(counts))
    big = rowMeans(apply(countsW, 2, rank) * (rowSums(counts) > 
                                                0) * (width > 10))
    usetree = cut_trees
    tree = usetree
    bigsum2 = function(tree, big) {
      lu = length(unique(tree))
      SS = rep(0, lu)
      for (i in 1:lu) {
        bb = big[tree == i]
        bb = bb[bb != 0]
        SS[i] = mean(head(-sort(-bb), 2))
      }
      which(SS == max(SS,na.rm = TRUE))
    }
    bs = bigsum2(usetree, big)
    major_cluster <- names(usetree[usetree %in% bs])
    CC = simpleGRanges[(1:NROW) %in% major_cluster]
  }
  else CC = simpleGRanges
  CC
}sydSeq/R/geneCounts.R100777      0      0         645 12405237430  10070 0geneCounts = function(anno){
columns = grep('counts.',IRanges::colnames(elementMetadata(anno)))
counts = IRanges::as.data.frame(elementMetadata(anno)[,columns])
geneSplit = IRanges::split(counts,as.character(elementMetadata(anno)$gene_id))
sumCounts = IRanges::lapply(geneSplit,function(x)colSums(x))
counts = do.call('rbind',sumCounts)
colnames(counts) = gsub("counts.","",colnames(counts))
as.matrix(counts)
}
sydSeq/R/makeUI.R100777      0      0        2316 12405237430   7146 0makeUI = function( simpleGRanges, gFeature ){
require(Biobase)
cat("Loading data from database...\n")
 	GR <- exons( gFeature ,columns = c('tx_name','exon_id','exon_name','gene_id'))
geneID <- IRanges::unlist(elementMetadata(GR)$gene_id)
names(geneID) <- IRanges::unlist(elementMetadata(GR)$exon_name)
exonsG <- IRanges::unlist(exonsBy(gFeature, by=c("tx"),use.names=TRUE))
transcriptID <- substr(names(exonsG[,]),1,nchar(strsplit(names((exonsG)[,]),'.1')[[1]]))
names(transcriptID) <- elementMetadata(exonsG)$exon_name

tranSplit <- split(transcriptID,geneID[names(transcriptID)])
noTran <- unlist(lapply(tranSplit,function(x)length(unique(x))))

findUI = function(x){
useSimple <- seqnames((simpleGRanges)) == x
useGenomic <- seqnames(exonsG) == x
geneNames <- elementMetadata(simpleGRanges[useSimple])$gene_id
simpleGRanges[useSimple][countOverlaps(ranges(simpleGRanges[useSimple]),ranges(exonsG[useGenomic]))==noTran[geneNames]]
}

chrNames <- as.character(unique(seqnames(simpleGRanges)))
names(chrNames) <- chrNames
chrNames <- as.list(chrNames)
UIanno <- lapply(chrNames,findUI)
eval(parse(text=(paste('c(',paste(paste('UIanno[[',1:length(UIanno),']]',sep = ''),collapse=','),')',sep = ''))))
}

sydSeq/R/mergeCounts.R100777      0      0         467 12405237430  10253 0mergeCounts = function(c,GRanges){
	counts <- c
	merged <- GRanges(seqnames = seqnames(GRanges),
		ranges = ranges(GRanges),
		strand = strand(GRanges), gene_id = GRanges@elementMetadata$gene_id,
		counts = as.data.frame(counts))
      	names(merged) <- GRanges@elementMetadata$gene_id
      	merged
}
sydSeq/R/pMim.R100777      0      0       11471 12552232511   6715 0pMim = function (DataMi, DataG, classes, targets, pathways, Zmi = NULL, 
                 Zg = NULL, stat = "cor", corP = "Fisher", nInt = 2, randomTarg = FALSE, 
                 randomPath = FALSE, verbose = TRUE) 
{
  if (verbose == TRUE) {
    if (sum(colnames(DataG) != names(classes)) > 0 | is.null(names(classes))) {
      cat("colnames of Gene Data do not match names of classes")
    }
  }
  if (verbose == TRUE) {
    if (sum(colnames(DataG) != names(classes)) > 0 | is.null(names(classes))) {
      cat("colnames of miRNA Data do not match names of classes")
    }
  }
  if (verbose == TRUE) {
    if (sum(rowMeans(DataG) < 5) > 0) {
      cat("You might want to consider filtering out genes with small average expression from your data.\n")
    }
  }
  if (verbose == TRUE) {
    if (sum(rowMeans(DataMi) < 5) > 0) {
      cat("You might want to consider filtering out miRNA with small average expression from your data.\n")
    }
  }
  y = names(classes)
  X = t(DataG)
  Y = t(DataMi)
 
    corYX = cor(Y[y, ], X[y, ])
    n = length(y)
    if (corP == "estimate") {
      TR = corYX
      TR[] = ((rank(corYX) - 0.5)/length(corYX))
      TR = qnorm(TR)
    }
    else {
      if (n > 3) {
        TR <- sqrt(n - 3) * atanh(corYX)
      }
      else {
        cat("Need more replicates or use cor=\"estimate\"")
        break
      }
    }
  design = model.matrix(~classes)
  
  if (is.null(Zmi)) {
    require(limma)
    design = model.matrix(~classes)
    Dat = DataMi[, names(classes)]
    fit = lmFit(Dat, design)
    ordinary.t <- fit$coef/fit$stdev.unscaled/pmax(pmax(fit$sigma, 
                                                        sqrt(rowMeans(Dat))), 1)
    Tmi = ordinary.t[, 2]
    Zmi = qnorm(pt((Tmi), fit$df.res))
  }
  
  if (is.null(Zg)) {
    require(limma)
    design = model.matrix(~classes)
    Dat = DataG[, names(classes)]
    fit = lmFit(Dat, design)
    ordinary.t <- fit$coef/fit$stdev.unscaled/pmax(pmax(fit$sigma, 
                                                        sqrt(rowMeans(Dat))), 1)
    Tg = ordinary.t[, 2]
    Zg = qnorm(pt((Tg), fit$df.res))
  }
  
  mi = colnames(Y)
  msgr = intersect(unlist(pathways), unlist(targets))
  mapMat = matrix(0, length(mi), length(msgr))
  rownames(mapMat) = mi
  colnames(mapMat) = msgr
  for (i in rownames(mapMat)) {
    if (randomTarg == TRUE) {
      mapMat[i, sample(colnames(mapMat), sum(msgr %in% 
                                               targets[[i]]))] = 1
    }
    else {
      mapMat[i, msgr %in% targets[[i]]] = 1
    }
  }
  use = names(which(rowSums(mapMat) > 0))
  mapMat = mapMat[use, ]
  mi = use
  pathname = names(pathways)
  pathMat = matrix(0, length(pathname), length(msgr))
  rownames(pathMat) = pathname
  colnames(pathMat) = msgr
  for (j in 1:length(pathname)) {
    if (randomPath == TRUE) {
      pathMat[j, sample(colnames(mapMat), sum(msgr %in% 
                                                pathways[[j]]))] = 1
    }
    else {
      pathMat[j, msgr %in% pathways[[j]]] = 1
    }
  }
  N = mapMat %*% t(pathMat)
  test = which(N >= nInt, 2)
  G = matrix(NA, length(mi), length(pathname))
  rownames(G) = mi
  colnames(G) = pathname
  Tmat = corYX
  Tmat[] = NA
  for (k in 1:nrow(test)) {
    i = mi[test[k, 1]]
    j = test[k, 2]
    pathway = pathMat[j, ]
    tr = TR[i, ]
    binding = mapMat[i, ]
    if (tolower(stat) == "de") {
      t = Zg[names(TR[i, ])]*sign(Zmi[i])
    }
    if (tolower(stat) == "cor") {
      t = TR[i, ]
    }
    Tmat[i, ] = t
    g = sum(t[names(which(binding * pathway == 1))])/sqrt(sum(binding * 
                                                                pathway))
    G[i, j] = pnorm(g)
  }
  G = t(apply(G, 1, function(x) {
    x[!is.na(x)] = p.adjust(x[!is.na(x)], "fdr")
    return(x)
  }))
  Score = pchisq(-2 * (log(1 - G) + log(1 - matrix(pnorm(-abs(Zmi))[rownames(G)], 
                                                   dim(G)[1], dim(G)[2], byrow = FALSE))), 4)
  pval = Score
  updown = matrix(sign(Zmi)[rownames(G)], dim(G)[1], dim(G)[2], 
                  byrow = FALSE)
  updown[updown == 1] = "Up"
  updown[updown == -1] = "Down"
  use = which(!is.na(pval), 2)
  res = cbind(rownames(use) ,  colnames(pval)[use[,2]] , updown[use] , pval[use])
  res = res[order(as.numeric(res[, 4])), ]
  
  colnames(res) = c("miRNA", "Pathways", "direction", "Score")
  
  direction = NULL
  if(!(is.null(Zmi)&is.null(Zg))){
    class1 = gsub('classes','',colnames(design)[2])
    class2 = classes[!classes%in%class1][1]
    direction = paste('"Up" is',class1,  'expression greater than', class2, 'expression in miRNA')
    cat(direction,'\n')
  }
  list(Results = res, Scores = Score, Zg = Zg, Zmi = Zmi, cor = corYX, corTransform = TR,Direction = direction,targets = targets, pathways = pathways)
}sydSeq/R/pMimHTML.R100777      0      0        3636 12552234101   7363 0pMimHTML = function(output,dir = 'mirpathways',cutoff = 0.05,filename = 'pMimResults.html',GeneSymbol = NULL,outputHTML=FALSE){
  require(googleVis)
  results = output$Results
  results = results[as.numeric(results[,'Score']) <= cutoff,]
  results = data.frame(results, miRNA.DE.pvalue = 2*pnorm(-abs(output$Zmi[results[,'miRNA']])),genes = rep(NA,dim(results)[1]),stringsAsFactors = FALSE)
  results[,'Score'] = as.numeric(results[,'Score'])
  results = results[,c(1,2,6,3,5,4)]
  colnames(results) = c('microRNA','Pathway','Genes','Direction','microRNA DE p-value','Intergrative Score')
  
  dir.create(dir)
  gene.files = NULL
  for(i in 1:dim(results)[1]){
    miRNA = as.character(results[i,1])
    path = as.character(results[i,2])
    genes = intersect(output$targets[[miRNA]],output$pathways[[path]])
    Correlation = output$cor[miRNA,genes]
    corP = pnorm(output$corTransform[miRNA,genes])
    DE = pnorm(output$Zg[genes]*sign(output$Zmi[miRNA]))
    if(!is.null(GeneSymbol))genes = GeneSymbol[genes]
    res = data.frame(Genes = genes,'One Sided DE P-value' = DE,Correlation = Correlation,'One Sided Correlation P-value' = corP)
    form = as.list(rep('#.###',4))
    names(form) = colnames(res)
    file = paste(miRNA,'-',path,'.html',sep = '')
    file = gsub('-','_',file)
    file = paste(dir,file,sep= '/')
    gene.files[i] = file
    res = res[order(res[,4]),]
    print(gvisTable(res,format = form,options = list(width = '60em')),type = 'html',file = file)
  }
  
  
  results[,'Genes'] = paste('<a href=','"',gene.files,'" target="_blank">','genes','</a>',sep="")
  
  library(googleVis)
  form = as.list(rep('#.###',6))
  names(form) = colnames(results)
  results = as.data.frame(results)
  print(gvisTable(results,format = form,options = list(width = '80em')),type = 'html',file = filename)
  
  if(outputHTML == TRUE)gvisTable(results,format = form,options = list(width = '100em'))
}
sydSeq/R/simplifyGFeatures.R100777      0      0        3051 12405237430  11432 0simplifyGFeatures = function( gFeature ){
require(GenomicFeatures)
require(GenomicRanges)
	cat("Loading data from database...\n")
 	GR <- exons( gFeature ,columns = c('exon_id','gene_id'))
	cat("Simplifying annotation...\n")
	ranges <- ranges(GR)
	names(ranges) <- IRanges::unlist(elementMetadata(GR)$gene_id)
 	chrList <- IRanges::split(ranges,as.character(seqnames(GR)))
	simpleList <- IRanges::lapply(chrList,simplifyIRanges)
	simpleAnno <- NULL
	for(i in 1:length(simpleList)){
		simpleAnno = rbind(cbind(ensembl_gene_id=names(simpleList[[i]]), start = as.numeric(simpleList[[i]]@start), end = as.numeric(simpleList[[i]]@start + simpleList[[i]]@width -1),chr = names(simpleList[i]), strand = '*'),simpleAnno)
		}
	simpleRanges <- GRanges(seqnames = simpleAnno[,'chr'],
		ranges = IRanges(start = as.numeric(simpleAnno[,'start']),
		end = as.numeric(simpleAnno[,'end'])),
		strand = simpleAnno[,'strand'], gene_id = simpleAnno[,'ensembl_gene_id'])
      names(simpleRanges) = simpleAnno[,'ensembl_gene_id']
	simpleRanges
	}


simplifyIRanges = function(iRange){
disjointRanges <- IRanges::disjoin(iRange)
un <- unique(names(iRange))
nam <- rep('non',length(disjointRanges))
for(i in 1:length(un)){
use = which(!is.na(IRanges::findOverlaps(disjointRanges,iRange[which(names(iRange)%in%un[i])],select = 'first')))
nam[use][nam[use]!= 'non'] <- 'multi'
nam[use][nam[use]== 'non'] <- as.character(un[i])
}
names(disjointRanges) <- nam
simpleRanges <- disjointRanges[names(disjointRanges)!='non'&names(disjointRanges)!='multi']
simpleRanges
}
sydSeq/R/TshrinkPlus.R100777      0      0        7774 12405237430  10276 0TshrinkPlus = function(Data, classes, meanPlus = NULL, variancePlus = NULL, verbose = TRUE,nn = 0.7){

if(verbose == TRUE){
  if(sum(colnames(Data) != names(classes))>0 |is.null(names(classes))){
    cat('colnames of Data do not match names of classes')
  }
  }

if(verbose == TRUE){
  if(sum(rowMeans(Data)<5)>0){
    cat('You might want to consider filtering out genes with small average expression from your data')
  }
  }


if(verbose == TRUE){
  if(sum((meanPlus)<5)>0&sum((meanPlus)>100)>0){
    cat('You might want to consider filtering out genes with small average expression from your additional data')
  }
}

int = rownames(Data)
if(!is.null(meanPlus)){
varStandard = StandardiseVariance(meanPlus,variancePlus)
int = intersect(rownames(Data), names(which(!is.na(varStandard ))))
}

if(is.null(meanPlus)&!is.null(variancePlus)){
  varStandard = variancePlus
  int = intersect(rownames(Data), names(which(!is.na(varStandard ))))
}

design = matrix(0,length(classes),length(unique(classes)))
design[,1] = 1
for(i in 2:length(unique(classes))){
design[classes == unique(classes)[i],i] = 1
}

if(length(unique(classes))>2){cat('Please use only two unique class labels... pretty please');break}

udes = unique(design)
designUniq = design
for(i in 1:nrow(udes)){
  designUniq[,i] = rowSums(design==matrix(udes[i,],dim(design)[1],dim(design)[2],byrow = TRUE))==ncol(design)
}
colUse = colSums(designUniq)>1
if(sum(colUse)==0)stop("Need more replication or simpler design")

if(is.null(meanPlus)){
  
  PREDVAR = as.matrix(Data)%*%designUniq
  VARIANCE = as.matrix(Data)%*%designUniq
  for(i in which(colUse)){
    rowM = rowMeans(Data[,designUniq[,i]==1])
    variance = apply(Data[,designUniq[,i]==1],1,var)
    VARIANCE[,i] = variance
    avgExpr = rowM
    predVar = avgExpr
    u = avgExpr>0
    fit = locfit((variance[u])~lp(log(avgExpr[u]),nn = nn),family = 'gamma')
    predVar4 = (fitted( fit )) 
    predVar[u] = predVar4
    predVar1 = pmax(predVar,avgExpr)
    PREDVAR[,i] = predVar1
  }
      
}


if(!is.null(meanPlus)){
  
  
  PREDVAR = as.matrix(Data[int,])%*%designUniq
  VARIANCE = as.matrix(Data[int,])%*%designUniq
  for(i in which(colUse)){
    rowM = rowMeans(Data[int,designUniq[,i]==1])
    variance = apply(Data[int,designUniq[,i]==1],1,var)
    VARIANCE[,i] = variance
    avgExpr = rowM
    predVar = avgExpr
    u = avgExpr>0
     qV = varStandard [int][u]
  
    fit = locfit((variance[u])~lp(log(avgExpr[u]),qV,nn = nn,scale = TRUE),maxk=300,family = 'gamma')
    predVar4 = (fitted( fit ))
    predVar[u] = predVar4
    predVar1 = pmax(predVar,avgExpr)
    PREDVAR[,i] = predVar1
  }
  }

u = names(which(rowMeans(Data[int,])>50))
N = matrix(colSums(designUniq),length(u),dim(designUniq)[2],byrow = TRUE)
lambda = 2/mean((N-1)*(VARIANCE[u,]/PREDVAR[u,]-1)^2,na.rm =TRUE)
N = matrix(colSums(designUniq),dim(VARIANCE)[1],2,byrow = TRUE)

VAR =  ((rowSums(PREDVAR*(N-1))/(sum(designUniq)-2))*lambda+ (rowSums(VARIANCE*(N-1))/(sum(designUniq)-2))*(1-lambda))



n2 = Inf
n1 = sum(designUniq)-2
df = ((1)*lambda+ ((rowSums(VARIANCE*(N-1))/(sum(designUniq)-2))/(rowSums(PREDVAR*(N-1))/(sum(designUniq)-2)))*(1-lambda))^2/(((rowSums(VARIANCE*(N-1))/(sum(designUniq)-2))/(rowSums(PREDVAR*(N-1))/(sum(designUniq)-2)))^2*(1-lambda)^2/n1+ lambda^2/n2)

df = df[!is.na(df)]
df = mean(df)

T = (rowSums(Data[int,classes == unique(classes)[1]]) - rowSums(Data[int,classes == unique(classes)[2]]))^2/(VAR*length(classes))
P = 2*pt(sqrt(T),df,lower.tail = FALSE)
P = P[!is.na(P)]

PP = rep(NA,dim(Data)[1])
names(PP) == rownames(Data)
PP[int] = P
pvalue = PP
list(pvalue = pvalue,df = df,lambda = lambda )

}

StandardiseVariance = function(Mean,Variance,nn = 0.7){
  use = Mean>0
  if(max(Mean)>30&min(Mean)>=0)Mean = log(Mean)
  # use = Mean>0#&Variance>0
  VARS = Variance[use]
  AVG = Mean[use]
  fit = locfit(VARS~lp(AVG,nn = nn),family = 'gamma')
  VARS = log(VARS/fitted(fit))
  Variance[!use] = NA
  Variance[use] = VARS
  Variance
}