Some very old Angomonas deanei data.

1 Adeanei fun!

load("RData")
rm(list=ls())
save(list=ls(all=TRUE), file="RData")
render("adeanei.rmd", output_format="pdf_document")
render("adeanei.rmd", output_format="html_document")
render("adeanei.rmd", 'knitrBootstrap::bootstrap_document')

annotation <- load_gff_annotations("adeanei.gff")

## Trying attempt: rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo=TRUE)

## Trying attempt: rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo=FALSE)

## Had a successful gff import with rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo=FALSE)

## Returning a df with 16 columns and 16888 rows.

##annotation = as.data.frame(import.gff3("adeanei.gff"))
genes <- annotation[annotation$type=="gene",]
rownames(genes) <- genes$ID
tooltip_data <- genes
tooltip_data <- tooltip_data[,c("ID","product")]
tooltip_data$tooltip <- paste(tooltip_data$ID, tooltip_data$product, sep=": ")
tooltip_data$tooltip <- gsub("\\+", " ", tooltip_data$tooltip)
tooltip_data <- tooltip_data[-1]
tooltip_data <- tooltip_data[-1]
colnames(tooltip_data) <- c("name.tooltip")
head(tooltip_data)

##                                             name.tooltip
## AGDE_00001 AGDE_00001: peptide alpha-N-acetyltransferase
## AGDE_00002              AGDE_00002: hypothetical protein
## AGDE_00003              AGDE_00003: hypothetical protein
## AGDE_00004      AGDE_00004: exosome-associated protein 2
## AGDE_00005              AGDE_00005: hypothetical protein
## AGDE_00006              AGDE_00006: hypothetical protein

all <- create_expt("all_samples.csv", gene_info=genes)

## Reading the sample metadata.

## The sample definitions comprises: 7 rows(samples) and 6 columns(metadata fields).

## Reading count tables.

## Reading count files with read.table().

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/small_rna/adeanei/preprocessing/small_rnav2/hpgl0293/hpgl0293_adeanei.count.xz contains 16893 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/small_rna/adeanei/preprocessing/small_rnav2/hpgl0294/hpgl0294_adeanei.count.xz contains 16893 rows and merges to 16893 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/small_rna/adeanei/preprocessing/small_rnav2/hpgl0295/hpgl0295_adeanei.count.xz contains 16893 rows and merges to 16893 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/small_rna/adeanei/preprocessing/small_rnav2/hpgl0561/hpgl0561_adeanei.count.xz contains 16893 rows and merges to 16893 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/small_rna/adeanei/preprocessing/small_rnav2/hpgl0562/hpgl0562_adeanei.count.xz contains 16893 rows and merges to 16893 rows.

## rnaseq/hpgl0323/hpgl0323_forward-trimmed-v1M1l20.count.xz contains 16893 rows and merges to 16893 rows.

## rnaseq/hpgl0563/hpgl0563_forward-trimmed-v1M1l20.count.xz contains 16893 rows and merges to 16893 rows.

## Finished reading count data.

## Matched 16888 annotations and counts.

## Bringing together the count matrix and gene information.

## Some annotations were lost in merging, setting them to 'undefined'.

## Saving the expressionset to 'expt.rda'.

## The final expressionset has 16888 rows and 7 columns.

small <- subset_expt(expt=all, subset="type=='small'")

## Using a subset expression.

## There were 7, now there are 5 samples.

rnaseq <- subset_expt(expt=all, subset="type=='rnaseq'")

## Using a subset expression.

## There were 7, now there are 2 samples.

all_norm <- normalize_expt(all, transform="log2", norm="quant", convert="cpm", filter=TRUE)

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

## log2(cpm(quant(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is 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

## 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: cbcb

## Removing 175 low-count genes (16713 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: not doing batch correction.

plot_pca(all_norm)$plot

small_norm <- normalize_expt(expt=small, transform="log2", norm="quant", convert="cpm", filter=TRUE)

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

## log2(cpm(quant(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is 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

## 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: cbcb

## Removing 235 low-count genes (16653 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: not doing batch correction.

rnaseq_norm <- normalize_expt(expt=rnaseq, transform="log2", norm="quant", convert="cpm", filter=TRUE)

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

## log2(cpm(quant(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is 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

## 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: cbcb

## Removing 8719 low-count genes (8169 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## Step 5: not doing batch correction.

plot_pca(small_norm)$plot

small_pair <- all_pairwise(small)

## Plotting a PCA before surrogates/batch inclusion.

## Not putting labels on the plot.

## Using limma's removeBatchEffect to visualize with(out) batch inclusion.

## Not putting labels on the plot.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

small_tables <- combine_de_tables(small_pair, excel="small_pairwise.xlsx")

## Deleting the file small_pairwise.xlsx before writing the tables.

## Writing a legend of columns.

## Printing a pca plot before/after surrogates/batch estimation.

## Working on table 1/1: wt_vs_alpha

## Adding venn plots for wt_vs_alpha.

## Limma expression coefficients for wt_vs_alpha; R^2: 0.424; equation: y = 0.615x + 1.31

## Deseq expression coefficients for wt_vs_alpha; R^2: 0.408; equation: y = 0.491x + 1.73

## Edger expression coefficients for wt_vs_alpha; R^2: 0.39; equation: y = 0.442x + 2.64

## Writing summary information, compare_plot is: TRUE.

## Performing save of small_pairwise.xlsx.

small_tables[["plots"]][["wt_vs_alpha"]][["deseq_ma_plots"]][["plot"]]

## Compare mapping parameters on rpkm/cpm
one_sample = create_expt(file="rpkm_samples.csv")
lengths = annotation[,c("seqnames","width")]
one_rpkm = normalize_expt(expt=one_sample, convert="rpkm", annotations=genes)
one_cpm = normalize_expt(expt=one_sample, convert="cpm")
one_sf_rpkm = normalize_expt(expt=one_sample, convert="rpkm", norm="sf", annotations=genes)
one_sf_cpm = normalize_expt(expt=one_sample, convert="cpm", norm="sf")

df = exprs(one_sample$expressionset)
dfrpkm = hpgl_norm(df=df, convert="rpkm", annotations=genes)

df = exprs(one_sample$expressionset)
df = merge(df, exprs(one_cpm$expressionset), by="row.names")
rownames(df) = df$Row.names
df = merge(df, exprs(one_sf_rpkm$expressionset), by="row.names")
rownames(df) = df$Row.names
df = merge(df, exprs(one_sf_cpm$expressionset), by="row.names")
rownames(df) = df$Row.names
df = df[-1]
df = df[-1]
df = df[-1]

tiff(file="fun.tiff")
hpgl_boxplot(df=df, scale="log2")
dev.off()

pander::pander(sessionInfo())
message(paste0("This is hpgltools commit: ", get_git_commit()))
this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to ", this_save))
tmp <- sm(saveme(filename=this_save))