Solanum betaceum sample estimation.

index.html preprocessing.html 01_annotation.html

1 Sample Estimation version: 20171102

Since I used kallisto for this work, I want to figure out tximport.

In 01_annotation.Rmd, I decided to open up a few lines of inquiry:

• k10_raw: All samples using the full set of putative annotations.
• k8_raw: Samples excluding hpgl1003 with the full set of annotations.
• k10_count: All samples, removing those with < ‘threshold (5)’ counts/gene.
• k8_count: 8 samples, ibid.
• k10_conf: All samples, but removing transcripts with annotations less confident than 1e-10 blastx e-value and less identity than 60 percent to some other known gene.
• k8_conf: 8 samples, ibid.
• k10_count_conf: A combination of the count filter and confidence filter.
• k8_count_conf: ibid, but with the 8 samples.

2 Examine these data sets

2.1 Raw data

Start with the full data set. This will take a rather long time to perform all the graphs, and since we are not likely to use it, I will limit myself to just plotting 1-2 graphs.

2.1.1 k10

k10_libsize <- plot_libsize(k10_raw)
k10_libsize$plot

k10_density <- sm(plot_density(k10_raw))
sm(k10_density$plot + ggplot2::scale_x_continuous(limits=c(1,100)))

k10_boxplot <- sm(plot_boxplot(k10_raw))
k10_boxplot

k10_corheat <- plot_corheat(k10_raw)

k10_tsne <- plot_tsne(k10_raw)
k10_tsne$plot

2.1.2 k8

k8_libsize <- plot_libsize(k8_raw)
k8_libsize$plot

k8_density <- sm(plot_density(k8_raw))
sm(k8_density$plot + ggplot2::scale_x_continuous(limits=c(1,100)))

k8_boxplot <- sm(plot_boxplot(k8_raw))
k8_boxplot

k8_corheat <- plot_corheat(k8_raw)

k8_tsne <- plot_tsne(k8_raw)
k8_tsne$plot

2.2 Count filtered data

When we get to the count filtered data, I suspect things will be more interesting. Therefore let us look at the full spectrum of plots before/after normalization.

2.2.1 k10

k10_count_raw_plots <- sm(graph_metrics(k10_count))
k10_count_norm <- sm(normalize_expt(k10_count, transform="log2", norm="quant", convert="cpm"))
k10_count_norm_plots <- sm(graph_metrics(k10_count_norm))
k10_count_normbatch <- sm(normalize_expt(k10_count, batch="svaseq", transform="log2", convert="cpm"))
k10_count_normbatch_plots <- graph_metrics(k10_count_normbatch)

## Closing the 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.

2.2.1.1 Look at some plots

Let us see how the k10 count data appears

k10_count_raw_plots$libsize

sm(k10_count_raw_plots$density + ggplot2::scale_x_continuous(limits=c(1,200)))

k10_count_raw_plots$boxplot

k10_count_norm_plots$corheat
k10_count_norm_plots$smc

k10_count_norm_plots$disheat
k10_count_norm_plots$smd

k10_count_norm_plots$pcaplot

## See if the additional information in the 'batch' factor made a difference
k10_count_normbatch_plots$pcaplot

2.2.2 k8

k8_count_raw_plots <- sm(graph_metrics(k8_count))
k8_count_norm <- sm(normalize_expt(k8_count, transform="log2", norm="quant", convert="cpm"))
k8_count_norm_plots <- sm(graph_metrics(k8_count_norm))
k8_count_normbatch <- sm(normalize_expt(k8_count, batch="svaseq", transform="log2", convert="cpm"))
k8_count_normbatch_plots <- sm(graph_metrics(k8_count_normbatch))

2.2.2.1 Look at some plots

Let us see how the k8 count data appears

k8_count_raw_plots$libsize

sm(k8_count_raw_plots$density + ggplot2::scale_x_continuous(limits=c(1,200)))

k8_count_raw_plots$boxplot

k8_count_norm_plots$corheat
k8_count_norm_plots$smc

k8_count_norm_plots$disheat
k8_count_norm_plots$smd

k8_count_norm_plots$pcaplot

## See if the additional information in the 'batch' factor made a difference
k8_count_normbatch_plots$pcaplot

k8_count_normbatch_plots$tsneplot

2.3 Conf filtered data

When we get to the conf filtered data, I suspect things will be more interesting. Therefore let us look at the full spectrum of plots before/after normalization.

2.3.1 k10

k10_conf_raw_plots <- sm(graph_metrics(k10_conf))
k10_conf_norm <- sm(normalize_expt(k10_conf, transform="log2", norm="quant", convert="cpm"))
k10_conf_norm_plots <- sm(graph_metrics(k10_conf_norm))
k10_conf_normbatch <- sm(normalize_expt(k10_conf, batch="svaseq", transform="log2",
                                        convert="cpm", filter=TRUE))
k10_conf_normbatch_plots <- sm(graph_metrics(k10_conf_normbatch))

2.3.1.1 Look at some plots

Let us see how the k10 conf data appears

k10_conf_raw_plots$libsize

sm(k10_conf_raw_plots$density + ggplot2::scale_x_continuous(limits=c(1,100)))

k10_conf_raw_plots$boxplot

k10_conf_norm_plots$corheat
k10_conf_norm_plots$smc

k10_conf_norm_plots$disheat
k10_conf_norm_plots$smd

k10_conf_norm_plots$pcaplot

## See if the additional information in the 'batch' factor made a difference
k10_conf_normbatch_plots$pcaplot

2.3.2 k8

k8_conf_raw_plots <- sm(graph_metrics(k8_conf))
k8_conf_norm <- sm(normalize_expt(k8_conf, transform="log2", norm="quant", convert="cpm"))
k8_conf_norm_plots <- sm(graph_metrics(k8_conf_norm))
k8_conf_normbatch <- sm(normalize_expt(k8_conf, batch="svaseq", transform="log2",
                                       convert="cpm", filter=TRUE))
k8_conf_normbatch_plots <- sm(graph_metrics(k8_conf_normbatch))

2.3.2.1 Look at some plots

Let us see how the k8 conf data appears

k8_conf_raw_plots$libsize

sm(k8_conf_raw_plots$density + ggplot2::scale_x_continuous(limits=c(1,100)))

k8_conf_raw_plots$boxplot

k8_conf_norm_plots$corheat
k8_conf_norm_plots$smc

k8_conf_norm_plots$disheat
k8_conf_norm_plots$smd

k8_conf_norm_plots$pcaplot

k8_conf_norm_plots$tsneplot

## See if the additional information in the 'batch' factor made a difference
k8_conf_normbatch_plots$pcaplot

k8_conf_normbatch_plots$tsneplot

2.4 Count_Conf filtered data

When we get to the count_conf filtered data, I suspect things will be more interesting. Therefore let us look at the full spectrum of plots before/after normalization.

2.4.1 k10

k10_conf_count_raw_plots <- sm(graph_metrics(k10_conf_count))
k10_conf_count_norm <- sm(normalize_expt(k10_conf_count, transform="log2",
                                         norm="quant", convert="cpm"))
k10_conf_count_norm_plots <- sm(graph_metrics(k10_conf_count_norm))
k10_conf_count_normbatch <- sm(normalize_expt(k10_conf_count, batch="svaseq", transform="log2",
                                              convert="cpm", filter=TRUE))
k10_conf_count_normbatch_plots <- sm(graph_metrics(k10_conf_count_normbatch))

2.4.1.1 Look at some plots

Let us see how the k10 count_conf data appears

k10_conf_count_raw_plots$libsize

sm(k10_conf_count_raw_plots$density + ggplot2::scale_x_continuous(limits=c(1,100)))

k10_conf_count_raw_plots$boxplot

k10_conf_count_norm_plots$corheat
k10_conf_count_norm_plots$smc

k10_conf_count_norm_plots$disheat
k10_conf_count_norm_plots$smd

k10_conf_count_norm_plots$pcaplot

## See if the additional information in the 'batch' factor made a difference
k10_conf_count_normbatch_plots$pcaplot

k10_conf_count_normbatch_plots$tsneplot

2.4.2 k8

k8_conf_count_raw_plots <- sm(graph_metrics(k8_conf_count))
k8_conf_count_norm <- sm(normalize_expt(k8_conf_count, transform="log2",
                                        norm="quant", convert="cpm"))
k8_conf_count_norm_plots <- sm(graph_metrics(k8_conf_count_norm))
k8_conf_count_normbatch <- sm(normalize_expt(k8_conf_count, batch="svaseq",
                                             transform="log2", convert="cpm"))
k8_conf_count_normbatch_plots <- sm(graph_metrics(k8_conf_count_normbatch))

2.4.2.1 Look at some plots

Let us see how the k8 count_conf data appears

k8_conf_count_raw_plots$libsize

sm(k8_conf_count_raw_plots$density + ggplot2::scale_x_continuous(limits=c(1,100)))

k8_conf_count_raw_plots$boxplot

k8_conf_count_norm_plots$corheat
k8_conf_count_norm_plots$smc

k8_conf_count_norm_plots$disheat
k8_conf_count_norm_plots$smd

k8_conf_count_norm_plots$pcaplot

## See if the additional information in the 'batch' factor made a difference
k8_conf_count_normbatch_plots$pcaplot

k8_conf_count_normbatch_plots$tsneplot

3 Make a nice stacked libsize barplot

library(ggplot2)
maximum <- plot_libsize(k8_raw, text=FALSE)
maximum_alpha <- maximum$table
maximum_alpha$alpha <- alpha(maximum_alpha$colors, 0.4)
conf <- plot_libsize(k8_conf, text=FALSE)
conf_alpha <- conf$table
conf_alpha$alpha <- alpha(conf_alpha$colors, 0.6)
count <- plot_libsize(k8_count, text=FALSE)
count_alpha <- count$table
count_alpha$alpha <- alpha(count_alpha$colors, 0.8)
conf_count <- plot_libsize(k8_conf_count, text=FALSE)
conf_count_alpha <- conf_count$table
conf_count_alpha$alpha <- alpha(conf_count_alpha$colors, 1.0)

all <- rbind(maximum_alpha, conf_alpha,
             count_alpha, conf_count_alpha)

## The following will stack the library sizes from each subset one on top of the other.
## This is likely not desired, but it does provide sufficient room to print the library
## sizes if one wishes.
stacked <- ggplot(all, aes(x=id)) +
  geom_bar(aes(weight=sum, alpha=alpha, fill=colors), color="black") +
  scale_fill_manual(values=c(levels(as.factor(all$colors)))) +
  theme(axis.text=ggplot2::element_text(size=10, colour="black"),
        axis.text.x=ggplot2::element_text(angle=90, vjust=0.5),
        legend.position="none")
stacked

## In this instance, they are stacked on top of 0, which I suspect is the goal.
columns <- ggplot(all, aes(x=id, y=sum)) +
  geom_col(position="identity", aes(fill=colors, alpha=alpha), color="black") +
  scale_fill_manual(values=c(levels(as.factor(all$colors)))) +
  theme(axis.text=ggplot2::element_text(size=10, colour="black"),
        axis.text.x=ggplot2::element_text(angle=90, vjust=0.5),
        legend.position="none")
columns

4 Write up the favorites

Let us assume for a moment that we will choose a couple of the above for looking at more closely. I suspect k8_count_conf, k10_count_conf, and k8_count and k8_conf are good candidates.

output <- paste0("excel/k8_conf_count_by_tx-v", ver, ".xlsx")
if (!file.exists(output)) {
  k8_conf_count_written <- sm(write_expt(k8_conf_count, excel=output))
}
output <- paste0("excel/k10_conf_count_by_tx-v", ver, ".xlsx")
if (!file.exists(output)) {
  k10_conf_count_written <- sm(write_expt(k10_conf_count, excel=output))
}
output <- paste0("excel/k8_count_by_tx-v", ver, ".xlsx")
if (!file.exists(output)) {
  k8_count_written <- sm(write_expt(k8_count, excel=output))
}
output <- paste0("excel/k8_conf_by_tx-v", ver, ".xlsx")
if (!file.exists(output)) {
  k8_conf_written <- sm(write_expt(k8_conf, excel=output))
}

5 Quantify tx/gene

I want to see if there are some changes in the number of reads/tx between the embryogenic/nonembryogenic. I am not entirely certain how to do this, so I will poke about and see what happens.

I do have a data structure telling me the mapping of tx->gene.

transcript_gene_map <- putative_annotations[, c("transcript_id", "gene_id")]
mappings <- data.table::as.data.table(transcript_gene_map)
mappings[["num"]] <- 1
num_tx <- mappings[, list(num_tx=sum(num)), by="gene_id"]

tx_hist <- plot_histogram(num_tx$num_tx)
tx_hist

start <- normalize_expt(k8_raw, convert="cpm")

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

## cpm(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

## Filter is false, this should likely be set to something, good
##  choices include cbcb, kofa, pofa (anything but FALSE).  If you want this to
##  stay FALSE, keep in mind that if other normalizations are performed, then the
##  resulting libsizes are likely to be strange (potentially negative!)

## Leaving the data in its current base format, keep in mind that
##  some metrics are easier to see when the data is log2 transformed, but
##  EdgeR/DESeq do not accept transformed data.

## 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: not doing count filtering.

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: not transforming the data.

## Step 5: not doing batch correction.

embryo <- subset_expt(start, subset="condition=='embryogenic'")
nonembryo <- subset_expt(start, subset="condition=='nonembryogenic'")
embryo_df <- as.data.frame(exprs(embryo))
embryo_df$em_mean <- rowMeans(embryo_df)
non_df <- as.data.frame(exprs(nonembryo))
non_df$non_mean <- rowMeans(non_df)
e <- as.data.frame(embryo_df[["em_mean"]])
rownames(e) <- rownames(embryo_df)
colnames(e) <- "em_mean"
n <- as.data.frame(non_df[["non_mean"]])
rownames(n) <- rownames(non_df)
colnames(n) <- "non_mean"
means <- merge(x=e, y=n, by="row.names")
rownames(means) <- means[["Row.names"]]
means <- means[, -1]
near_zero <- means <= 0.001
embryo_near_zero <- sum(near_zero[, "em_mean"])
embryo_near_zero

## [1] 35080

non_near_zero <- sum(near_zero[, "non_mean"])
non_near_zero

## [1] 52196

plot_boxplot(embryo)

## 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 391753 zero count features.

plot_boxplot(nonembryo)

## 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 508520 zero count features.

6 Note to self

I decided to not save the state of this file, because it is huge. The differential expression analyses do not need this, but need only the result of the annotation worksheet.

pander::pander(sessionInfo())

R version 3.4.4 RC (2018-03-09 r74380)

**Platform:** x86_64-pc-linux-gnu (64-bit)

locale: LC_CTYPE=en_US.utf8, LC_NUMERIC=C, LC_TIME=en_US.utf8, LC_COLLATE=en_US.utf8, LC_MONETARY=en_US.utf8, LC_MESSAGES=en_US.utf8, LC_PAPER=en_US.utf8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.utf8 and LC_IDENTIFICATION=C

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: hpgltools(v.2017.10), ruv(v.0.9.6), data.table(v.1.10.4-3), beanplot(v.1.2), plyr(v.1.8.4) and ggplot2(v.2.2.1)

loaded via a namespace (and not attached): nlme(v.3.1-131.1), bitops(v.1.0-6), matrixStats(v.0.53.1), devtools(v.1.13.5), bit64(v.0.9-7), RColorBrewer(v.1.1-2), rprojroot(v.1.3-2), tools(v.3.4.4), backports(v.1.1.2), R6(v.2.2.2), KernSmooth(v.2.23-15), DBI(v.0.7), lazyeval(v.0.2.1), BiocGenerics(v.0.24.0), mgcv(v.1.8-23), colorspace(v.1.3-2), withr(v.2.1.1), tidyselect(v.0.2.4), bit(v.1.1-12), compiler(v.3.4.4), preprocessCore(v.1.40.0), Biobase(v.2.38.0), xml2(v.1.2.0), labeling(v.0.3), scales(v.0.5.0), readr(v.1.1.1), genefilter(v.1.60.0), quadprog(v.1.5-5), commonmark(v.1.4), stringr(v.1.3.0), digest(v.0.6.15), rmarkdown(v.1.8), base64enc(v.0.1-3), pkgconfig(v.2.0.1), htmltools(v.0.3.6), limma(v.3.34.9), rlang(v.0.2.0), RSQLite(v.2.0), bindr(v.0.1), BiocParallel(v.1.12.0), gtools(v.3.5.0), dplyr(v.0.7.4), RCurl(v.1.95-4.10), magrittr(v.1.5), Matrix(v.1.2-12), Rcpp(v.0.12.15), munsell(v.0.4.3), S4Vectors(v.0.16.0), stringi(v.1.1.6), yaml(v.2.1.16), edgeR(v.3.20.8), MASS(v.7.3-49), zlibbioc(v.1.24.0), rhdf5(v.2.22.0), Rtsne(v.0.13), grid(v.3.4.4), blob(v.1.1.0), parallel(v.3.4.4), ggrepel(v.0.7.0), lattice(v.0.20-35), splines(v.3.4.4), annotate(v.1.56.1), pander(v.0.6.1), hms(v.0.4.1), locfit(v.1.5-9.1), knitr(v.1.20), pillar(v.1.2.1), corpcor(v.1.6.9), reshape2(v.1.4.3), codetools(v.0.2-15), stats4(v.3.4.4), XML(v.3.98-1.10), glue(v.1.2.0), evaluate(v.0.10.1), foreach(v.1.4.4), gtable(v.0.2.0), purrr(v.0.2.4), tidyr(v.0.8.0), assertthat(v.0.2.0), openxlsx(v.4.0.17), xtable(v.1.8-2), roxygen2(v.6.0.1), survival(v.2.41-3), tibble(v.1.4.2), iterators(v.1.0.9), AnnotationDbi(v.1.40.0), memoise(v.1.1.0), IRanges(v.2.12.0), tximport(v.1.6.0), bindrcpp(v.0.2), sva(v.3.26.0) and directlabels(v.2017.03.31)

message(paste0("This is hpgltools commit: ", get_git_commit()))

## If you wish to reproduce this exact build of hpgltools, invoke the following:

## > git clone http://github.com/abelew/hpgltools.git

## > git reset e3a7f13d63d6906d825dbfad37273ddd71e44e86

## R> packrat::restore()

## This is hpgltools commit: Mon Mar 12 12:11:17 2018 -0400: e3a7f13d63d6906d825dbfad37273ddd71e44e86

this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to ", this_save))

## Saving to 02_sample_estimation_tx-v20171102.rda.xz

## tmp <- sm(saveme(filename=this_save))