Package {DiPALM}

DiPALM - Differential Pattern Analysis via Linear Modeling

Description

Individual gene expression patterns are encoded into a series of eigenvector patterns (WGCNA package). Using the framework of linear model-based differential expression comparisons (limma package), time-course expression patterns for genes in different conditions are compared and analyzed for significant pattern changes.

Author(s)

Ryan C. Sartor & Kathleen Greenham

Maintainer: Ryan C. Sartor <sartorry@gmail.com>

See Also

limma package blockwiseModules from the WGCNA package

Determine p-values based on permutation test results

Description

Used to calculate the false discovery rate associated with a specified cutoff. It uses a vector of actual test scores and a vector of permuted scores. The result can be considered a corrected p-value or q-value.

Usage

AdjustPvalue(tVal, tVec, pVec)

Arguments

Details

This function assumes all test result values above the given threshold are true positives (TP) and all permuted result values above the threshold are false positives (FP). It calculates FDR as (FP)/(TP). When applied to individual test results, this works to estimate an adjusted p-value (or q-value) for that result.

Value

A numeric value representing the FDR associated with the given cutoff (tVal).

Author(s)

Ryan C. Sartor

Examples

# Returns a vector of adjusted p-values from a vector of test results
data("testData")
AdjPval_kMEs<-sapply(testData$testResults[1:100],function(x) 
AdjustPvalue(tVal = x, tVec = testData$testResults, pVec = testData$permutedResults))

Generate a blank plot window

Description

Creates a blank plotting window of desired size.

Usage

BlankPlot(xrng = c(0, 1), yrng = c(0, 1), Main = "", xlab = "", ylab = "", ...)

Arguments

Details

This function is used to initiate a blank plotting window in R. xrng and yrng are used to initiate the size of the plotting region.

Value

Generates a blank plot

Author(s)

Ryan C. Sartor

See Also

plot

Examples

BlankPlot(xrng=c(0,10),yrng=c(0,10),Main="Test",xlab="xlab",ylab="ylab")

Fit linear models to time-course expression patterns using the limma package

Description

This function fits a linear model (using lmFit form the limma package) to each set of kME values returned from BuildModMembership.

Usage

BuildLimmaLM(dataMat, designMat, contrastStr)

Arguments

Details

This function uses the limma package functionality to carry out an analysis of differential patterning. Asking if a gene displays a significantly different expression pattern across time when looking at one condition compared to another. The concept of differential patterning is similar to differential expression. However, the expression level of most genes fluctuates significantly throughout the day. Most of this is due to diurnal or circadian rhythms. These regular changes in expression severely complicate differential expression analysis and give rise to the need for examining differential patterning. By "encoding" the expression pattern of each gene into a vector of module membership (kME) values using the BuildModMembership function, the DiPALM analysis provides a way to use the limma method to examine differential patterning.

Value

This function returns an MArrayLM-class object. This is a custom class from the limma package. It is an object that contains all information and results of individual linear models built for each gene in the input matrix.

Author(s)

Ryan C. Sartor

See Also

BuildLimmaLM , lmFit , makeContrasts , contrasts.fit , eBayes

Examples

data("testData")
# First calculate the module memberships
kMEsList<-BuildModMembership(MeMat=testData$moduleEigengenes, TCsLst=testData$timeCourseList)

# Then build the limma models
LimmaMods<-lapply(kMEsList, function(x) 
BuildLimmaLM(dataMat = x, designMat = testData$modelDesign, contrastStr = testData$modelContrast))

Calculate module membership (kME) for all genes in a time-course data set

Description

This function calculates module membership (kME) between each gene's expression vector throughout a time-course and the module eigengenes of the data set (from the WGCNA package, using blockwiseModules). Module membership is defined as the Pearson correlation between the expression of a single gene and the expression of a single eigengene.

Usage

BuildModMembership(MeMat, TCsLst)

Arguments

Details

The matrix of eigengene vectors, MeMat (see blockwiseModules from the WGCNA package) is used as a convenient way to summarize all the unique expression patterns (across time) that are present in a time-course data set. Any large gene expression data set from a complex organism will contain tens of thousands of genes. However, extensive coordinated regulation is present in all biological systems. Because of this co-regulation, usually only dozens of distinct expression patterns exist and most genes can be categorized by one (or a few) of them. By using a relatively small set of descriptive expression vectors (eigengenes), we can "encode" the pattern of expression for any single gene into a length N vector where N is the number of eigengenes and the values of the vector are the Pearson correlation of the gene to each eigengene (kME). These kME values can then be compared using the limma package to determine differential patterning similar to how expression values are compared to determine differential expression. This function is the first step to calculate all kMEs for each gene in a list of time-course matrices relative to the set of eigengenes for the whole data set.

Value

This function returns a list of kME matrices. The list contains a separate element for each eigengene. Each element is a matrix with a row for each gene and a column for each member of the input TCsLst.

Author(s)

Ryan C. Sartor

See Also

blockwiseModules

Examples

data("testData")
kMEsList<-BuildModMembership(MeMat=testData$moduleEigengenes, TCsLst=testData$timeCourseList)

Plot multiple time-course expression patterns for one gene

Description

A plotting function that generates a line plot for a single gene in multiple time-course expression sets. Each time-course is overlaid on the same plot.

Usage

PlotTCs(TClst, tgene, main = "", scale = TRUE, xlab = "Time", 
    ylab = "", xAxsLabs = colnames(TClst[[1]]), ledgeX = "top", 
    colAdj = 0.4, tcols = c("red", "red", "blue", "blue"), tltys = c(1, 1, 2, 2))

Arguments

Details

This function is used to view the expression pattern of a single gene across multiple time-course samples.

Value

The function generates a plot with multiple lines.

Author(s)

Ryan C. Sartor

Examples

data(testData)
PlotTCs(TClst = testData$timeCourseList, tgene = "BraA05g36370R" ,
scale = TRUE,tcols = c("red","red","blue","blue"), tltys = c(1,2,1,2), ledgeX = "topleft")

Plot multiple averaged time-course expression patterns for a group of genes

Description

A plotting function that generates a ribbon plot representing the summarized expression of a set of genes. A separate ribbon is generated for each time-course. The ribbon represents + & - 1 standard deviation from the mean at each point.

Usage

PlotTCsRibbon(TClst, tgenes, main = "", xlab = "Time", ylab = "",
 xAxsLabs = colnames(TClst[[1]]),scale = TRUE, alpha = 0.1, tcols, tltys,
 splits = rep(1, length(TClst)))

Arguments

Details

This can be used to plot the summarized expression pattern of an entire cluster of genes and compare between multiple time courses. Line plots are generated representing the mean expression value at each time point for the group. Ribbons are plotted along the lines in semi-transparent colors. The ribbons represent + and - one standard deviation from the mean expression level of the group. #

Value

This function generates multiple line plots representing the average expression of a group of genes with ribbons representing standard deviation.

Author(s)

Ryan C. Sartor

Examples

data("testData")

# Calculate the module memberships
kMEsList<-BuildModMembership(MeMat=testData$moduleEigengenes, TCsLst=testData$timeCourseList)

# Build the limma models
LimmaMods<-lapply(kMEsList, function(x) 
BuildLimmaLM(dataMat = x, designMat = testData$modelDesign, contrastStr = testData$modelContrast))

# Get the adjusted Pvalues for each gene
AdjPval_kMEs<-sapply(testData$testResults[rownames(LimmaMods[[1]])],function(x) 
AdjustPvalue(tVal = x, tVec = testData$testResults, pVec = testData$permutedResults))

# Cluster the genes base on the limma results
LimmaSums<-do.call(cbind,lapply(LimmaMods,function(x) x$t))
LimmaModsSig<-LimmaSums[names(AdjPval_kMEs)[which(AdjPval_kMEs<0.01)],]
patternCor<-cor(t(LimmaModsSig))
patternTree<-hclust(as.dist(1-patternCor),method = "complete")
patternClusters<-cutree(tree = patternTree, k = 25)
patternClusters<-tapply(rownames(patternCor),INDEX = patternClusters,function(x) x)

#Replicates are grouped together (same color + lty values)
PlotTCsRibbon(TClst = testData$timeCourseList, main="Drought Time-Course", 
xAxsLabs = c(seq(1,23,4),seq(1,23,4)), xlab="ZT Time (hours)", tgenes = patternClusters[[15]],
scale = TRUE, tcols = c("red","red","blue","blue"), tltys = c(1,1,1,1))

#Replicates are in different groups (same color but different lty values)
PlotTCsRibbon(TClst = testData$timeCourseList, main="Drought Time-Course",
xAxsLabs = c(seq(1,23,4),seq(1,23,4)), xlab="ZT Time (hours)", tgenes = patternClusters[[15]],
scale = TRUE, tcols = c("red","red","blue","blue"), tltys = c(1,2,1,2))

#Replicates are represented on differnt plots (using the splits argument)
PlotTCsRibbon(TClst = testData$timeCourseList, main=c("Rep1","Rep2"), 
xAxsLabs = c(seq(1,23,4),seq(1,23,4)), xlab="ZT Time (hours)", tgenes = patternClusters[[15]],
scale = TRUE, tcols = c("red","red","blue","blue"), tltys = c(1,2,1,2), splits = c(1,2,1,2))

Substitute values in for the time, treatment, and replicate portion of a column format.

Description

Given a column format as a string this function substitutes each instance of '[TREATMENT]', '[REPLICATE]' and '[TIME]' with the corresponding given values.

Usage

col_format_sub(
	col.format = "[TREATMENT].[REPLICATE].[TIME]",
	tr = ".*",
	re = ".*",
	ti = ".*",
	escape = TRUE
)

Arguments

Value

A new string with the placeholders '[TREATMENT]', '[REPLICATE]' and '[TIME]' substituted.

Author(s)

William Gustafson

See Also

re_esc

Examples

column.names<-c(
	'CTL.ZT1.1',
	'CTL.ZT1.2',
	'FRZ.ZT1.1',
	'FRZ.ZT1.2',
	'CTL.ZT4.1',
	'CTL.ZT4.2',
	'FRZ.ZT4.1',
	'FRZ.ZT4.2'
)
col.format<-'[TREATMENT].ZT[TIME].[REPLICATE]'
#grab time 1 samples
time1.regex<-re_esc(col_format_sub(col.format=col.format, ti='1'));
column.names[grep(time1.regex,column.names)]

#remove replicate
time.trt.regex<-re_esc( col_format_sub(col.format=col.format,ti='(.*)',tr='(.*)') )
column.names<-sub(time.trt.regex,'\1.ZT\2',column.names)

Run a single instance of the DiPALM pipeline with replicates permuted.

Description

Applys the DiPALM pipeline to a given set of gene expression data with replicates permuted within each timepoint. The user may want to instead use the frontend function main which runs this function many times in parallel.

Usage

dipalm(
	expr.file,
	 wgcna.power,
	 threads,
	 col.format,
	 ctrl.treatment,
	 p.value,
	 .log = open_log()
)

Arguments

Value

A list with the following components:

Author(s)

William Gustafson

Examples

expr.file<-tempfile();
data('smallTestData')
write.csv(smallTestData,file=expr.file);

dipalm(
	expr.file,
	7,
	parallel::detectCores(),
	'[TREATMENT].[REPLICATE].[TIME]',
	'CTL',
	0.05,
	open_log()
);

Example Data: Data for use with the DiPALM vignette

Description

A list of data objects that are used to run through the DiPALM example vignette.

Usage

data("exampleData")

Format

A List of 3 objects:

1. $rawCounts - A data frame [42383 rows X 48 columns] containing raw counts for mRNA expression data on Brassica rapa R500.

   • rows : Each row represents the expression of one gene across many different samples.

   • columns: Each column represents a different sample. Both drought-treated and properly watered samples are present, with two replicates of each time-course. Tissue was sampled at 4-hour intervals and column names are encoded as follows: S[replicate number] [D-Drought / W-Watered] [Time - sample] Example: S2D10 is part of the second replcate of the drought data and is the 10th time-point (ZT13, 13th hour after dawn on day two).

2. $moduleEigengenes - A dataframe [12 rows and 90 columns] representing 90 different expression vectors that describe descrete expression patterns found in the whole dataset ($timeCourseList).

   • Rows: each row represents a single time-point.

   • Columns: each column represents a distinct expression pattern seen repeatedly in the full data.

3. $geneAnnotations - A data frame [42383 rows and 5 columns] giving the locations within the Brassica rapa R500 genome where each gene is found. This is in SAF format which is used by the Linux Subread software package.

   • rows: Each row represents a single gene.

   • $GeneID: (character string) Unique gene accession.

   • $Chr: (character string) Chromosome number.

   • $Start: (integer) The position where the gene starts (first basepair).

   • $End: (integer) The position where the gene ends (last basepair).

   • $Strand: (character) is the gene on the positive [+] or negative [-] strand.

Details

Raw count data was summarized by mapping RNA-seq counts to the R500 Brassica rapa genome using the hisat2 aligner. Mapped reads were counted using the featureCounts function from the Linux Subread software package.

Source

Data is from: Greenham et al., eLife 2017;6:e29655

Plot multiple density distributions

Description

This function generates one or more smoothed density distributions using ggplot2 functionality

Usage

ggPlotMultiDensities(denslist, main = "", xlab = "", ylab = "Normalized Frequency",
 scale = T, cols = c("red", "blue", "grey50", "black", "skyblue", "orange"),
 ledgx = "topright", ltype = NULL, lwidth = NULL, xrng = NULL, yrng = NULL, Ledge = T)

Arguments

Details

This function takes in a list of one or more numeric vectors. Each vector is considered a separate population and a smoothed density plot (similar to a histogram) will be generated for each and plotted together. The actual numeric y-axis values are meaningless. A density plot is normalized so that the total area under the curve is equal to 1. When scale = T, the total area under the density plots may not equal one, but instead the total area is scaled to the relative proportions of the total population sizes of each population.

Value

Generates multiple density plots on the same graph.

Author(s)

Ryan C. Sartor

See Also

ggplot2

Examples

data(testData)
require(ggplot2)

# Two populations of the same size, unscaled
ggPlotMultiDensities(denslist = list(Test=testData$testResults,Permuted=testData$permutedResults),
main = "Pattern Change Scores", xlab = "Differential Pattern Score",lwidth = 1)

# Two populations of different sizes, scaled
ggPlotMultiDensities(denslist = list(Test=testData$testResults[1:3000],
Permuted=testData$permutedResults), scale=TRUE , main = "Pattern Change Scores (Scaled)",
xlab = "Differential Pattern Score",lwidth = 1)

Joins outputs from multiple invocations of the function dipalm.

Description

Given a collection of RData files containing the output of multiple invocations of the function dipalm this function combines them into a single set of results.

Usage

joinOutputs(files, p, n, .log = open_log())

Arguments

Value

A list with the following components.

Author(s)

William Gustafson

Examples

expr.file<-tempfile();
dipalm.file.1<-tempfile();
dipalm.file.2<-tempfile();
data('smallTestData')
write.csv(smallTestData,file=expr.file);
tmp<-dipalm(
	expr.file,
	7,
	parallel::detectCores(),
	'[TREATMENT].[REPLICATE].[TIME]',
	'CTL',
	0.05,
	open_log()
);
#save each component of the list returned by the function dipalm
do.call(save,as.list(c(names(tmp),envir=as.environment(tmp),file=dipalm.file.1)))
tmp<-dipalm(
	expr.file,
	7,
	parallel::detectCores(),
	'[TREATMENT].[REPLICATE].[TIME]',
	'CTL',
	0.05,
	open_log()
);
#save each component of the list returned by the function dipalm
do.call(save,as.list(c(names(tmp),envir=as.environment(tmp),file=dipalm.file.2)))
results<-joinOutputs(c(dipalm.file.1,dipalm.file.2),0.05,2);

Runs the full DiPALM pipeline with replicate permutation in parallel.

Description

The DiPALM pipeline uses timecourse gene expression data to test whether genes exhibit a differential pattern of expression between two treatments. The method is sensitive to the labeling of replicates which is typically arbitrary and thus unrelated between different timepoints. This function runs the pipeline a specified number of times in parallel and on each individual invocation of the pipeline permutes replicates separately within each timepoint.

Usage

main(
	expr.file,
	threads,
	wgcna.power,
	iterations,
	p.value,
	sig.ratio,
	col.format = "[TREATMENT].[REPLICATE].[TIME]",
	ctrl.treatment = c(),
	file.patt = "%d.dipalm.RData",
	tmpdir = ".",
	logdir = "logs",
	clear.old.files = FALSE
)

Arguments

Value

Returns the result of calling joinOutputs on the files containing the iteration outputs. This is a list with the following values:

Author(s)

William Gustafson

Examples

data("smallTestData");
expr.file<-tempfile();
dipalm.dir<-tempdir();
logdir<-tempdir();
write.csv(smallTestData,file=expr.file);

main(
	expr.file,
	2,
	7,
	2,
	0.05,
	0.9,
	'[TREATMENT].[REPLICATE].[TIME]',
	'CTL',
	tmpdir=dipalm.dir,
	logdir=logdir
);

Break a data matrix into time courses separated by treatment and replicate.

Description

Given a data set split into treatment groups, replicates and timepoints this function separates the dataset by treatment groups and replicate labels.

Usage

make_TCs(
	gene.expr,
	time_order = NULL,
	col.format = "[TREATMENT].[REPLICATE].[TIME]",
	.log = open_log()
)

Arguments

Value

A list indexed by all pairs of a treatment and a replicate label from the input dataset. Each component of the returned list is a matrix containing the samples from the corresponding treatment and replicate label pair. The column names of these matrices are the timepoints from the input dataset.

Author(s)

William Gustafson

See Also

open_log

Examples

gene.expr<-matrix(1:16,nrow=2)
colnames(gene.expr)<-c(
	'CTL.1.0',
	'FRZ.1.0',
	'CTL.1.4',
	'FRZ.1.4',
	'CTL.2.0',
	'FRZ.2.0',
	'CTL.2.4',
	'FRZ.2.4'
)
make_TCs(gene.expr,col.format='[TREATMENT].[REPLICATE].[TIME]',.log=open_log())

Regex substitution applied to a vector or list.

Description

A convenient frontend for the function sub that applies the function to each element of a vector or list and returns the result in the same form.

Usage

msub(regex, repl, ...)

Arguments

Value

A vector of strings resulting from applying the regular expression substutition to each component of the given vector.

Author(s)

William Gustafson

See Also

sub

Examples

msub("^CTL\\.(.*)\\.(.*)$","\\1.\\2",c("CTL.ZT1.1","FRZ.ZT1.1","CTL.ZT4.2","FRZ.ZT4.2"));

Open a log

Description

Returns a function, which is a wrapper around the function cat, to be used for logging.

Usage

open_log(log_file = NULL)

Arguments

Value

A function that can be used to write to the log file passed to open_log. This function is a wrapper for cat that sets the argument file to the value passed for the parameter log_file. The function also writes a newline provided the argument sep was not set to be a newline (by default sep is a single space). Calling the returned function with NULL as the first argument will close the log file.

If the value of log_file is NULL the returned function will print to standard output instead of writing to a file.

Author(s)

William Gustafson

See Also

cat

Examples

logfile<-tempfile()
.log<-open_log(logfile)
.log('This string is written to the file `log`.');
.log(NULL); #close the file `log`
.log<-open_log()
.log('This string is written to standard output.');

Format a string as a regular expression literal.

Description

Given a string returns escapes special characters used in regular expressions.

Usage

re_esc(s)

Arguments

Value

A string with special characters escaped.

Author(s)

William Gustafson

Examples

re_esc('x.csv')
re_esc('abc')
re_esc('a\\(.*)\\')

Samples replicates from a dataset without replacement.

Description

Given a dataset as a matrix or data.frame whose columns are samples grouped by treatments, times and replicates this function permutes the replicates which are typically arbitrary labels. Treatments and times of each sample are fixed.

Usage

sampleReplicates(gene.expr, col.format, .log)

Arguments

Value

A matrix containing the permuted dataset.

Author(s)

William Gustafson

See Also

col_format_sub open_log

Examples

set.seed(132)
gene.expr<-matrix(1:24,nrow=3)
colnames(gene.expr)<-c(paste0('CTL.ZT1.',1:4),paste0('CTL.ZT4.',1:4))
sampleReplicates(gene.expr,'[TREATMENT].ZT[TIME].[REPLICATE]',open_log());

Small test data for the DiPALM pipeline.

Description

A selection of 250 Arabidopsis Thaliana transcript profiles from a timecourse study on cold acclimation. When applying the DiPALM pipeline roughly 100 genes should be called significant at a p-value of 0.05. Furthermore, the genes called significant are expected to be concentrated within the last 125 genes.

Usage

data("smallTestData")

Format

A data frame with 250 observations on the 48 variables. Observations, which correspond to genes, are labeled by locus ids, e.g. AT5G64650. Variables, which correspond to samples, are labeled in the format [TREATMENT].[REPLICATE].[TIME]. There are two treatments, CTL and FRZ, four replicates 1,2,3,4 and 6 times 1,5,9,13,17,21.

Examples

file<-tempfile();
data("smallTestData")
write.csv(smallTestData, file=file)
results<-dipalm(file, 7, 2, "[TREATMENT].[REPLICATE].[TIME]","CTL",0.05);
sum(results$AdjkMEs<0.05);
## maybe str(x) ; plot(x) ...

Test Data: Data for function testing

Description

A list of objects that can be used to test individual DiPALM functions.

Usage

data("testData")

Format

A List of 6 objects:

1. $timeCourseList - A List of 4 matrices containing mRNA expression data on Brassica rapa R500.

   • $Drought.R1 (Plants exposed to drought, Replicate 1) A matrix of time-course gene expression data [18428 Genes X 12 Time-points]

   • $Drought.R2 (Plants exposed to drought, Replicate 2) A matrix of time-course gene expression data [18428 Genes X 12 Time-points]

   • $Watered.R1 (Properly watered plants, Replicate 1) A matrix of time-course gene expression data [18428 Genes X 12 Time-points]

   • $Watered.R2 (Properly watered plants, Replicate 2) A matrix of time-course gene expression data [18428 Genes X 12 Time-points]

2. $moduleEigengenes - A dataframe [12 rows and 90 columns] representing 90 different expression vectors that describe descrete expression patterns found in the whole dataset ($timeCourseList).

   • Rows: each row represents a single time-point

   • Columns: each column represents a distinct expression pattern seen repeatedly in the full data

3. $modelDesign - A matrix [4 rows and 2 columns] used as a design matrix for a simple linear model. This matrix provides a mapping between the input instances (timeCourseList) and the two fixed effects fitted in the mode (Drought and Watered). This matrix was created with model.matrix.

4. $modelContrast - A character string specifying the contrast to be evaluated using the fitted model for each gene. This example is used to compare watered to drought conditions.

5. $testResults - A named numeric vector with 18428 values from DiPALM test results on actual data.

   • names: Gene accessions

   • values: The summed absolute values of t-values generated from limma model contrasts related to each eigengene pattern. This can be thought of as a DiPALM score for differential patterning.

6. $permutedResults- A named numeric vector with 18428 values from DiPALM test results on permuted data. This vector acts as a null distribution in order to estimate significance of actual test results.

   • names: Gene accessions

   • values: The summed absolute values of t-values generated from limma model contrasts related to each eigengene pattern using permuted expression data. This can be thought of as a null distribution of DiPALM scores for differential patterning.

Source

Data is from: Greenham et al., eLife 2017;6:e29655