Solanum betaceum annotation data

index.html preprocessing.html

1 Annotation version: 20171102

Any annotation data will need to come from our trinotate information… An interesting caveat, kallisto uses a transcript database. tximport is the library used to read the counts/transcript into R, it is able to take a set oF mappings between gene/transcript in order to provide gene counts rather than transcript.

putative_annotations <- load_trinotate_annotations("reference/trinotate.csv")
annot_by_tx <- putative_annotations
rownames(annot_by_tx) <- make.names(putative_annotations[["transcript_id"]], unique=TRUE)

## I split the annotations into one set keyed by gene IDs and one by transcript IDs.
## This way it is possible to make an expressionset of transcripts or genes by
## removing/keeping the tx_gene_map argument.

## Why would you want to do this? you might ask...  Well, sadly the annotation database
## Does not always have filled-in data for the first transcript for a given gene ID,
## therefore the annotation information for those genes will be lost when we merge
## the count tables / annotation tables by gene ID.
## This is of course also true for transcript IDs, but for a vastly smaller number of them.
transcript_gene_map <- putative_annotations[, c("transcript_id", "gene_id")]
k10tx_raw <- create_expt(metadata="sample_sheets/all_samples.xlsx",
                         gene_info=annot_by_tx)

## Reading the sample metadata.

## The sample definitions comprises: 10, 7 rows, columns.

## Reading count tables.

## Finished reading count tables.

## Matched 234330 annotations and counts.

## Bringing together the count matrix and gene information.

putative_annotations <- fData(k10tx_raw)  ## A cheap way to get back the unique gene IDs.
dim(exprs(k10tx_raw))

## [1] 234330     10

## On further examination, we decided to remove both samples 1003 and 1008
## which Sandra pointed out are in fact from the same flow cytometry isolation.
k8tx_raw <- create_expt(metadata="sample_sheets/kept_samples.xlsx",
                        gene_info=annot_by_tx)

## Reading the sample metadata.

## The sample definitions comprises: 8, 7 rows, columns.

## Reading count tables.

## Finished reading count tables.

## Matched 234330 annotations and counts.

## Bringing together the count matrix and gene information.

dim(exprs(k8tx_raw))

## [1] 234330      8

## If we do not use the tx_gene_map information, we get 234,330 transcripts in the data set.
## Instead, we get 58,486 genes.

## Perform a count/gene cutoff here
threshold <- 4
method <- "cbcb"  ## other methods include 'genefilter' and 'simple'
k10tx_count <- sm(normalize_expt(k10tx_raw, filter=method, thresh=threshold))
dim(exprs(k10tx_count))

## [1] 42143    10

k8tx_count <- sm(normalize_expt(k8tx_raw, filter=method, thresh=threshold))
dim(exprs(k8tx_count))

## [1] 38571     8

## First get rid of the rRNA
sb_rrna <- putative_annotations[["rrna_subunit"]] != ""
rrna_droppers <- rownames(putative_annotations)[sb_rrna]
k10tx_norrna <- exclude_genes_expt(expt=k10tx_raw, ids=rrna_droppers)

## Before removal, there were 234330 entries.

## Now there are 234309 entries.

## Percent kept: 99.990, 99.988, 99.991, 99.993, 99.992, 99.989, 99.989, 99.987, 99.993, 99.995

## Percent removed: 0.010, 0.012, 0.009, 0.007, 0.008, 0.011, 0.011, 0.013, 0.007, 0.005

k8tx_norrna <- exclude_genes_expt(expt=k8tx_raw, ids=rrna_droppers)

## Before removal, there were 234330 entries.

## Now there are 234309 entries.

## Percent kept: 99.990, 99.988, 99.991, 99.992, 99.989, 99.989, 99.987, 99.995

## Percent removed: 0.010, 0.012, 0.009, 0.008, 0.011, 0.011, 0.013, 0.005

## Keep only the blastx hits with a low e-value and high identity.
filtered_keepers <- putative_annotations[["blastx_evalue"]] <= 1e-10 &
  putative_annotations[["blastx_identity"]] >= 60
keepers <- rownames(putative_annotations)[filtered_keepers]
## This excludes all but 33,451 transcripts.

## Now keep only those transcript IDs which fall into the above categories.
k10tx_conf <- exclude_genes_expt(expt=k10tx_norrna, ids=keepers, method="keep")

## Before removal, there were 234309 entries.

## Now there are 33451 entries.

## Percent kept: 25.523, 25.012, 26.033, 25.603, 26.420, 26.603, 27.573, 24.233, 24.821, 27.304

## Percent removed: 74.477, 74.988, 73.967, 74.397, 73.580, 73.397, 72.427, 75.767, 75.179, 72.696

dim(exprs(k10tx_conf))

## [1] 33451    10

k8tx_conf <- exclude_genes_expt(expt=k8tx_norrna, ids=keepers, method="keep")

## Before removal, there were 234309 entries.

## Now there are 33451 entries.

## Percent kept: 25.523, 25.012, 26.033, 26.420, 26.603, 27.573, 24.233, 27.304

## Percent removed: 74.477, 74.988, 73.967, 73.580, 73.397, 72.427, 75.767, 72.696

dim(exprs(k8tx_conf))

## [1] 33451     8

## Perform both the confidence and count filtration
k10tx_conf_count <- sm(normalize_expt(k10tx_conf, filter=method, thresh=threshold))
dim(exprs(k10tx_conf_count))

## [1] 18018    10

k8tx_conf_count <- sm(normalize_expt(k8tx_conf, filter=method, thresh=threshold))
dim(exprs(k8tx_conf_count))

## [1] 16826     8

rrnatx_expt <- exclude_genes_expt(expt=k10tx_raw, ids=sb_rrna, method="keep")

## Before removal, there were 234330 entries.

## Now there are 21 entries.

## Percent kept: 0.010, 0.012, 0.009, 0.007, 0.008, 0.011, 0.011, 0.013, 0.007, 0.005

## Percent removed: 99.990, 99.988, 99.991, 99.993, 99.992, 99.989, 99.989, 99.987, 99.993, 99.995

## This has only 21 putative rRNA transcripts.

2 Plot changes in number oF putative genes

genes_by_subset <- data.frame(
  samples = c(rep("ten", 4), rep("eight", 4)),
  colors = c(rep("#1B9E77", 4), rep("#7570B3", 4)),
  alpha = c("#1B9E77FF", "#1B9E77CC", "#1B0E77AA",
            "#1B0E7777", "#7570B3FF", "#7570B3CC", "#7570B3AA", "#7570B377"),
  type = rep(c("raw", "count", "conf", "conf_count"), 2),
  genes = c(nrow(exprs(k10tx_raw)), nrow(exprs(k10tx_count)),
            nrow(exprs(k10tx_conf)), nrow(exprs(k10tx_conf_count)),
            nrow(exprs(k8tx_raw)), nrow(exprs(k8tx_count)),
            nrow(exprs(k8tx_conf)), nrow(exprs(k8tx_conf_count))))

library(ggplot2)
columns <- ggplot(genes_by_subset, aes(x=samples, y=genes)) +
  geom_col(position="identity", aes(fill=colors, alpha=alpha), color="black") +
  scale_fill_manual(values=c(levels(as.factor(genes_by_subset$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

c1 <- genes_by_subset[["samples"]] == "ten"
c2 <- genes_by_subset[["samples"]] == "eight"
tmp <- cbind(genes_by_subset[c1, "genes"], genes_by_subset[c2, "genes"])
colnames(tmp) <- c("ten", "eight")
rownames(tmp) <- c("raw", "count", "conf", "conf_count")
tmp

##               ten  eight
## raw        234330 234330
## count       42143  38571
## conf        33451  33451
## conf_count  18018  16826

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: ggplot2(v.2.2.1) and hpgltools(v.2018.03)

loaded via a namespace (and not attached): Rcpp(v.0.12.15), RColorBrewer(v.1.1-2), bindr(v.0.1.1), pillar(v.1.2.1), compiler(v.3.4.4), plyr(v.1.8.4), zlibbioc(v.1.24.0), base64enc(v.0.1-3), iterators(v.1.0.9), tools(v.3.4.4), digest(v.0.6.15), rhdf5(v.2.22.0), evaluate(v.0.10.1), memoise(v.1.1.0), tibble(v.1.4.2), gtable(v.0.2.0), pkgconfig(v.2.0.1), rlang(v.0.2.0), openxlsx(v.4.0.17), foreach(v.1.4.4), commonmark(v.1.4), yaml(v.2.1.18), parallel(v.3.4.4), bindrcpp(v.0.2), dplyr(v.0.7.4), withr(v.2.1.1), stringr(v.1.3.0), knitr(v.1.20), roxygen2(v.6.0.1), xml2(v.1.2.0), hms(v.0.4.2), devtools(v.1.13.5), tidyselect(v.0.2.4), rprojroot(v.1.3-2), grid(v.3.4.4), glue(v.1.2.0), data.table(v.1.10.4-3), Biobase(v.2.38.0), R6(v.2.2.2), rmarkdown(v.1.9), pander(v.0.6.1), readr(v.1.1.1), purrr(v.0.2.4), tidyr(v.0.8.0), magrittr(v.1.5), scales(v.0.5.0), backports(v.1.1.2), codetools(v.0.2-15), htmltools(v.0.3.6), BiocGenerics(v.0.24.0), tximport(v.1.6.0), assertthat(v.0.2.0), colorspace(v.1.3-2), labeling(v.0.3), stringi(v.1.1.7), lazyeval(v.0.2.1) and munsell(v.0.4.3)

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 01_annotation_tx-v20171102.rda.xz

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