Intermediate.R
In myTAI: Evolutionary Transcriptomics

## ---- 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=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  data(PhyloExpressionSetExample)
#  # Visualize the TAI profile of correctly assigned Phylostrata
#  PlotSignature( ExpressionSet = PhyloExpressionSetExample,
#                 p.value       = FALSE )

## ----eval=FALSE, 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=FALSE, echo=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # Variance of the TAI profile based on correctly assigned Phylostrata
#  var(TAI(PhyloExpressionSetExample))

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

## ----eval=FALSE, 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=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # 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=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # 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=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # install.packages("fitdistrplus")
#  
#  # plot a Cullen and Frey graph
#  fitdistrplus::descdist(variance_vector)

## ---- fig.width= 7, fig.height= 5,eval=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # 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=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # 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=FALSE, 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=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 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 )
#  

## ----eval=FALSE, 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 Reverse Hourglass Test
#  ReverseHourglassTest( 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 Reverse Hourglass Test and plot the test statistic
#  ReverseHourglassTest( 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=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  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=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  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=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  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=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  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" )