COPDSexualDimorphism (version 1.3.0)

sdcd: Sexual dimorphic and COPD differential analysis

Description

Given linear models (limma) from a stratified analysis, sdcd compares the coefficients of the main effects across sexes. It then reports on the markers with significant differences in the coefficients. Typically used in conjunction with limma.

Usage

sdcd(male.fit, female.fit, coeff, genes, fdr.cutoff = 0.25, stat = c("z", "t"), file.prefix = "male.female.copd", class.names = c("male", "female"), write.file = TRUE)
sdcd.vmr(male.fit, female.fit, coeff, genes, fdr.cutoff = 0.25, stat = c("z", "t"), annotate = FALSE, annotate.with = c("genes","NCBI"), file.prefix = "male.female.copd", class.names = c("male", "female"), write.file = TRUE)
sdcd.core(male.fit, female.fit, coeff, stat = c("z", "t"))

Arguments

male.fit, female.fit

Objects of type limma, as generated by the function eBayes. OR a list with fields: coefficients (matrix), stdev.unscaled (matrix), sigma (numeric vector), df.residual (numeric), and df.prior (numeric).

coeff

Coefficients of the main effect of interest. sex for COPD-strratified analysis and COPD for sex-stratified analysis. This should correspond to a column name of the matrices coefficients and stdev.unscaled in male.fit and female.fit.

genes

Annotation of the gene expression probes.

fdr.cutoff

Numeric cutoff for FDR q-values

stat

Choices between "z" and "t". For "z", a z-test is used to assess significance of the difference between the regression coefficients. For "t", the t-statistics, as opposed to the coefficients themselves, are contrasted.

annotate

For SDCD methylation analyis, a boolean to determin if the methylated regions (VMRs) should be annotated by genes within 10kb or not.

annotate.with

For SDCD methylation analyis, when annotate == TRUE this option indicates how to annotate the methylated regions. The "genes" option uses the input object while the "NCBI" option uses function GetNeighGenes in the package NCBI2R.

file.prefix

Prefix for output file name.

class.names

An array of character strings of length two representing the two strata.

write.file

A boolean that determined weather the results should be written out as files.

Value

A data.frame with gene information and the following columns:

*.beta

Coefficients of the main effects in the two strata

*.sd

Standard deviation of the coefficients of the main effects in the two strata

*.t

T-statistics associated with the coefficients of the main effects in the two strata

*.p.value

P-value associated with the coefficients of the main effects in the two strata

beta.diff

The difference between the coefficients in the two strata

beta.diff.pooled.sd

Standard deviation of beta.diff

stratum1.stratum2.p

P-value associated with beta.diff

stratum1.stratum2.p.adj

Benjamini-Hochberg FDR corrected p-value

Details

This is to be used in conjuction with the package limma. Linear model fits are passed to the sdcd function. The main function is for gene expression data while the function sdcd.vmr is for methylation data. The main functionality is the same between the two, but the data annotation and output are slightly different. See intended usage in the vignette.

References

Sathirapongsasuti JF, Glass K, Huttenhower C, Quackenbush J, DeMeo DL. Integrative Genomics of Sexual Dimorphism in COPD. (In Prep).

See Also

lgrc.expr

Examples

data(lgrc.expr.meta)
data(lgrc.expr)
data(lgrc.genes)

## Sex-stratified
design.mtx = cbind(ctrl=1,
		copd=as.integer(grepl("COPD",colnames(expr))),
		age=expr.meta$age,
		pkyr=expr.meta$pkyrs)

good.idx = apply(design.mtx,1,function(x){!any(is.na(x))}) & (expr.meta$gender == "1-Male")
male.fit = lmFit(log(expr)[,good.idx], design.mtx[good.idx,])
male.fit = eBayes(male.fit)

good.idx = apply(design.mtx,1,function(x){!any(is.na(x))}) & (expr.meta$gender == "2-Female")
female.fit = lmFit(log(expr)[,good.idx], design.mtx[good.idx,])
female.fit = eBayes(female.fit)

male.female.copd.beta.diff.genes = sdcd(male.fit, female.fit, "copd", lgrc.genes, fdr.cutoff=0.25, file.prefix="male.female.copd")

## COPD-stratified
design.mtx = cbind(ctrl=1,
		gender=expr.meta$gender,
		age=expr.meta$age,
		pkyr=expr.meta$pkyrs)

good.idx = apply(design.mtx,1,function(x){!any(is.na(x))}) & grepl("COPD",colnames(expr))
copd.fit = lmFit(log(expr)[,good.idx], design.mtx[good.idx,])
copd.fit = eBayes(copd.fit)

good.idx = apply(design.mtx,1,function(x){!any(is.na(x))}) & grepl("CTRL",colnames(expr))
ctrl.fit = lmFit(log(expr)[,good.idx], design.mtx[good.idx,])
ctrl.fit = eBayes(ctrl.fit)

copd.ctrl.gender.beta.diff.genes = sdcd(copd.fit, ctrl.fit, "gender", lgrc.genes, fdr.cutoff=0.25, file.prefix="copd.ctrl.gender", class.names=c("copd","ctrl"))

## Combine
sdcd.genes = merge(copd.ctrl.gender.beta.diff.genes, male.female.copd.beta.diff.genes, by=setdiff(intersect(names(copd.ctrl.gender.beta.diff.genes), names(male.female.copd.beta.diff.genes)),c("beta.diff","beta.diff.pooled.sd")))
sdcd.genes = unique(sdcd.genes)
print(paste("There are ", nrow(sdcd.genes), " SDCD genes"))