Gallery: TCDD Sex Difference Dose Response

Correlation of gene fold-changes in male and female mice at 4 different TCDD dose groups was visualized using a scatter plot. The vast majority of gene fold-changes were highly correlation in male and female animal. But a subpopulation of data points that are off the diagonal represent genes that were significantly altered in one sex but not the other. Pearson correlation test was conducted between gene alteration in male and female group of each dose. In general, the correlation increased with higher TCDD dose treatment.

### scatterplot_script.R ##########################################################################
# Purpose: To generate scatter plot comparing gene alterations across male and female cohorts
 
### PREAMBLE ######################################################################################
# loading the affy library into R
library(BoutrosLab.plotting.general);
 
# unix-specific parameters
setwd('~/BoutrosLab/Toxicology/TCDD/mouse/Liver/mRNA/male_vs_female_dose-response/');
# general parameters
date <- Sys.Date();
cfg <- read.config.file();
 
# change '~/isilon/'
if (!file.exists(cfg$output.directory)) {
        print('gsub file names');
        cfg <- lapply( cfg, function(x) {
                gsub('~/isilon/', '/.mounts/labs/boutroslab/', x);
                });
        }
 
### READ IN THE FIT FILES #########################################################################
# set directory
setwd(cfg$output.directory);
 
# read in annotated fit results
fit.data <- read.table(
        '2015-01-08_Sex-TCDD_DoseResponse_fit_data.txt',
        header = TRUE,
        sep = "\t"
        );
 
# read in phenodata
phenodata.374 <- read.table(
        file = cfg$phenodata.file.374,
        sep = "\t",
        header = TRUE,
        row.names = 1
        );
 
# get phenodata
phenodata.386 <- read.table(
        file = cfg$phenodata.file.386,
        sep = "\t",
        header = TRUE,
        row.names = 1
        );
 
phenodata <- as.data.frame(rbind(phenodata.386, phenodata.374[,names(phenodata.386)]));
phenodata$Sex <- as.character('Male');
phenodata[grep('S386',rownames(phenodata)),]$Sex  <- 'Female';
 
### ORGANIZE DATA #################################################################################
# determine correlation between sexes at each dose
correlation.text <- c(
        as.expression(substitute(
                R['125'] == paste(correlation, phantom('|')[phantom('|')], phantom('|')^phantom('|')),
                env = list(correlation = round(cor(fit.data$Male.T125, fit.data$Female.T125, method = 'pearson'), 2))
                )),
        as.expression(substitute(
                R['250'] == paste(correlation, phantom('|')[phantom('|')], phantom('|')^phantom('|')),
                env = list(correlation = round(cor(fit.data$Male.T250, fit.data$Female.T250, method = 'pearson'), 2))
                )),
        as.expression(substitute(
                R['500'] == paste(correlation, phantom('|')[phantom('|')], phantom('|')^phantom('|')),
                env = list(correlation = round(cor(fit.data$Male.T500, fit.data$Female.T500, method = 'pearson'), 2))
                )),
        as.expression(substitute(
                R['1000'] == paste(correlation, phantom('|')[phantom('|')], phantom('|')^phantom('|')),
                env = list(correlation = round(cor(fit.data$Male.T1000, fit.data$Female.T1000, method = 'pearson'), 2))
                ))
        );
 
# maybe add p-values??
 
# use only fold-change values
plot.data <- data.frame(
        'Male.Coef'     = stack(fit.data[,c('Male.T125','Male.T250','Male.T500','Male.T1000')])$values,
        'Female.Coef'   = stack(fit.data[,c('Female.T125','Female.T250','Female.T500','Female.T1000')])$values,
        'Male.groups'   = stack(fit.data[,c('Male.T125','Male.T250','Male.T500','Male.T1000')])$ind,
        'Female.groups' = stack(fit.data[,c('Female.T125','Female.T250','Female.T500','Female.T1000')])$ind,
        'EntrezID'      = rep(fit.data$EntrezID,4),
        'Gene'          = rep(fit.data$Gene,4)
        );
 
# select 4 shades of green colors to represent doses
light.green <- rgb(red=(66/255),green=(255/255),blue=(66/255));
dark.green <- rgb(red=(0/255),green=(48/255),blue=(0/255));
ramp.green <- colorRamp(c(light.green,dark.green));
dose.colours <- rgb(ramp.green(seq(0,1,length=4)),max=255); 
 
# select 4 shades of red colors to highlight genes
light.red <- rgb(red=(255/255),green=(85/255),blue=(82/255));
dark.red <- rgb(red=(178/255),green=(0/255),blue=(0/255));
ramp.red <- colorRamp(c(light.red,dark.red));
highlight.colours <- rgb(ramp.red(seq(0,1,length=4)),max=255); 
 
# create correlation plot
create.scatterplot(
        Female.Coef ~ Male.Coef,
        plot.data,
        filename = generate.filename(cfg$project.stem, 'SexCorrelationPlot','tiff'),
        resolution = 1000,
        groups = plot.data$Male.groups,
        col = dose.colours,
        type = c('p','g'),
        pch = c(15,17,18,19),
        xlab.label = 'Male',
        ylab.label = 'Female',
        xlab.cex = 2.5,
        ylab.cex = 2.5,
        xaxis.cex = 2,
        yaxis.cex = 2,
        xlimits = c(-6,9),
        ylimits = c(-6,9),
        xat = seq(-6,9,3),
        yat = seq(-6,9,3),
        add.xyline = TRUE,
        xyline.col = 'grey',
        xyline.lwd = 1,
        xyline.lty = 'dashed',
        xgrid.at = 0,
        ygrid.at = 0,
        key = list(
                points = list(
                        pch = c(15,17,18,19),
                        col = dose.colours,
                        cex = 1.25
                        ),
                text = list(
                        lab = correlation.text,
                        cex = 1.25,
                        col = 'black'
                        ),
                x = 0.039,
                y = 0.99
                )
        );
 
### WRITE SESSION INFO TO FILE ####################################################################
save.session.profile(generate.filename(cfg$project.stem,'SexCorrelation_SessionInfo','txt'));

Created by Pretty R at inside-R.org