S. cerevisiae sample estimation, merged sample edition.

1 S. cerevisiae sample estimation of all samples (old and new)

This document should make clear the suitability of our yeast data for differential expression analyses. It should also give some ideas about the depth and distribution of the data.

2 Gathering samples

merged_expt <- sm(create_expt(metadata="sample_sheets/all_samples.xlsx",
                              gene_info=sc_all_annotations,
                              file_column="allfile"))
## Don't forget, I need to change the condition names.
merged_nor <- subset_expt(merged_expt, subset="batch!='r'")
merged_nos <- subset_expt(merged_expt, subset="batch!='s'")

3 Visualizing raw data

There are lots of toys we have learned to use to play with with raw data and explore stuff like batch effects or non-canonical distributions or skewed counts. hpgltools provides some functionality to make this process easier. The graphs shown below and many more are generated with the wrapper ‘graph_metrics()’ but that takes away the chance to explain the graphs as I generate them.

merged_filt <- sm(normalize_expt(merged_expt, filter=TRUE))

merged_libsize <- plot_libsize(merged_expt)

##tt <- pp(file="illustrator_input/libsize.pdf", image=merged_libsize$plot)
##dev.off()
merged_libsize$plot

merged_nonzero <- plot_nonzero(merged_filt)
tt <- pp(file="illustrator_input/nonzero.pdf", image=merged_nonzero$plot)

## Wrote the image to: illustrator_input/nonzero.pdf

merged_density <- plot_density(merged_filt)

## This data will benefit from being displayed on the log scale.

## If this is not desired, set scale='raw'

## Some entries are 0.  We are on log scale, setting them to 0.5.

## Changed 1152 zero count features.

tt <- pp(file="illustrator_input/raw_density.pdf", image=merged_density$plot)

## Wrote the image to: illustrator_input/raw_density.pdf

merged_boxplot <- plot_boxplot(merged_filt)

## This data will benefit from being displayed on the log scale.

## If this is not desired, set scale='raw'

## Some entries are 0.  We are on log scale, adding 1 to the data.

## Changed 1152 zero count features.

tt <- pp(file="illustrator_input/raw_boxplot.pdf", image=merged_boxplot)

## Wrote the image to: illustrator_input/raw_boxplot.pdf

topn <- plot_topn(merged_filt)
tt <- pp(file="illustrator_input/topn.pdf", image=topn$plot)

## Wrote the image to: illustrator_input/topn.pdf

prepost <- plot_libsize_prepost(merged_expt)
tt <- pp(file="illustrator_input/libsize_changes_lowfilter.pdf", image=prepost$lowgene_plot)

## Wrote the image to: illustrator_input/libsize_changes_lowfilter.pdf

4 Normalize and visualize

merged_norm <- normalize_expt(merged_expt, transform="log2", convert="cpm", filter=TRUE)

## This function will replace the expt$expressionset slot with:

## log2(cpm(hpgl(data)))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: hpgl

## Removing 1030 low-count genes (6095 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 1152 values equal to 0, adding 1 to the matrix.

## Step 5: not doing batch correction.

merged_metrics <- graph_metrics(merged_norm)

## Closing the pdf, X11cairo, cairo_pdf, cairo_pdf, cairo_pdf, cairo_pdf, png, cairo_pdf, cairo_pdf, cairo_pdf, cairo_pdf, cairo_pdf plotting device(s) before printing plots.

## Graphing number of non-zero genes with respect to CPM by library.

## Graphing library sizes.

## Graphing a boxplot.

## Graphing a correlation heatmap.

## Graphing a standard median correlation.

## Performing correlation.

## Graphing a distance heatmap.

## Graphing a standard median distance.

## Performing distance.

## Graphing a PCA plot.

## Graphing a T-SNE plot.

## Plotting a density plot.

## Plotting the representation of the top-n genes.

## Printing a color to condition legend.

The data should now be normalized, lets view some metrics post-facto.

I think something is very wrong, this should be printing pictures but seems to not be. Hmm I ran it manually and everything looks fine. I wonder what is going on? Lets re-render this and see what happens? I wonder if I left a graphics device open and didn’t realize it?

tt <- pp(file="illustrator_input/norm_corheat.pdf", image=merged_metrics$corheat)

## Wrote the image to: illustrator_input/norm_corheat.pdf

dev.off()

## pdf 
##   2

merged_metrics$corheat

tt <- pp(file="illustrator_input/norm_disheat.pdf", image=merged_metrics$disheat)

## Wrote the image to: illustrator_input/norm_disheat.pdf

## It appears that just the normalization is sufficient to split the samples completely by type and deeply separate them from the heterologous samples

tt <- pp(file="illustrator_input/norm_smc.pdf", image=merged_metrics$smc)

## Wrote the image to: illustrator_input/norm_smc.pdf

tt <- pp(file="illustrator_input/norm_smd.pdf", image=merged_metrics$smd)

## Wrote the image to: illustrator_input/norm_smd.pdf

## The samples are very well behaved, none fall below the red line.

tt <- pp(file="illustrator_input/norm_pca_test.pdf", image=merged_metrics$pcaplot)

## Wrote the image to: illustrator_input/norm_pca_test.pdf

tt <- pp(file="illustrator_input/norm_tsne.pdf", image=merged_metrics$tsneplot)

## Wrote the image to: illustrator_input/norm_tsne.pdf

## The homogeneous wt/mutant are nicely separated, and what is more, the exogeneous samples also split wt/mutant, that might prove to be quite useful.

tt <- pp(file="illustrator_input/legend.pdf", image=merged_metrics$legend)

## Wrote the image to: illustrator_input/legend.pdf

4.1 Print some variance partition information

This is an interesting aside which came up last week for some other data. It might be good to include the % variance correlated with ‘condition’ as a column in the annotation data for the expressionset. Thus, when we do the differential expression analysis later, we can look and see if a ‘significant’ gene has variance which is actually correlated with condition. I started to implement this idea, but it isn’t finished, perhaps it should be included in this?

vp <- varpart(merged_expt, predictor=NULL, factors=c("condition", "batch"))
tt <- pp(file="images/variance_partitions.png", image=vp$partition_plot)
merged_sva <- sm(normalize_expt(merged_expt, transform="log2",
                                filter=TRUE, batch="sva", low_to_zero=TRUE))
vpsva <- varpart(merged_sva, predictor=NULL, factors=c("condition", "batch"))
tt <- pp(file="images/variance_partitions_sva.png", image=vpsva$partition_plot)

5 Now look at the data without batch ‘r’

nor_norm <- sm(normalize_expt(merged_nor, transform="log2", convert="cpm", norm="quant", filter=TRUE))
nor_plots <- sm(graph_metrics(nor_norm))
nor_batchnorm <- sm(normalize_expt(merged_nor, transform="log2", convert="cpm", norm="quant", filter=TRUE,
                                batch="fsva"))
nor_batchplots <- sm(graph_metrics(nor_batchnorm))

## Before fsva
tt <- pp(file="illustrator_input/nor_corheat.pdf", image=nor_plots$corheat)

## Wrote the image to: illustrator_input/nor_corheat.pdf

tt <- pp(file="illustrator_input/nor_pca.pdf", image=nor_plots$pcaplot)

## Wrote the image to: illustrator_input/nor_pca.pdf

## After fsva
tt <- pp(file="illustrator_input/norbatch_corheat.pdf", image=nor_batchplots$corheat)

## Wrote the image to: illustrator_input/norbatch_corheat.pdf

tt <- pp(file="illustrator_input/norbatch_pca.pdf", image=nor_batchplots$pcaplot)

## Wrote the image to: illustrator_input/norbatch_pca.pdf

6 Try many batch estimation methods to see what is up.

I do not think any of this information is currently being used, so I am going to stop running it for the moment but keep it here in case we decide to revive it.

merged_pcabatch <- sm(normalize_expt(merged_expt, transform="log2",
                                filter=TRUE, batch="pca", low_to_zero=TRUE))
merged_pcabatch_metrics <- sm(graph_metrics(merged_pcabatch))

merged_nor_pcabatch <- sm(normalize_expt(merged_nor, transform="log2",
                                filter=TRUE, batch="pca", low_to_zero=TRUE))
merged_nor_pcabatch_metrics <- sm(graph_metrics(merged_nor_pcabatch))

merged_nos_pcabatch <- sm(normalize_expt(merged_nos, transform="log2",
                                filter=TRUE, batch="pca", low_to_zero=TRUE))
merged_nos_pcabatch_metrics <- sm(graph_metrics(merged_nos_pcabatch))

merged_sva <- sm(normalize_expt(merged_expt, transform="log2",
                                filter=TRUE, batch="sva", low_to_zero=TRUE))
merged_sva_metrics <- sm(graph_metrics(merged_sva))

merged_nor_sva <- sm(normalize_expt(merged_nor, transform="log2",
                                filter=TRUE, batch="sva", low_to_zero=TRUE))
merged_nor_sva_metrics <- sm(graph_metrics(merged_nor_sva))

merged_nos_sva <- sm(normalize_expt(merged_nos, transform="log2",
                                filter=TRUE, batch="sva", low_to_zero=TRUE))
merged_nos_sva_metrics <- sm(graph_metrics(merged_nos_sva))

merged_combat <- sm(normalize_expt(merged_expt, transform="log2",
                                   filter=TRUE, batch="combat_scale", low_to_zero=TRUE))
merged_combat_metrics <- sm(graph_metrics(merged_combat))

merged_nor_combat <- sm(normalize_expt(merged_nor, transform="log2",
                                   filter=TRUE, batch="combat_scale", low_to_zero=TRUE))
merged_nor_combat_metrics <- sm(graph_metrics(merged_nor_combat))

merged_nos_combat <- sm(normalize_expt(merged_nos, transform="log2",
                                   filter=TRUE, batch="combat_scale", low_to_zero=TRUE))
merged_nos_combat_metrics <- sm(graph_metrics(merged_nos_combat))

merged_limma <- sm(normalize_expt(merged_expt, transform="log2",
                                   filter=TRUE, batch="limmaresid", low_to_zero=TRUE))
merged_limma_metrics <- sm(graph_metrics(merged_limma))

merged_nor_limma <- sm(normalize_expt(merged_nor, transform="log2",
                                   filter=TRUE, batch="limmaresid", low_to_zero=TRUE))
merged_nor_limma_metrics <- sm(graph_metrics(merged_nor_limma))

merged_nos_limma <- sm(normalize_expt(merged_nos, transform="log2",
                                   filter=TRUE, batch="limmaresid", low_to_zero=TRUE))
merged_nos_limma_metrics <- sm(graph_metrics(merged_nos_limma))

merged_fsva <- sm(normalize_expt(merged_expt, transform="log2",
                              filter=TRUE, batch="fsva", low_to_zero=TRUE))
merged_fsva_metrics <- sm(graph_metrics(merged_fsva))

merged_nor_fsva <- sm(normalize_expt(merged_nor, transform="log2",
                              filter=TRUE, batch="fsva", low_to_zero=TRUE))
merged_nor_fsva_metrics <- sm(graph_metrics(merged_nor_fsva))

merged_nos_fsva <- sm(normalize_expt(merged_nos, transform="log2",
                              filter=TRUE, batch="ssva", low_to_zero=TRUE))
merged_nos_fsva_metrics <- sm(graph_metrics(merged_nos_fsva))

6.0.1 The resulting plots

Why is is suddenly printing stuff now and not before?

merged_pcabatch_metrics$pcaplot

merged_nor_pcabatch_metrics$pcaplot

merged_nos_pcabatch_metrics$pcaplot

merged_sva_metrics$pcaplot

merged_nor_sva_metrics$pcaplot

merged_nos_sva_metrics$pcaplot

merged_combat_metrics$pcaplot

merged_nor_combat_metrics$pcaplot

merged_nos_combat_metrics$pcaplot

merged_limma_metrics$pcaplot

merged_nor_limma_metrics$pcaplot

merged_nos_limma_metrics$pcaplot

merged_fsva_metrics$pcaplot

merged_nor_fsva_metrics$pcaplot

merged_nos_fsva_metrics$pcaplot

6.1 Look at specific normalizations via pca/tsne

merged_nor_fsva_metrics$pcaplot

7 Write the expt

## hmm I was going to silence this line, but looking at the report it seems to be missing some pictures.
fun_nor <- write_expt(merged_nor, excel=paste0("excel/samples_written_merged_nor-v", ver, ".xlsx"),
                      filter=TRUE, norm="raw", convert="cpm", batch="fsva", transform="log2")