inst/doc/Intermediate.R

## ---- echo = FALSE, message = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
library(myTAI)
options(width = 750)
knitr::opts_chunk$set(
  comment = "#>",
  error = FALSE,
  tidy = FALSE)

## ----eval=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  library(myTAI)
#  
#  # load an example PhyloExpressionSet stored in the myTAI package
#  data(PhyloExpressionSetExample)
#  
#  # look at the standardized data set format
#  head(PhyloExpressionSetExample, 3)

## ----eval=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # TAI profile of correctly assigned Phylostrata
#  TAI(PhyloExpressionSetExample)

## ---- fig.width= 7, fig.height= 5,eval=TRUE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
data(PhyloExpressionSetExample)
# Visualize the TAI profile of correctly assigned Phylostrata
PlotSignature( ExpressionSet = PhyloExpressionSetExample, 
               p.value       = FALSE ) 

## ----eval=TRUE, echo=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# randomly permute the Phylostratum assignment of all genes 
randomPhyloExpressionSetExample <- data.frame( sample(PhyloExpressionSetExample[ , "Phylostratum"]), 
                                               PhyloExpressionSetExample[, 2:9] )

# TAI profile based on randomly assigned Phylostrata
TAI(randomPhyloExpressionSetExample) 

## ---- fig.width= 7, fig.height= 5,eval=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # Visualize the TAI profile based on randomly assigned Phylostrata
#  PlotSignature( ExpressionSet =  randomPhyloExpressionSetExample,
#                 p.value       = FALSE )

## ----eval=TRUE, echo=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Variance of the TAI profile based on correctly assigned Phylostrata
var(TAI(PhyloExpressionSetExample)) 

## ----eval=TRUE, echo=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Variance of the TAI profile based on randomly assigned Phylostrata
var(TAI(randomPhyloExpressionSetExample)) 

## ----eval=TRUE, echo=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Compute 1000 TAI profiles based on randomly assigned Phylostrata
randomTAIs <- bootMatrix( ExpressionSet = PhyloExpressionSetExample, 
                          permutations  = 1000 )

head(randomTAIs)

## ---- fig.width= 7, fig.height= 5,eval=TRUE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# compute the variance of the random TAI profile for each row
variance_vector <- apply(randomTAIs, 1 , var)

# and visualize the distribution of variances
hist(variance_vector, breaks = 100)

## ---- fig.width= 7, fig.height= 5,eval=TRUE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# variance of the TAI profile based on correctly assigned Phylostrata
var_real <- var(TAI(PhyloExpressionSetExample)) 

# visualize the distribution of variances
hist( x      = c(variance_vector,var_real), 
      breaks = 100, 
      xlab   = "variance", 
      main   = "Histogram of variance_vector" )

# and plot a red line at the position where we can find the 
# real variance
abline(v = var_real, lwd = 5, col = "red")

## ---- fig.width= 7, fig.height= 5,eval=TRUE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# install.packages("fitdistrplus")

# plot a Cullen and Frey graph
fitdistrplus::descdist(variance_vector)

## ---- fig.width= 7, fig.height= 5,eval=TRUE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# estimate the parameters: shape and rate using 'moment matching estimation'
gamma_MME <- fitdistrplus::fitdist(variance_vector,distr = "gamma", method = "mme")
# estimate shape:
shape <- gamma_MME$estimate[1]
# estimate the rate:
rate <- gamma_MME$estimate[2]

# define an expression written as function as input for the curve() function
gamma_distr <- function(x){ return(dgamma(x = x, shape = shape, rate = rate)) }

# plot the density function and the histogram of variance_vector
curve( expr = gamma_distr,
       xlim = c(min(variance_vector),max(c(variance_vector,var_real))),
       col  = "steelblue",
       lwd  = 5,
       xlab = "Variances", 
       ylab = "Frequency" )

# plot the histogram of variance_vector
histogram <- hist(variance_vector,prob = TRUE,add = TRUE, breaks = 100)
rug(variance_vector)

# plot a red line at the position where we can find the real variance
abline(v = var_real, lwd = 5, col = "red")

## ----eval=TRUE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# p-value of var_real
pgamma(var_real, shape = shape,rate = rate, lower.tail = FALSE)

## ----eval=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # Perform the FlatLineTest
#  FlatLineTest( ExpressionSet = PhyloExpressionSetExample,
#                permutations  = 1000 )

## ---- fig.width= 7, fig.height= 3,eval=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  
#  # perform the FlatLineTest and investigate the goodness of the test statistic
#  FlatLineTest( ExpressionSet = PhyloExpressionSetExample,
#                permutations  = 1000,
#                plotHistogram = TRUE )
#  

## ----eval=TRUE, echo=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
library(myTAI)

# load an example PhyloExpressionSet stored in the myTAI package
data(PhyloExpressionSetExample) 

# look at the standardized data set format
head(PhyloExpressionSetExample, 3) 

## ----eval=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # TAI profile of correctly assigned Phylostrata
#  TAI(PhyloExpressionSetExample)

## ---- fig.width= 7, fig.height= 5,eval=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # visualize the TAI profile of correctly assigned Phylostrata
#  PlotSignature( ExpressionSet = PhyloExpressionSetExample,
#                 p.value       = FALSE )

## ---- fig.width= 7, fig.height= 5,eval=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # Perform the Reductive Hourglass Test
#  ReductiveHourglassTest( ExpressionSet = PhyloExpressionSetExample,
#                          modules       = list(early = 1:2, mid = 3:5, late = 6:7),
#                          lillie.test   = TRUE )
#  

## ---- fig.width= 7, fig.height= 6,eval=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # perform the Reductive Hourglass Test and plot the test statistic
#  ReductiveHourglassTest( ExpressionSet = PhyloExpressionSetExample,
#                          modules       = list(early = 1:2, mid = 3:5, late = 6:7),
#                          plotHistogram = TRUE,
#                          lillie.test   = TRUE )
#  

## ----eval=TRUE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
library(myTAI)

# load an example PhyloExpressionSet stored in the myTAI package
data(PhyloExpressionSetExample) 

# look at the standardized data set format
head(PhyloExpressionSetExample, 3) 

## ----eval=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # TAI profile of correctly assigned Phylostrata
#  TAI(PhyloExpressionSetExample)

## ---- fig.width= 7, fig.height= 5,eval=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # Visualize the TAI profile of correctly assigned Phylostrata
#  PlotSignature( ExpressionSet = PhyloExpressionSetExample,
#                 p.value       = FALSE )
#  

## ---- fig.width= 7, fig.height= 5,eval=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  
#  # Perform the Reductive Early Conservation Test
#  EarlyConservationTest( ExpressionSet = PhyloExpressionSetExample,
#                         modules       = list(early = 1:2, mid = 3:5, late = 6:7),
#                         lillie.test   = TRUE )
#  

## ---- fig.width= 7, fig.height= 6,eval=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # perform the Reductive Early Conservation Test and plot the test statistic
#  EarlyConservationTest( ExpressionSet = PhyloExpressionSetExample,
#                         modules       = list(early = 1:2, mid = 3:5, late = 6:7),
#                         plotHistogram = TRUE,
#                         lillie.test   = TRUE )
#  

## ---- fig.width= 7, fig.height= 5-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
library(myTAI)

data(PhyloExpressionSetExample)

# a simple example is to transform the gene expression levels of a given PhyloExpressionSet
# using a sqrt or log2 transformation

PES.sqrt <- tf(PhyloExpressionSetExample, sqrt)

head(PES.sqrt)

## ---- fig.width= 7, fig.height= 5, eval=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  
#  PES.log2 <- tf(PhyloExpressionSetExample, log2)
#  
#  head(PES.log2)

## ---- fig.width= 7, fig.height= 5,eval=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # in case a given PhyloExpressionSet already stores gene expression levels
#  # that are log2 transformed and need to be re-transformed to absolute
#  # expression levels, to perform subsequent phylotranscriptomics analyses
#  # (that are defined for absolute expression levels),
#  # one can re-transform a PhyloExpressionSet like this:
#  
#  PES.absolute <- tf(PES.log2 , function(x) 2^x)
#  
#  # which should be the same as  PhyloExpressionSetExample :
#  head(PhyloExpressionSetExample)
#  head(PES.absolute)

## ---- fig.width= 7, fig.height= 5, eval=TRUE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
data(PhyloExpressionSetExample)
# plotting the TAI using log2 transformed expression levels
# and performing the Flat Line Test to obtain the p-value
PlotSignature( ExpressionSet = tf(PhyloExpressionSetExample, log2),
               TestStatistic = "FlatLineTest",
               xlab          = "Ontogeny", 
               ylab          = "TAI" )



## ---- fig.width= 7, fig.height= 5, eval=TRUE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
data(PhyloExpressionSetExample)
# plotting the TAI using sqrt transformed expression levels
# and performing the Flat Line Test to obtain the p-value
PlotSignature( ExpressionSet = tf(PhyloExpressionSetExample, sqrt),
               TestStatistic = "FlatLineTest",
               xlab          = "Ontogeny", 
               ylab          = "TAI" )

## ---- fig.width= 7, fig.height= 5, eval=TRUE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
data(DivergenceExpressionSetExample)
# plotting the TDI using log2 transformed expression levels
# and performing the Flat Line Test to obtain the p-value
PlotSignature( ExpressionSet = tf(DivergenceExpressionSetExample, log2),
               TestStatistic = "FlatLineTest",
               xlab          = "Ontogeny", 
               ylab          = "TDI" )

## ---- fig.width= 7, fig.height= 5, eval=TRUE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
data(DivergenceExpressionSetExample)
# plotting the TDI using sqrt transformed expression levels
# and performing the Flat Line Test to obtain the p-value
PlotSignature( ExpressionSet = tf(DivergenceExpressionSetExample, sqrt),
               TestStatistic = "FlatLineTest",
               xlab          = "Ontogeny", 
               ylab          = "TDI" )

Try the myTAI package in your browser

Any scripts or data that you put into this service are public.

myTAI documentation built on July 4, 2017, 9:09 a.m.