This worksheet aims to lay out the tasks I performed when analyzing some TNSeq data of a group B streptococcus.
It looks to me, that despite the oddities in processing the raw reads, there is nice coverage and some obviously essential genes. The next question: did any change status as more calprotectin was added?
The ID for strain a909 at microbesonline.org is: 205921.
Let us load up annotations from my gff file along with the microbesonline.
As a heads up, the count tables are using IDs which look like: SAK_RS00185. This appears to be the ‘sysName’ column from microbesonline and the locus_tag column from the gff annotations. In addition, there are a bunch of unused columns in both data sets which we likely want to prune.
Ahh, that is incorrect, the microbesonline ‘sysName’ is the same as ‘old_locus_tag’ column.
There are three relatively closely related strains which may be sufficiently similar to use in this analysis. The actual strain is cjb111, but that has not yet been quite finished, as far as I can tell. Therefore I will repeat most(all?) tasks with strains a909 and vr2603 to see if they may be more useful.
load_microbesonline_annotations(species="A909") a909_microbes <-
## Found 1 entry.
## Streptococcus agalactiae A909Firmicutesyes2006-01-18yes102136205921
## The species being downloaded is: Streptococcus agalactiae A909
## Downloading: http://www.microbesonline.org/cgi-bin/genomeInfo.cgi?tId=205921;export=tab
load_gff_annotations("reference/sagalactiae_a909_all.gff") a909_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 34 columns and 4426 rows.
as.data.frame(a909_microbes)
a909_microbes <-rownames(a909_gff) <- make.names(a909_gff[["locus_tag"]], unique=TRUE)
## I am going to only pay attention to the first annotation for each locus tag from microbesonline.
"sysName"]] <- make.names(a909_microbes[["sysName"]], unique=TRUE)
a909_microbes[[ merge(a909_gff, a909_microbes, by.x="old_locus_tag", by.y="sysName")
a909_annot <-rownames(a909_annot) <- make.names(a909_annot[["locus_tag"]], unique=TRUE)
## Rename the merged start/strand columns
colnames(a909_annot)[3] <- "start"
colnames(a909_annot)[6] <- "strand"
## And drop the duplicates
c(39, 41)] <- NULL a909_annot[,
load_microbesonline_annotations(species="CJB111") cjb111_microbes <-
## Found 1 entry.
## Streptococcus agalactiae CJB111Firmicutesyes2007-05-08noNANA2208342617
## The species being downloaded is: Streptococcus agalactiae CJB111
## Downloading: http://www.microbesonline.org/cgi-bin/genomeInfo.cgi?tId=342617;export=tab
load_gff_annotations("reference/sagalactiae_cjb111.gff") cjb111_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 14 columns and 2208 rows.
as.data.frame(cjb111_microbes)
cjb111_microbes <-rownames(cjb111_gff) <- make.names(cjb111_gff[["locus_tag"]], unique=TRUE)
## I am going to only pay attention to the first annotation for each locus tag from microbesonline.
"sysName"]] <- make.names(cjb111_microbes[["sysName"]], unique=TRUE)
cjb111_microbes[[ merge(cjb111_gff, cjb111_microbes, by.x="locus_tag", by.y="sysName")
cjb111_annot <-rownames(cjb111_annot) <- make.names(cjb111_annot[["locus_tag"]], unique=TRUE)
## Rename the merged start/strand columns
colnames(cjb111_annot)[3] <- "start"
colnames(cjb111_annot)[6] <- "strand"
c(19, 21)] <- NULL cjb111_annot[,
I think this might actually be 2603vr, I get confused, and under a few specific circumstances R acts strange when things start with numbers.
load_microbesonline_annotations(species="2603V") vr2603_microbes <-
## Found 1 entry.
## Streptococcus agalactiae 2603V/RFirmicutesyesyes102273208435
## The species being downloaded is: Streptococcus agalactiae 2603V/R
## Downloading: http://www.microbesonline.org/cgi-bin/genomeInfo.cgi?tId=208435;export=tab
load_gff_annotations("reference/sagalactiae_2603vr.gff") vr2603_gff <-
## Trying attempt: rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo=TRUE)
## Had a successful gff import with rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo=TRUE)
## Returning a df with 28 columns and 4611 rows.
as.data.frame(vr2603_microbes)
vr2603_microbes <-rownames(vr2603_gff) <- make.names(vr2603_gff[["locus_tag"]], unique=TRUE)
## I am going to only pay attention to the first annotation for each locus tag from microbesonline.
"sysName"]] <- make.names(vr2603_microbes[["sysName"]], unique=TRUE)
vr2603_microbes[[ merge(vr2603_gff, vr2603_microbes, by.x="locus_tag", by.y="sysName")
vr2603_annot <-rownames(vr2603_annot) <- make.names(vr2603_annot[["ID"]], unique=TRUE)
## Rename the merged start/strand columns
colnames(vr2603_annot)[3] <- "start"
colnames(vr2603_annot)[6] <- "strand"
c(33, 35)] <- NULL vr2603_annot[,
The following block merges the various counts, annotations, and experimental metadata.
Just as with the annotations, I will create one expressionset for each strain.
create_expt(metadata="sample_sheets/sagalactiae_samples.xlsx",
a909_expt <-batch=FALSE, gene_info=a909_annot, file_column="a909_filename")
## Reading the sample metadata.
## The sample definitions comprises: 9 rows(samples) and 30 columns(metadata fields).
## Reading count tables.
## Reading count files with read.table().
## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/tnseq/sagalacticae_2019/preprocessing/01/outputs/bowtie_sagalactiae_a909/trimmed_ca-v0M1.count.xz contains 2212 rows.
## preprocessing/02/outputs/bowtie_sagalactiae_a909/trimmed_ca-v0M1.count.xz contains 2212 rows and merges to 2212 rows.
## preprocessing/03/outputs/bowtie_sagalactiae_a909/trimmed_ca-v0M1.count.xz contains 2212 rows and merges to 2212 rows.
## preprocessing/04/outputs/bowtie_sagalactiae_a909/trimmed_ca-v0M1.count.xz contains 2212 rows and merges to 2212 rows.
## preprocessing/05/outputs/bowtie_sagalactiae_a909/trimmed_ca-v0M1.count.xz contains 2212 rows and merges to 2212 rows.
## preprocessing/06/outputs/bowtie_sagalactiae_a909/trimmed_ca-v0M1.count.xz contains 2212 rows and merges to 2212 rows.
## preprocessing/07/outputs/bowtie_sagalactiae_a909/trimmed_ca-v0M1.count.xz contains 2212 rows and merges to 2212 rows.
## preprocessing/08/outputs/bowtie_sagalactiae_a909/trimmed_ca-v0M1.count.xz contains 2212 rows and merges to 2212 rows.
## preprocessing/09/outputs/bowtie_sagalactiae_a909/trimmed_ca-v0M1.count.xz contains 2212 rows and merges to 2212 rows.
## Finished reading count data.
## Matched 2048 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 2207 rows and 9 columns.
write_expt(a909_expt,
a909_written <-excel=glue::glue("excel/{rundate}-a909_counts-v{ver}.xlsx"))
## Writing the first sheet, containing a legend and some summary data.
## Writing the raw reads.
## Graphing the raw reads.
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete
## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete
## Attempting mixed linear model with: ~ (1|condition) + (1|batch)
## Fitting the expressionset to the model, this is slow.
## Dividing work into 100 chunks...
## Error in serialize(data, node$con, xdr = FALSE) :
## error writing to connection
## A couple of common errors:
## An error like 'vtv downdated' may be because there are too many 0s, filter the data and rerun.
## An error like 'number of levels of each grouping factor must be < number of observations' means
## that the factor used is not appropriate for the analysis - it really only works for factors
## which are shared among multiple samples.
## Retrying with only condition in the model.
## Loading required package: Matrix
##
## Total:10 s
## Placing factor: condition at the beginning of the model.
## Writing the normalized reads.
## Graphing the normalized reads.
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete
## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete
## Attempting mixed linear model with: ~ (1|condition) + (1|batch)
## Fitting the expressionset to the model, this is slow.
## Dividing work into 100 chunks...
## Error in serialize(data, node$con, xdr = FALSE) :
## error writing to connection
## A couple of common errors:
## An error like 'vtv downdated' may be because there are too many 0s, filter the data and rerun.
## An error like 'number of levels of each grouping factor must be < number of observations' means
## that the factor used is not appropriate for the analysis - it really only works for factors
## which are shared among multiple samples.
## Retrying with only condition in the model.
##
## Total:9 s
## Placing factor: condition at the beginning of the model.
## Writing the median reads by factor.
create_expt(metadata="sample_sheets/sagalactiae_samples.xlsx",
cjb111_expt <-gene_info=cjb111_annot, file_column="cjb111_filename")
## Reading the sample metadata.
## The sample definitions comprises: 9 rows(samples) and 30 columns(metadata fields).
## Reading count tables.
## Reading count files with read.table().
## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/tnseq/sagalacticae_2019/preprocessing/01/outputs/bowtie_sagalactiae_cjb111/trimmed_ca-v0M1.count.xz contains 2213 rows.
## preprocessing/02/outputs/bowtie_sagalactiae_cjb111/trimmed_ca-v0M1.count.xz contains 2213 rows and merges to 2213 rows.
## preprocessing/03/outputs/bowtie_sagalactiae_cjb111/trimmed_ca-v0M1.count.xz contains 2213 rows and merges to 2213 rows.
## preprocessing/04/outputs/bowtie_sagalactiae_cjb111/trimmed_ca-v0M1.count.xz contains 2213 rows and merges to 2213 rows.
## preprocessing/05/outputs/bowtie_sagalactiae_cjb111/trimmed_ca-v0M1.count.xz contains 2213 rows and merges to 2213 rows.
## preprocessing/06/outputs/bowtie_sagalactiae_cjb111/trimmed_ca-v0M1.count.xz contains 2213 rows and merges to 2213 rows.
## preprocessing/07/outputs/bowtie_sagalactiae_cjb111/trimmed_ca-v0M1.count.xz contains 2213 rows and merges to 2213 rows.
## preprocessing/08/outputs/bowtie_sagalactiae_cjb111/trimmed_ca-v0M1.count.xz contains 2213 rows and merges to 2213 rows.
## preprocessing/09/outputs/bowtie_sagalactiae_cjb111/trimmed_ca-v0M1.count.xz contains 2213 rows and merges to 2213 rows.
## Finished reading count data.
## Matched 2208 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 2208 rows and 9 columns.
write_expt(cjb111_expt,
cjb111_written <-excel=glue::glue("excel/{rundate}-cjb111_counts-v{ver}.xlsx"))
## Writing the first sheet, containing a legend and some summary data.
## Writing the raw reads.
## Graphing the raw reads.
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete
## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete
## Attempting mixed linear model with: ~ (1|condition) + (1|batch)
## Fitting the expressionset to the model, this is slow.
## Dividing work into 100 chunks...
## Error in serialize(data, node$con, xdr = FALSE) :
## error writing to connection
## A couple of common errors:
## An error like 'vtv downdated' may be because there are too many 0s, filter the data and rerun.
## An error like 'number of levels of each grouping factor must be < number of observations' means
## that the factor used is not appropriate for the analysis - it really only works for factors
## which are shared among multiple samples.
## Retrying with only condition in the model.
##
## Total:10 s
## Placing factor: condition at the beginning of the model.
## Writing the normalized reads.
## Graphing the normalized reads.
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete
## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete
## Attempting mixed linear model with: ~ (1|condition) + (1|batch)
## Fitting the expressionset to the model, this is slow.
## Dividing work into 100 chunks...
## Error in serialize(data, node$con, xdr = FALSE) :
## error writing to connection
## A couple of common errors:
## An error like 'vtv downdated' may be because there are too many 0s, filter the data and rerun.
## An error like 'number of levels of each grouping factor must be < number of observations' means
## that the factor used is not appropriate for the analysis - it really only works for factors
## which are shared among multiple samples.
## Retrying with only condition in the model.
##
## Total:10 s
## Placing factor: condition at the beginning of the model.
## Writing the median reads by factor.
create_expt(metadata="sample_sheets/sagalactiae_samples.xlsx",
vr2603_expt <-gene_info=vr2603_annot, file_column="vr2603_filename")
## Reading the sample metadata.
## The sample definitions comprises: 9 rows(samples) and 30 columns(metadata fields).
## Reading count tables.
## Reading count files with read.table().
## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/tnseq/sagalacticae_2019/preprocessing/01/outputs/bowtie_sagalactiae_2603vr/trimmed_ca-v0M1.count.xz contains 2284 rows.
## preprocessing/02/outputs/bowtie_sagalactiae_2603vr/trimmed_ca-v0M1.count.xz contains 2284 rows and merges to 2284 rows.
## preprocessing/03/outputs/bowtie_sagalactiae_2603vr/trimmed_ca-v0M1.count.xz contains 2284 rows and merges to 2284 rows.
## preprocessing/04/outputs/bowtie_sagalactiae_2603vr/trimmed_ca-v0M1.count.xz contains 2284 rows and merges to 2284 rows.
## preprocessing/05/outputs/bowtie_sagalactiae_2603vr/trimmed_ca-v0M1.count.xz contains 2284 rows and merges to 2284 rows.
## preprocessing/06/outputs/bowtie_sagalactiae_2603vr/trimmed_ca-v0M1.count.xz contains 2284 rows and merges to 2284 rows.
## preprocessing/07/outputs/bowtie_sagalactiae_2603vr/trimmed_ca-v0M1.count.xz contains 2284 rows and merges to 2284 rows.
## preprocessing/08/outputs/bowtie_sagalactiae_2603vr/trimmed_ca-v0M1.count.xz contains 2284 rows and merges to 2284 rows.
## preprocessing/09/outputs/bowtie_sagalactiae_2603vr/trimmed_ca-v0M1.count.xz contains 2284 rows and merges to 2284 rows.
## Finished reading count data.
## Matched 2193 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 2279 rows and 9 columns.
write_expt(vr2603_expt,
vr2603_written <-excel=glue::glue("excel/{rundate}-vr2603_counts-v{ver}.xlsx"))
## Writing the first sheet, containing a legend and some summary data.
## Writing the raw reads.
## Graphing the raw reads.
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete
## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete
## Attempting mixed linear model with: ~ (1|condition) + (1|batch)
## Fitting the expressionset to the model, this is slow.
## Dividing work into 100 chunks...
## Error in serialize(data, node$con, xdr = FALSE) :
## error writing to connection
## A couple of common errors:
## An error like 'vtv downdated' may be because there are too many 0s, filter the data and rerun.
## An error like 'number of levels of each grouping factor must be < number of observations' means
## that the factor used is not appropriate for the analysis - it really only works for factors
## which are shared among multiple samples.
## Retrying with only condition in the model.
##
## Total:11 s
## Placing factor: condition at the beginning of the model.
## Writing the normalized reads.
## Graphing the normalized reads.
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete
## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete
## Attempting mixed linear model with: ~ (1|condition) + (1|batch)
## Fitting the expressionset to the model, this is slow.
## Dividing work into 100 chunks...
## Error in serialize(data, node$con, xdr = FALSE) :
## error writing to connection
## A couple of common errors:
## An error like 'vtv downdated' may be because there are too many 0s, filter the data and rerun.
## An error like 'number of levels of each grouping factor must be < number of observations' means
## that the factor used is not appropriate for the analysis - it really only works for factors
## which are shared among multiple samples.
## Retrying with only condition in the model.
##
## Total:9 s
## Placing factor: condition at the beginning of the model.
## Writing the median reads by factor.
"legend_plot"]] a909_written[[
"raw_libsize"]] a909_written[[
"raw_density"]] a909_written[[
## awesome
"norm_disheat"]] a909_written[[
"norm_corheat"]] a909_written[[
"norm_pca"]] a909_written[[
"raw_libsize"]] cjb111_written[[
"raw_density"]] cjb111_written[[
## awesome
"norm_disheat"]] cjb111_written[[
"norm_corheat"]] cjb111_written[[
"norm_pca"]] cjb111_written[[
"raw_libsize"]] vr2603_written[[
"raw_density"]] vr2603_written[[
## awesome
"norm_disheat"]] vr2603_written[[
"norm_corheat"]] vr2603_written[[
"norm_pca"]] vr2603_written[[
I think this looks reasonable, though it makes me slightly wonder if 04 and 09 are switched. But as long as we are willing to state that the primary difference is between calprotectin and control, then I would suggest against considering it. I think it is reasonable to assume the samples are not switched and this is just how they are. If however, the primary goal is to investigate changing concentrations of calprotectin, then I would want to check into this distribution of samples or make the statement that these two concentrations have no significant difference unless we get more samples to look at.
I moved this above the differential “expression”/“fitness” analysis so that we can add the results from it as annotation data to the DE tables if requested.
01 <- tnseq_saturation(
saturation_"preprocessing/01/outputs/essentiality_sagalactiae_a909/trimmed_ca-v0M1.wig",
adjust=2)
01$plot
saturation_01$hits_summary
saturation_ plot_essentiality(
ess_plts <-"preprocessing/01/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v0M1_gene_tas_m2.csv")
"zbar"]]
ess_plts[[
02 <- tnseq_saturation(
saturation_"preprocessing/02/outputs/essentiality_sagalactiae_a909/trimmed_ca-v0M1.wig")
02$plot
saturation_02$hits_summary
saturation_ plot_essentiality(
ess_plts <-"preprocessing/02/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v0M1_gene_tas_m2.csv")
"zbar"]]
ess_plts[[
03 <- tnseq_saturation(
saturation_"preprocessing/03/outputs/essentiality_sagalactiae_a909/trimmed_ca-v0M1.wig")
03$plot
saturation_03$hits_summary
saturation_ plot_essentiality(
ess_plts <-"preprocessing/03/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v0M1_gene_tas_m2.csv")
"zbar"]]
ess_plts[[
04 <- tnseq_saturation(
saturation_"preprocessing/04/outputs/essentiality_sagalactiae_a909/trimmed_ca-v0M1.wig")
04$plot
saturation_04$hits_summary
saturation_ plot_essentiality(
ess_plts <-"preprocessing/04/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v0M1_gene_tas_m2.csv")
"zbar"]]
ess_plts[[
05 <- tnseq_saturation(
saturation_"preprocessing/05/outputs/essentiality_sagalactiae_a909/trimmed_ca-v0M1.wig")
05$plot
saturation_05$hits_summary
saturation_ plot_essentiality(
ess_plts <-"preprocessing/06/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v0M1_gene_tas_m2.csv")
"zbar"]]
ess_plts[[
06 <- tnseq_saturation(
saturation_"preprocessing/06/outputs/essentiality_sagalactiae_a909/trimmed_ca-v0M1.wig")
06$plot
saturation_06$hits_summary
saturation_ plot_essentiality(
ess_plts <-"preprocessing/06/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v0M1_gene_tas_m2.csv")
"zbar"]]
ess_plts[[
07 <- tnseq_saturation(
saturation_"preprocessing/07/outputs/essentiality_sagalactiae_a909/trimmed_ca-v0M1.wig")
07$plot
saturation_07$hits_summary
saturation_ plot_essentiality(
ess_plts <-"preprocessing/07/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v0M1_gene_tas_m2.csv")
"zbar"]]
ess_plts[[
08 <- tnseq_saturation(
saturation_"preprocessing/08/outputs/essentiality_sagalactiae_a909/trimmed_ca-v0M1.wig")
08$plot
saturation_08$hits_summary
saturation_ plot_essentiality(
ess_plts <-"preprocessing/08/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v0M1_gene_tas_m2.csv")
"zbar"]]
ess_plts[[
09 <- tnseq_saturation(
saturation_"preprocessing/09/outputs/essentiality_sagalactiae_a909/trimmed_ca-v0M1.wig")
09$plot
saturation_09$hits_summary
saturation_ plot_essentiality(
ess_plts <-"preprocessing/09/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v0M1_gene_tas_m2.csv")
"zbar"]]
ess_plts[[
tnseq_saturation(
saturation_control <-"preprocessing/combined_control/outputs/essentiality_sagalactiae_a909/trimmed_ca-v1m1.wig")
$plot
saturation_control plot_essentiality(
ess_plts <-"preprocessing/combined_control/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v1m1_gene_tas_m2.csv")
"zbar"]]
ess_plts[[
sm(tnseq_multi_saturation(meta=pData(a909_expt), meta_column="a909esswig"))
plt <-$ggstats plt
tnseq_multi_saturation(meta=pData(cjb111_expt), meta_column="cjb111esswig") plt <-
## Registered S3 method overwritten by 'broom.mixed':
## method from
## tidy.gamlss broom
## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4
## Registered S3 method overwritten by 'DescTools':
## method from
## reorder.factor gplots
## Warning: Number of labels is greater than default palette color count.
## Try using another color `palette` (and/or `package`).
##
"plot"]] + ggplot2::scale_y_continuous(limits=c(0, 12000)) plt[[
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## I expect the saturation of this to be a bit better than 2603V/R.
tnseq_multi_saturation(meta=pData(vr2603_expt), meta_column="vr2603esswig") plt <-
## Warning: Number of labels is greater than default palette color count.
## Try using another color `palette` (and/or `package`).
##
"plot"]] + ggplot2::scale_y_continuous(limits=c(0, 12000)) plt[[
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## hmm this is not definitive. It looks like cjb111 has more TAs with ~ 1 hit.
## but R/V has more hits with ~ 16<x<64 hits.
tnseq_multi_saturation(meta=pData(a909_expt), meta_column="a909esswig") plt <-
## Warning: Number of labels is greater than default palette color count.
## Try using another color `palette` (and/or `package`).
##
"plot"]] + ggplot2::scale_y_continuous(limits=c(0, 20000)) plt[[
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## I am not quite sure how to interpret this one, we have much more discrete
## numbers of reads than the others.
For differential expression, I am going to assume until I hear otherwise, that my batch assignments are not correct and that the 1,2,3 assignments of the sample names do not actually delineate separate batches. Though, if they do delineate separate batches, it might be taken as a (very)small degree of evidence that 04 and 09 were switched.
list(
combined_essentiality <-"control" = "preprocessing/combined_control/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v1m1_gene_tas_m2.csv",
"low" = "preprocessing/combined_low/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v1m1_gene_tas_m2.csv",
"high" = "preprocessing/combined_high/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v1m1_gene_tas_m2.csv",
"previous" = "preprocessing/combined_previous/outputs/essentiality_sagalactiae_a909/mh_ess-combined_ca_ta-v0M1_gene_tas_m2.csv")
data.frame()
ess_table <-for (f in 1:length(combined_essentiality)) {
names(combined_essentiality)[f]
name <- c("orf", "k", "n", "r", "s", "zbar", "call")
column_names <- paste0(name, "_", column_names)
names <- readr::read_tsv(combined_essentiality[[f]], comment="#", col_names=names)
r <-colnames(r)[1] <- "orf"
if (f == 1) {
r
ess_table <-else {
} merge(ess_table, r, by="orf")
ess_table <-
} }
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## control_orf = col_character(),
## control_k = col_double(),
## control_n = col_double(),
## control_r = col_double(),
## control_s = col_double(),
## control_zbar = col_double(),
## control_call = col_character()
## )
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## low_orf = col_character(),
## low_k = col_double(),
## low_n = col_double(),
## low_r = col_double(),
## low_s = col_double(),
## low_zbar = col_double(),
## low_call = col_character()
## )
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## high_orf = col_character(),
## high_k = col_double(),
## high_n = col_double(),
## high_r = col_double(),
## high_s = col_double(),
## high_zbar = col_double(),
## high_call = col_character()
## )
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## previous_orf = col_character(),
## previous_k = col_double(),
## previous_n = col_double(),
## previous_r = col_double(),
## previous_s = col_double(),
## previous_zbar = col_double(),
## previous_call = col_character()
## )
rownames(ess_table) <- gsub(x=ess_table[["orf"]], pattern="^cds_", replacement="")
"orf"]] <- NULL
ess_table[[
all_pairwise(a909_expt, model_batch=FALSE, parallel=FALSE) a909_de <-
## Plotting a PCA before surrogate/batch inclusion.
## Not putting labels on the PC plot.
## Assuming no batch in model for testing pca.
## Not putting labels on the PC plot.
## Starting basic_pairwise().
## Starting basic pairwise comparison.
## Basic step 0/3: Filtering data.
## Basic step 0/3: Normalizing data.
## Basic step 0/3: Converting data.
## Basic step 0/3: Transforming data.
## Basic step 1/3: Creating mean and variance tables.
## Basic step 2/3: Performing 6 comparisons.
## Basic step 3/3: Creating faux DE Tables.
## Basic: Returning tables.
## Starting deseq_pairwise().
## Starting DESeq2 pairwise comparisons.
## The data should be suitable for EdgeR/DESeq/EBSeq. If they freak out, check the state of the count table and ensure that it is in integer counts.
## Choosing the non-intercept containing model.
## DESeq2 step 1/5: Including only condition in the deseq model.
## converting counts to integer mode
## it appears that the last variable in the design formula, 'condition',
## has a factor level, 'control', which is not the reference level. we recommend
## to use factor(...,levels=...) or relevel() to set this as the reference level
## before proceeding. for more information, please see the 'Note on factor levels'
## in vignette('DESeq2').
## it appears that the last variable in the design formula, 'condition',
## has a factor level, 'control', which is not the reference level. we recommend
## to use factor(...,levels=...) or relevel() to set this as the reference level
## before proceeding. for more information, please see the 'Note on factor levels'
## in vignette('DESeq2').
## DESeq2 step 2/5: Estimate size factors.
## DESeq2 step 3/5: Estimate dispersions.
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## Using a parametric fitting seems to have worked.
## DESeq2 step 4/5: nbinomWaldTest.
## Starting ebseq_pairwise().
## The data should be suitable for EdgeR/DESeq/EBSeq. If they freak out, check the state of the count table and ensure that it is in integer counts.
## Starting EBSeq pairwise subset.
## Choosing the non-intercept containing model.
## Starting EBTest of cal_high vs. cal_low.
## Copying ppee values as ajusted p-values until I figure out how to deal with them.
## Starting EBTest of cal_high vs. control.
## Copying ppee values as ajusted p-values until I figure out how to deal with them.
## Starting EBTest of cal_low vs. control.
## Copying ppee values as ajusted p-values until I figure out how to deal with them.
## Starting edger_pairwise().
## Starting edgeR pairwise comparisons.
## The data should be suitable for EdgeR/DESeq/EBSeq. If they freak out, check the state of the count table and ensure that it is in integer counts.
## Choosing the non-intercept containing model.
## EdgeR step 1/9: Importing and normalizing data.
## EdgeR step 2/9: Estimating the common dispersion.
## EdgeR step 3/9: Estimating dispersion across genes.
## EdgeR step 4/9: Estimating GLM Common dispersion.
## EdgeR step 5/9: Estimating GLM Trended dispersion.
## EdgeR step 6/9: Estimating GLM Tagged dispersion.
## EdgeR step 7/9: Running glmFit, switch to glmQLFit by changing the argument 'edger_test'.
## EdgeR step 8/9: Making pairwise contrasts.
## Starting limma_pairwise().
## Starting limma pairwise comparison.
## libsize was not specified, this parameter has profound effects on limma's result.
## Using the libsize from expt$libsize.
## Limma step 1/6: choosing model.
## Choosing the non-intercept containing model.
## Limma step 2/6: running limma::voom(), switch with the argument 'which_voom'.
## Using normalize.method=quantile for voom.
## Limma step 3/6: running lmFit with method: ls.
## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)
## Limma step 5/6: Running eBayes with robust=FALSE and trend=FALSE.
## Limma step 6/6: Writing limma outputs.
## Limma step 6/6: 1/3: Creating table: cal_low_vs_cal_high. Adjust=BH
## Limma step 6/6: 2/3: Creating table: control_vs_cal_high. Adjust=BH
## Limma step 6/6: 3/3: Creating table: control_vs_cal_low. Adjust=BH
## Limma step 6/6: 1/3: Creating table: cal_high. Adjust=BH
## Limma step 6/6: 2/3: Creating table: cal_low. Adjust=BH
## Limma step 6/6: 3/3: Creating table: control. Adjust=BH
## Comparing analyses.
## Used reverse contrast for deseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts
## = extra_contrasts): The merge of limma, cal_low_vs_cal_high and deseq,
## cal_low_vs_cal_high failed.
## Used reverse contrast for deseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts
## = extra_contrasts): The merge of limma, control_vs_cal_high and deseq,
## control_vs_cal_high failed.
## Used reverse contrast for deseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts
## = extra_contrasts): The merge of limma, control_vs_cal_low and deseq,
## control_vs_cal_low failed.
## Used reverse contrast for edger.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, cal_low_vs_cal_high=cal_low-cal_high, and
## edger, cal_low_vs_cal_high=cal_low-cal_high, failed.
## Used reverse contrast for edger.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_high=control-cal_high, and
## edger, control_vs_cal_high=control-cal_high, failed.
## Used reverse contrast for edger.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_low=control-cal_low, and
## edger, control_vs_cal_low=control-cal_low, failed.
## Used reverse contrast for ebseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, cal_low_vs_cal_high=cal_low-cal_high, and
## ebseq, cal_low_vs_cal_high=cal_low-cal_high, failed.
## Used reverse contrast for ebseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_high=control-cal_high, and
## ebseq, control_vs_cal_high=control-cal_high, failed.
## Used reverse contrast for ebseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_low=control-cal_low, and
## ebseq, control_vs_cal_low=control-cal_low, failed.
## Used reverse contrast for basic.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, cal_low_vs_cal_high=cal_low-cal_high, and
## basic, cal_low_vs_cal_high=cal_low-cal_high, failed.
## Used reverse contrast for basic.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_high=control-cal_high, and
## basic, control_vs_cal_high=control-cal_high, failed.
## Used reverse contrast for basic.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_low=control-cal_low, and
## basic, control_vs_cal_low=control-cal_low, failed.
list(
a909_contrasts <-"low_vs_control" = c("cal_low", "control"),
"high_vs_control" = c("cal_high", "control"))
combine_de_tables(
a909_tables <-extra_annot=ess_table,
keepers=a909_contrasts,
a909_de, excel=glue::glue("excel/{rundate}-a909_tables-v{ver}.xlsx"))
## Writing a legend of columns.
## Printing a pca plot before/after surrogates/batch estimation.
## Working on 1/2: low_vs_control which is: cal_low/control.
## Found inverse table with control_vs_cal_low
## Used the inverse table, might need to -1 the logFC and stat.
## Warning in combine_single_de_table(li = limma, ed = edger, eb = ebseq, de =
## deseq, : The deseq table seems to be missing.
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Unable to find the table in the set of possible tables.
## The possible tables are: cal_low_vs_cal_high=cal_low-cal_high,, control_vs_cal_high=control-cal_high,, control_vs_cal_low=control-cal_low,
## Working on 2/2: high_vs_control which is: cal_high/control.
## Found inverse table with control_vs_cal_high
## Used the inverse table, might need to -1 the logFC and stat.
## Warning in combine_single_de_table(li = limma, ed = edger, eb = ebseq, de =
## deseq, : The deseq table seems to be missing.
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Unable to find the table in the set of possible tables.
## The possible tables are: cal_low_vs_cal_high=cal_low-cal_high,, control_vs_cal_high=control-cal_high,, control_vs_cal_low=control-cal_low,
## Adding venn plots for low_vs_control.
## Limma expression coefficients for low_vs_control; R^2: 0.994; equation: y = 0.996x + 0.0272
## Deseq expression coefficients for low_vs_control; R^2: 0.987; equation: y = 0.985x + 0.167
## Edger expression coefficients for low_vs_control; R^2: 0.994; equation: y = 0.996x + 0.0363
## Adding venn plots for high_vs_control.
## Limma expression coefficients for high_vs_control; R^2: 0.994; equation: y = 0.999x - 0.0103
## Deseq expression coefficients for high_vs_control; R^2: 0.988; equation: y = 0.993x + 0.0786
## Edger expression coefficients for high_vs_control; R^2: 0.995; equation: y = 0.997x + 0.0218
## Writing summary information, compare_plot is: TRUE.
## Performing save of excel/20201015-a909_tables-v20191105.xlsx.
##a909_sig <- extract_significant_genes(
## a909_tables,
## excel=glue::glue("excel/a909_sig-v{ver}.xlsx"))
all_pairwise(vr2603_expt, model_batch=FALSE, parallel=FALSE) vr2603_de <-
## Plotting a PCA before surrogate/batch inclusion.
## Not putting labels on the PC plot.
## Assuming no batch in model for testing pca.
## Not putting labels on the PC plot.
## Starting basic_pairwise().
## Starting basic pairwise comparison.
## Basic step 0/3: Filtering data.
## Basic step 0/3: Normalizing data.
## Basic step 0/3: Converting data.
## Basic step 0/3: Transforming data.
## Basic step 1/3: Creating mean and variance tables.
## Basic step 2/3: Performing 6 comparisons.
## Basic step 3/3: Creating faux DE Tables.
## Basic: Returning tables.
## Starting deseq_pairwise().
## Starting DESeq2 pairwise comparisons.
## The data should be suitable for EdgeR/DESeq/EBSeq. If they freak out, check the state of the count table and ensure that it is in integer counts.
## Choosing the non-intercept containing model.
## DESeq2 step 1/5: Including only condition in the deseq model.
## converting counts to integer mode
## it appears that the last variable in the design formula, 'condition',
## has a factor level, 'control', which is not the reference level. we recommend
## to use factor(...,levels=...) or relevel() to set this as the reference level
## before proceeding. for more information, please see the 'Note on factor levels'
## in vignette('DESeq2').
## it appears that the last variable in the design formula, 'condition',
## has a factor level, 'control', which is not the reference level. we recommend
## to use factor(...,levels=...) or relevel() to set this as the reference level
## before proceeding. for more information, please see the 'Note on factor levels'
## in vignette('DESeq2').
## DESeq2 step 2/5: Estimate size factors.
## DESeq2 step 3/5: Estimate dispersions.
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## Using a parametric fitting seems to have worked.
## DESeq2 step 4/5: nbinomWaldTest.
## Starting ebseq_pairwise().
## The data should be suitable for EdgeR/DESeq/EBSeq. If they freak out, check the state of the count table and ensure that it is in integer counts.
## Starting EBSeq pairwise subset.
## Choosing the non-intercept containing model.
## Starting EBTest of cal_high vs. cal_low.
## Copying ppee values as ajusted p-values until I figure out how to deal with them.
## Starting EBTest of cal_high vs. control.
## Copying ppee values as ajusted p-values until I figure out how to deal with them.
## Starting EBTest of cal_low vs. control.
## Copying ppee values as ajusted p-values until I figure out how to deal with them.
## Starting edger_pairwise().
## Starting edgeR pairwise comparisons.
## The data should be suitable for EdgeR/DESeq/EBSeq. If they freak out, check the state of the count table and ensure that it is in integer counts.
## Choosing the non-intercept containing model.
## EdgeR step 1/9: Importing and normalizing data.
## EdgeR step 2/9: Estimating the common dispersion.
## EdgeR step 3/9: Estimating dispersion across genes.
## EdgeR step 4/9: Estimating GLM Common dispersion.
## EdgeR step 5/9: Estimating GLM Trended dispersion.
## EdgeR step 6/9: Estimating GLM Tagged dispersion.
## EdgeR step 7/9: Running glmFit, switch to glmQLFit by changing the argument 'edger_test'.
## EdgeR step 8/9: Making pairwise contrasts.
## Starting limma_pairwise().
## Starting limma pairwise comparison.
## libsize was not specified, this parameter has profound effects on limma's result.
## Using the libsize from expt$libsize.
## Limma step 1/6: choosing model.
## Choosing the non-intercept containing model.
## Limma step 2/6: running limma::voom(), switch with the argument 'which_voom'.
## Using normalize.method=quantile for voom.
## Limma step 3/6: running lmFit with method: ls.
## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)
## Limma step 5/6: Running eBayes with robust=FALSE and trend=FALSE.
## Limma step 6/6: Writing limma outputs.
## Limma step 6/6: 1/3: Creating table: cal_low_vs_cal_high. Adjust=BH
## Limma step 6/6: 2/3: Creating table: control_vs_cal_high. Adjust=BH
## Limma step 6/6: 3/3: Creating table: control_vs_cal_low. Adjust=BH
## Limma step 6/6: 1/3: Creating table: cal_high. Adjust=BH
## Limma step 6/6: 2/3: Creating table: cal_low. Adjust=BH
## Limma step 6/6: 3/3: Creating table: control. Adjust=BH
## Comparing analyses.
## Used reverse contrast for deseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts
## = extra_contrasts): The merge of limma, cal_low_vs_cal_high and deseq,
## cal_low_vs_cal_high failed.
## Used reverse contrast for deseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts
## = extra_contrasts): The merge of limma, control_vs_cal_high and deseq,
## control_vs_cal_high failed.
## Used reverse contrast for deseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts
## = extra_contrasts): The merge of limma, control_vs_cal_low and deseq,
## control_vs_cal_low failed.
## Used reverse contrast for edger.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, cal_low_vs_cal_high=cal_low-cal_high, and
## edger, cal_low_vs_cal_high=cal_low-cal_high, failed.
## Used reverse contrast for edger.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_high=control-cal_high, and
## edger, control_vs_cal_high=control-cal_high, failed.
## Used reverse contrast for edger.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_low=control-cal_low, and
## edger, control_vs_cal_low=control-cal_low, failed.
## Used reverse contrast for ebseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, cal_low_vs_cal_high=cal_low-cal_high, and
## ebseq, cal_low_vs_cal_high=cal_low-cal_high, failed.
## Used reverse contrast for ebseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_high=control-cal_high, and
## ebseq, control_vs_cal_high=control-cal_high, failed.
## Used reverse contrast for ebseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_low=control-cal_low, and
## ebseq, control_vs_cal_low=control-cal_low, failed.
## Used reverse contrast for basic.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, cal_low_vs_cal_high=cal_low-cal_high, and
## basic, cal_low_vs_cal_high=cal_low-cal_high, failed.
## Used reverse contrast for basic.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_high=control-cal_high, and
## basic, control_vs_cal_high=control-cal_high, failed.
## Used reverse contrast for basic.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_low=control-cal_low, and
## basic, control_vs_cal_low=control-cal_low, failed.
combine_de_tables(
vr2603_tables <-keepers=a909_contrasts,
vr2603_de, excel=glue::glue("excel/{rundate}-vr2603_tables-v{ver}.xlsx"))
## Writing a legend of columns.
## Printing a pca plot before/after surrogates/batch estimation.
## Working on 1/2: low_vs_control which is: cal_low/control.
## Found inverse table with control_vs_cal_low
## Used the inverse table, might need to -1 the logFC and stat.
## Warning in combine_single_de_table(li = limma, ed = edger, eb = ebseq, de =
## deseq, : The deseq table seems to be missing.
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Unable to find the table in the set of possible tables.
## The possible tables are: cal_low_vs_cal_high=cal_low-cal_high,, control_vs_cal_high=control-cal_high,, control_vs_cal_low=control-cal_low,
## Working on 2/2: high_vs_control which is: cal_high/control.
## Found inverse table with control_vs_cal_high
## Used the inverse table, might need to -1 the logFC and stat.
## Warning in combine_single_de_table(li = limma, ed = edger, eb = ebseq, de =
## deseq, : The deseq table seems to be missing.
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Unable to find the table in the set of possible tables.
## The possible tables are: cal_low_vs_cal_high=cal_low-cal_high,, control_vs_cal_high=control-cal_high,, control_vs_cal_low=control-cal_low,
## Adding venn plots for low_vs_control.
## Limma expression coefficients for low_vs_control; R^2: 0.994; equation: y = 0.997x + 0.0108
## Deseq expression coefficients for low_vs_control; R^2: 0.989; equation: y = 0.998x - 0.00338
## Edger expression coefficients for low_vs_control; R^2: 0.995; equation: y = 0.999x + 0.00155
## Adding venn plots for high_vs_control.
## Limma expression coefficients for high_vs_control; R^2: 0.994; equation: y = 1x - 0.0297
## Deseq expression coefficients for high_vs_control; R^2: 0.99; equation: y = 0.993x + 0.0655
## Edger expression coefficients for high_vs_control; R^2: 0.995; equation: y = 0.999x + 0.00492
## Writing summary information, compare_plot is: TRUE.
## Performing save of excel/20201015-vr2603_tables-v20191105.xlsx.
##vr2630_sig <- extract_significant_genes(
## vr2603_tables,
## excel=glue::glue("excel/{rundate}-vr2603_sig-v{ver}.xlsx"))
all_pairwise(cjb111_expt, model_batch=FALSE, parallel=FALSE) cjb111_de <-
## Plotting a PCA before surrogate/batch inclusion.
## Not putting labels on the PC plot.
## Assuming no batch in model for testing pca.
## Not putting labels on the PC plot.
## Starting basic_pairwise().
## Starting basic pairwise comparison.
## Basic step 0/3: Filtering data.
## Basic step 0/3: Normalizing data.
## Basic step 0/3: Converting data.
## Basic step 0/3: Transforming data.
## Basic step 1/3: Creating mean and variance tables.
## Basic step 2/3: Performing 6 comparisons.
## Basic step 3/3: Creating faux DE Tables.
## Basic: Returning tables.
## Starting deseq_pairwise().
## Starting DESeq2 pairwise comparisons.
## The data should be suitable for EdgeR/DESeq/EBSeq. If they freak out, check the state of the count table and ensure that it is in integer counts.
## Choosing the non-intercept containing model.
## DESeq2 step 1/5: Including only condition in the deseq model.
## converting counts to integer mode
## it appears that the last variable in the design formula, 'condition',
## has a factor level, 'control', which is not the reference level. we recommend
## to use factor(...,levels=...) or relevel() to set this as the reference level
## before proceeding. for more information, please see the 'Note on factor levels'
## in vignette('DESeq2').
## it appears that the last variable in the design formula, 'condition',
## has a factor level, 'control', which is not the reference level. we recommend
## to use factor(...,levels=...) or relevel() to set this as the reference level
## before proceeding. for more information, please see the 'Note on factor levels'
## in vignette('DESeq2').
## DESeq2 step 2/5: Estimate size factors.
## DESeq2 step 3/5: Estimate dispersions.
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## Using a parametric fitting seems to have worked.
## DESeq2 step 4/5: nbinomWaldTest.
## Starting ebseq_pairwise().
## The data should be suitable for EdgeR/DESeq/EBSeq. If they freak out, check the state of the count table and ensure that it is in integer counts.
## Starting EBSeq pairwise subset.
## Choosing the non-intercept containing model.
## Starting EBTest of cal_high vs. cal_low.
## Copying ppee values as ajusted p-values until I figure out how to deal with them.
## Starting EBTest of cal_high vs. control.
## Copying ppee values as ajusted p-values until I figure out how to deal with them.
## Starting EBTest of cal_low vs. control.
## Copying ppee values as ajusted p-values until I figure out how to deal with them.
## Starting edger_pairwise().
## Starting edgeR pairwise comparisons.
## The data should be suitable for EdgeR/DESeq/EBSeq. If they freak out, check the state of the count table and ensure that it is in integer counts.
## Choosing the non-intercept containing model.
## EdgeR step 1/9: Importing and normalizing data.
## EdgeR step 2/9: Estimating the common dispersion.
## EdgeR step 3/9: Estimating dispersion across genes.
## EdgeR step 4/9: Estimating GLM Common dispersion.
## EdgeR step 5/9: Estimating GLM Trended dispersion.
## EdgeR step 6/9: Estimating GLM Tagged dispersion.
## EdgeR step 7/9: Running glmFit, switch to glmQLFit by changing the argument 'edger_test'.
## EdgeR step 8/9: Making pairwise contrasts.
## Starting limma_pairwise().
## Starting limma pairwise comparison.
## libsize was not specified, this parameter has profound effects on limma's result.
## Using the libsize from expt$libsize.
## Limma step 1/6: choosing model.
## Choosing the non-intercept containing model.
## Limma step 2/6: running limma::voom(), switch with the argument 'which_voom'.
## Using normalize.method=quantile for voom.
## Limma step 3/6: running lmFit with method: ls.
## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)
## Limma step 5/6: Running eBayes with robust=FALSE and trend=FALSE.
## Limma step 6/6: Writing limma outputs.
## Limma step 6/6: 1/3: Creating table: cal_low_vs_cal_high. Adjust=BH
## Limma step 6/6: 2/3: Creating table: control_vs_cal_high. Adjust=BH
## Limma step 6/6: 3/3: Creating table: control_vs_cal_low. Adjust=BH
## Limma step 6/6: 1/3: Creating table: cal_high. Adjust=BH
## Limma step 6/6: 2/3: Creating table: cal_low. Adjust=BH
## Limma step 6/6: 3/3: Creating table: control. Adjust=BH
## Comparing analyses.
## Used reverse contrast for deseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts
## = extra_contrasts): The merge of limma, cal_low_vs_cal_high and deseq,
## cal_low_vs_cal_high failed.
## Used reverse contrast for deseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts
## = extra_contrasts): The merge of limma, control_vs_cal_high and deseq,
## control_vs_cal_high failed.
## Used reverse contrast for deseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts
## = extra_contrasts): The merge of limma, control_vs_cal_low and deseq,
## control_vs_cal_low failed.
## Used reverse contrast for edger.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, cal_low_vs_cal_high=cal_low-cal_high, and
## edger, cal_low_vs_cal_high=cal_low-cal_high, failed.
## Used reverse contrast for edger.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_high=control-cal_high, and
## edger, control_vs_cal_high=control-cal_high, failed.
## Used reverse contrast for edger.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_low=control-cal_low, and
## edger, control_vs_cal_low=control-cal_low, failed.
## Used reverse contrast for ebseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, cal_low_vs_cal_high=cal_low-cal_high, and
## ebseq, cal_low_vs_cal_high=cal_low-cal_high, failed.
## Used reverse contrast for ebseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_high=control-cal_high, and
## ebseq, control_vs_cal_high=control-cal_high, failed.
## Used reverse contrast for ebseq.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_low=control-cal_low, and
## ebseq, control_vs_cal_low=control-cal_low, failed.
## Used reverse contrast for basic.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, cal_low_vs_cal_high=cal_low-cal_high, and
## basic, cal_low_vs_cal_high=cal_low-cal_high, failed.
## Used reverse contrast for basic.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_high=control-cal_high, and
## basic, control_vs_cal_high=control-cal_high, failed.
## Used reverse contrast for basic.
## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The merge of deseq, control_vs_cal_low=control-cal_low, and
## basic, control_vs_cal_low=control-cal_low, failed.
combine_de_tables(
cjb111_tables <-keepers=a909_contrasts,
cjb111_de, excel=glue::glue("excel/{rundate}-cjb111_tables-v{ver}.xlsx"))
## Writing a legend of columns.
## Printing a pca plot before/after surrogates/batch estimation.
## Working on 1/2: low_vs_control which is: cal_low/control.
## Found inverse table with control_vs_cal_low
## Used the inverse table, might need to -1 the logFC and stat.
## Warning in combine_single_de_table(li = limma, ed = edger, eb = ebseq, de =
## deseq, : The deseq table seems to be missing.
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Unable to find the table in the set of possible tables.
## The possible tables are: cal_low_vs_cal_high=cal_low-cal_high,, control_vs_cal_high=control-cal_high,, control_vs_cal_low=control-cal_low,
## Working on 2/2: high_vs_control which is: cal_high/control.
## Found inverse table with control_vs_cal_high
## Used the inverse table, might need to -1 the logFC and stat.
## Warning in combine_single_de_table(li = limma, ed = edger, eb = ebseq, de =
## deseq, : The deseq table seems to be missing.
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights)
## Unable to find the table in the set of possible tables.
## The possible tables are: cal_low_vs_cal_high=cal_low-cal_high,, control_vs_cal_high=control-cal_high,, control_vs_cal_low=control-cal_low,
## Adding venn plots for low_vs_control.
## Limma expression coefficients for low_vs_control; R^2: 0.993; equation: y = 0.997x + 0.0131
## Deseq expression coefficients for low_vs_control; R^2: 0.985; equation: y = 0.984x + 0.172
## Edger expression coefficients for low_vs_control; R^2: 0.994; equation: y = 0.996x + 0.0349
## Adding venn plots for high_vs_control.
## Limma expression coefficients for high_vs_control; R^2: 0.993; equation: y = 1x - 0.0172
## Deseq expression coefficients for high_vs_control; R^2: 0.987; equation: y = 0.992x + 0.0875
## Edger expression coefficients for high_vs_control; R^2: 0.994; equation: y = 0.997x + 0.0239
## Writing summary information, compare_plot is: TRUE.
## Performing save of excel/20201015-cjb111_tables-v20191105.xlsx.
##cjb111_sig <- extract_significant_genes(
## a909_tables,
## excel=glue::glue("excel/{rundate}-cjb111_sig-v{ver}.xlsx"))
create_expt("sample_sheets/sagalactiae_combined_samples.xlsx",
combined_expt <-gene_info=a909_annot)
## Reading the sample metadata.
## The sample definitions comprises: 3 rows(samples) and 11 columns(metadata fields).
## Reading count tables.
## Reading count files with read.table().
## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/tnseq/sagalacticae_2019/preprocessing/combined_control/outputs/bowtie_sagalactiae_a909/trimmed_ca-v1m1.count.xz contains 2212 rows.
## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/tnseq/sagalacticae_2019/preprocessing/combined_low/outputs/bowtie_sagalactiae_a909/trimmed_ca-v1m1.count.xz contains 2212 rows and merges to 2212 rows.
## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/tnseq/sagalacticae_2019/preprocessing/combined_high/outputs/bowtie_sagalactiae_a909/trimmed_ca-v1m1.count.xz contains 2212 rows and merges to 2212 rows.
## Finished reading count data.
## Matched 2048 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 2207 rows and 3 columns.
normalize_expt(combined_expt, convert="rpkm",
combined_norm <-transform="log2", na_to_zero=TRUE)
## This function will replace the expt$expressionset slot with:
## log2(rpkm(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
## 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 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 rpkm.
## Step 4: transforming the data with log2.
## transform_counts: Found 436 values equal to 0, adding 1 to the matrix.
## Step 5: not doing batch correction.
exprs(combined_norm)
combined_exprs <-
ess_table[, c("control_call", "low_call", "high_call")] %>%
ess_circos <- dplyr::mutate(control_num = dplyr::case_when(
== "E" ~ 2,
control_call == "NE" ~ 0,
control_call == "S" ~ -1,
control_call == "U" ~ 1),
control_call low_num = dplyr::case_when(
== "E" ~ 2,
low_call == "NE" ~ 0,
low_call == "S" ~ -1,
low_call == "U" ~ 1),
low_call high_num = dplyr::case_when(
== "E" ~ 2,
high_call == "NE" ~ 0,
high_call == "S" ~ -1,
high_call == "U" ~ 1))
high_call "control_call"]] <- as.factor(ess_circos[["control_call"]])
ess_circos[["low_call"]] <- as.factor(ess_circos[["low_call"]])
ess_circos[["high_call"]] <- as.factor(ess_circos[["high_call"]])
ess_circos[[rownames(ess_circos) <- rownames(ess_table)
c("990000", "008800", "000000", "0000AA")
colors <-names(colors) <- c("E", "NE", "S", "U")
a909_tables[["data"]][["low_vs_control"]]
low_df <- a909_tables[["data"]][["high_vs_control"]]
high_df <-
circos_prefix(a909_annot, name="a909") circos_cfg <-
## This assumes you have a colors.conf in circos/colors/ and fonts.conf in circos/fonts/
## It also assumes you have conf/ideogram.conf, conf/ticks.conf, and conf/housekeeping.conf
## It will write circos/conf/a909.conf with a reasonable first approximation config file.
## Wrote karyotype to circos/conf/ideograms/a909.conf
## This should match the karyotype= line in a909.conf
## Wrote ticks to circos/conf/ticks_a909.conf
circos_karyotype(circos_cfg, fasta="reference/sagalactiae_a909.fasta") a909_kary <-
## Wrote karyotype to circos/conf/karyotypes/a909.conf
## This should match the karyotype= line in a909.conf
circos_plus_minus(circos_cfg, width=0.06, thickness=40) a909_plus_minus <-
## Writing data file: circos/data/a909_plus_go.txt with the + strand GO data.
## Writing data file: circos/data/a909_minus_go.txt with the - strand GO data.
## Wrote the +/- config files. Appending their inclusion to the master file.
## Returning the inner width: 0.88. Use it as the outer for the next ring.
circos_hist(circos_cfg, combined_exprs, colname="control",
a909_control_rpkm <-basename="control_rpkm", outer=a909_plus_minus,
fill_color="vvdpblue", width=0.06)
## Writing data file: circos/data/a909_control_rpkmcontrol_hist.txt with the control_rpkmcontrol column.
## Returning the inner width: 0.82. Use it as the outer for the next ring.
circos_hist(circos_cfg, combined_exprs, colname="low",
a909_low_rpkm <-basename="low_rpkm", outer=a909_control_rpkm,
fill_color="vdpblue", width=0.06)
## Writing data file: circos/data/a909_low_rpkmlow_hist.txt with the low_rpkmlow column.
## Returning the inner width: 0.76. Use it as the outer for the next ring.
circos_hist(circos_cfg, combined_exprs, colname="high",
a909_high_rpkm <-basename="high_rpkm", outer=a909_low_rpkm,
fill_color="dpblue", width=0.06)
## Writing data file: circos/data/a909_high_rpkmhigh_hist.txt with the high_rpkmhigh column.
## Returning the inner width: 0.7. Use it as the outer for the next ring.
circos_tile(circos_cfg, ess_circos, colname="control_call",
a909_control_tile <-colors=colors, basename="control_tile",
outer=a909_high_rpkm, width=0.06)
## Writing data file: circos/data/a909control_call_tile.txt with the control_tilecontrol_call column.
## Returning the inner width: 0.64. Use it as the outer for the next ring.
circos_hist(circos_cfg, low_df, colname="deseq_logfc", basename="low",
a909_low <-outer=0.60, fill_color="vvdpgreen", width=0.06, thickness=0.1)
## Writing data file: circos/data/a909_lowdeseq_logfc_hist.txt with the lowdeseq_logfc column.
## Returning the inner width: 0.54. Use it as the outer for the next ring.
circos_hist(circos_cfg, high_df, colname="deseq_logfc", basename="high",
a909_high <-outer=a909_low, fill_color="vdpgreen", width=0.06, thickness=0.1)
## Writing data file: circos/data/a909_highdeseq_logfc_hist.txt with the highdeseq_logfc column.
## Returning the inner width: 0.48. Use it as the outer for the next ring.
circos_suffix(circos_cfg)
a909_suffix <- circos_make(circos_cfg, target="a909") made <-
Yoann sent an email this morning with the following request:
“I have now the table summarizing the data and we would like to summarize it with a Circos plot. I was hoping you could do your amazing magic on this. The data is in the Excel file I am attaching to the email. The figure would look very similar to what you did not too long ago, but with extra data. The idea will be to have the genome information and then going down into the circle we would show the THY data (in column BE), than maybe a little space and we would show the RPMI-based data (BL, BS and BZ), then again a little space and we would have a final circle integrating all previous data for the final call (CI).”
I copied the xls file to the ‘from_yoann/’ directory and exported the relevant portions to a csv file named ‘gbs_essentiality.csv’.
I just realized that the material Yoann sent me is precisely what I already did in the above section which created the ess_circos table.
readr::read_csv("from_yoann/gbs_essentiality.csv") yoann_table <-
## Warning: Missing column names filled in: 'X51' [51], 'X58' [58], 'X65' [65],
## 'X72' [72], 'X79' [79], 'X80' [80], 'X81' [81], 'X82' [82], 'X83' [83],
## 'X84' [84], 'X86' [86], 'X88' [88]
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_character(),
## X51 = col_logical(),
## previousk = col_double(),
## previousn = col_double(),
## previousr = col_double(),
## previouss = col_double(),
## previouszbar = col_double(),
## X58 = col_logical(),
## controlk = col_double(),
## controln = col_double(),
## controlr = col_double(),
## controls = col_double(),
## controlzbar = col_double(),
## X65 = col_logical(),
## lowk = col_double(),
## lown = col_double(),
## lowr = col_double(),
## lows = col_double(),
## lowzbar = col_double(),
## X72 = col_logical(),
## highk = col_double()
## # ... with 10 more columns
## )
## ℹ Use `spec()` for the full column specifications.
## The new column which is desired is actually the 'Final'.
## So, I will merge this with my existing ess_circos data frame and add a color for 'NC'.
## The wanted column names are:
## BE: previouscall BL: controlcall BS: lowcall BZ: highcall CI: 'Final call' (85)
c("controlcall", "Final")
wanted <- as.data.frame(yoann_table[, wanted])
yoann_wanted <-rownames(yoann_wanted) <- yoann_table[["row.names"]]
merge(ess_circos, yoann_wanted, by.x="row.names", by.y="row.names")
yoann_circos <-"Final"]] <- as.factor(yoann_circos[["Final"]])
yoann_circos[[rownames(yoann_circos) <- yoann_circos[["Row.names"]]
"Row.names"]] <- NULL
yoann_circos[[ c("990000", "008800", "444444", "0000AA", "000000")
colors <-names(colors) <- c("E", "NE", "S", "U", "NC")
circos_prefix(a909_annot, name="a909_yoann") yoann_cfg <-
## This assumes you have a colors.conf in circos/colors/ and fonts.conf in circos/fonts/
## It also assumes you have conf/ideogram.conf, conf/ticks.conf, and conf/housekeeping.conf
## It will write circos/conf/a909_yoann.conf with a reasonable first approximation config file.
## Wrote karyotype to circos/conf/ideograms/a909_yoann.conf
## This should match the karyotype= line in a909_yoann.conf
## Wrote ticks to circos/conf/ticks_a909_yoann.conf
circos_karyotype(yoann_cfg, fasta="reference/sagalactiae_a909.fasta") yoann_kary <-
## Wrote karyotype to circos/conf/karyotypes/a909_yoann.conf
## This should match the karyotype= line in a909_yoann.conf
circos_plus_minus(yoann_cfg, width=0.06, thickness=40) yoann_plus_minus <-
## Writing data file: circos/data/a909_yoann_plus_go.txt with the + strand GO data.
## Writing data file: circos/data/a909_yoann_minus_go.txt with the - strand GO data.
## Wrote the +/- config files. Appending their inclusion to the master file.
## Returning the inner width: 0.88. Use it as the outer for the next ring.
circos_tile(yoann_cfg, ess_circos, colname="control_call",
yoann_control_tile <-colors=colors, basename="control_tile",
outer=yoann_plus_minus, width=0.06)
## Writing data file: circos/data/a909_yoanncontrol_call_tile.txt with the control_tilecontrol_call column.
## Returning the inner width: 0.82. Use it as the outer for the next ring.
circos_tile(yoann_cfg, ess_circos, colname="low_call",
yoann_low_tile <-colors=colors, basename="low_tile",
outer=yoann_control_tile, width=0.06)
## Writing data file: circos/data/a909_yoannlow_call_tile.txt with the low_tilelow_call column.
## Returning the inner width: 0.76. Use it as the outer for the next ring.
circos_tile(yoann_cfg, ess_circos, colname="high_call",
yoann_high_tile <-colors=colors, basename="high_tile",
outer=yoann_low_tile, width=0.06)
## Writing data file: circos/data/a909_yoannhigh_call_tile.txt with the high_tilehigh_call column.
## Returning the inner width: 0.7. Use it as the outer for the next ring.
circos_tile(yoann_cfg, yoann_circos, colname="Final",
yoann_final_tile <-colors=colors, basename="final_tile",
outer=yoann_high_tile, width=0.06)
## Writing data file: circos/data/a909_yoannFinal_tile.txt with the final_tileFinal column.
## Returning the inner width: 0.64. Use it as the outer for the next ring.
circos_suffix(yoann_cfg)
yoann_suffix <- circos_make(yoann_cfg, target="a909_yoann") made <-
First, between the COG categories and the THY data, I would like to add the saturation level of the library: this would be the blue lines from the Circos figure you produce a while back (see Circos plot attached to this email).
Second, I would like to add a new ring of essentiality that summarizes the data from the rings (1), (2), (3) and (4) abovementioned. The data is listed in the attached Excel document and is found in column CI, called Lindsey’s call. I would like the data to look a little different so it “pops” as it is a summary. I remember back in the day, we did something similar with some of my data in Group A strep, and we had the non-essential genes showing in yellow.
In order to do the above, I saved the excel file Yoann sent me, exported the relevant column to a 2 column cvs file, containing the gene IDs and calls. Now I will copy/paste the previous set of rings and add the new requests.
202006 <- readr::read_csv("from_yoann/lindsey_calls.csv") yoann_
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## ID = col_character(),
## Final = col_character()
## )
202006 <- as.data.frame(yoann_202006)
yoann_rownames(yoann_202006) <- yoann_202006[["ID"]]
colnames(yoann_202006) <- c("ID202006", "lindsey_call")
c("e31212", "f0c505", "444444", "0000AA", "000000")
lindsey_colors <-names(lindsey_colors) <- c("E", "NE", "S", "U", "NC")
202006 <- merge(yoann_circos, yoann_202006, by.x="row.names", by.y="row.names")
yoann_circos_rownames(yoann_circos_202006) <- yoann_circos_202006[["Row.names"]]
202006[["Row.names"]] <- NULL
yoann_circos_
202006 <- circos_prefix(a909_annot, name="a909_yoann_202006") yoann_cfg_
## This assumes you have a colors.conf in circos/colors/ and fonts.conf in circos/fonts/
## It also assumes you have conf/ideogram.conf, conf/ticks.conf, and conf/housekeeping.conf
## It will write circos/conf/a909_yoann_202006.conf with a reasonable first approximation config file.
## Wrote karyotype to circos/conf/ideograms/a909_yoann_202006.conf
## This should match the karyotype= line in a909_yoann_202006.conf
## Wrote ticks to circos/conf/ticks_a909_yoann_202006.conf
202006 <- circos_karyotype(yoann_cfg_202006, fasta="reference/sagalactiae_a909.fasta") yoann_kary_
## Wrote karyotype to circos/conf/karyotypes/a909_yoann_202006.conf
## This should match the karyotype= line in a909_yoann_202006.conf
202006 <- circos_plus_minus(yoann_cfg_202006, width=0.05, thickness=40) yoann_plus_minus_
## Writing data file: circos/data/a909_yoann_202006_plus_go.txt with the + strand GO data.
## Writing data file: circos/data/a909_yoann_202006_minus_go.txt with the - strand GO data.
## Wrote the +/- config files. Appending their inclusion to the master file.
## Returning the inner width: 0.9. Use it as the outer for the next ring.
202006 <- circos_hist(yoann_cfg_202006, combined_exprs, colname="control",
a909_control_rpkm_basename="control_rpkm", outer=yoann_plus_minus_202006,
fill_color="vvdpblue", width=0.05)
## Writing data file: circos/data/a909_yoann_202006_control_rpkmcontrol_hist.txt with the control_rpkmcontrol column.
## Returning the inner width: 0.85. Use it as the outer for the next ring.
202006 <- circos_hist(yoann_cfg_202006, combined_exprs, colname="low",
a909_low_rpkm_basename="low_rpkm", outer=a909_control_rpkm_202006,
fill_color="vdpblue", width=0.05)
## Writing data file: circos/data/a909_yoann_202006_low_rpkmlow_hist.txt with the low_rpkmlow column.
## Returning the inner width: 0.8. Use it as the outer for the next ring.
202006 <- circos_hist(yoann_cfg_202006, combined_exprs, colname="high",
a909_high_rpkm_basename="high_rpkm", outer=a909_low_rpkm_202006,
fill_color="dpblue", width=0.05)
## Writing data file: circos/data/a909_yoann_202006_high_rpkmhigh_hist.txt with the high_rpkmhigh column.
## Returning the inner width: 0.75. Use it as the outer for the next ring.
202006 <- circos_tile(yoann_cfg_202006, yoann_circos_202006, colname="lindsey_call",
lindsey_tile_colors=lindsey_colors, basename="summary_tile", thickness=60,
stroke_thickness=0.001,
outer=a909_high_rpkm_202006, width=0.05)
## Writing data file: circos/data/a909_yoann_202006lindsey_call_tile.txt with the summary_tilelindsey_call column.
## Returning the inner width: 0.7. Use it as the outer for the next ring.
202006 <- circos_tile(yoann_cfg_202006, ess_circos, colname="control_call",
yoann_control_tile_colors=colors, basename="control_tile", thickness=60,
stroke_thickness=0.001,
outer=lindsey_tile_202006, width=0.05)
## Writing data file: circos/data/a909_yoann_202006control_call_tile.txt with the control_tilecontrol_call column.
## Returning the inner width: 0.65. Use it as the outer for the next ring.
202006 <- circos_tile(yoann_cfg_202006, ess_circos, colname="low_call",
yoann_low_tile_colors=colors, basename="low_tile",
thickness=60, stroke_thickness=0.001,
outer=yoann_control_tile_202006, width=0.05)
## Writing data file: circos/data/a909_yoann_202006low_call_tile.txt with the low_tilelow_call column.
## Returning the inner width: 0.6. Use it as the outer for the next ring.
202006 <- circos_tile(yoann_cfg_202006, ess_circos, colname="high_call",
yoann_high_tile_colors=colors, basename="high_tile",
thickness=60, stroke_thickness=0.001,
outer=yoann_low_tile_202006, width=0.05)
## Writing data file: circos/data/a909_yoann_202006high_call_tile.txt with the high_tilehigh_call column.
## Returning the inner width: 0.55. Use it as the outer for the next ring.
202006 <- circos_tile(yoann_cfg_202006, yoann_circos_202006, colname="Final",
yoann_final_tile_colors=colors, basename="final_tile",
thickness=60, stroke_thickness=0.001,
outer=yoann_high_tile_202006, width=0.05)
## Writing data file: circos/data/a909_yoann_202006Final_tile.txt with the final_tileFinal column.
## Returning the inner width: 0.5. Use it as the outer for the next ring.
202006 <- circos_suffix(yoann_cfg_202006)
yoann_suffix_ circos_make(yoann_cfg_202006, target="a909_yoann_202006") made <-
Here is this morning’s email:
For the information we want to display. - rings 1 and 2, from the outside: the genes by COG, which I think is already what you have there. - ring 3: the mutation saturation, how many Krmit insertions in the initial library in THY, which again I think is what you have. - ring 4: Bayesian essentiality data for the THY medium (this is column BE in the crazy table I sent you this morning). - ring 5: Bayesian essentiality data for the mRPMI medium (this is column BL in the crazy table). - ring 6: Bayesian essentiality data for the mRPMI medium + low Calprotectin (this is column BS in the crazy table). - ring 7: Bayesian essentiality data for the mRPMI medium + high Calprotectin (this is column BZ in the crazy table). - I was wondering if you could include a little space here (maybe it’s too much, don’t spend too much time on that).
Since all the stuff is apparently in the table Yoann called ‘crazy table’, I exported those columns into a csv with that name…
read.csv("from_yoann/crazy_table.csv")
crazy_table <-rownames(crazy_table) <- crazy_table[["row.names"]]
"row.names"]] <- NULL
crazy_table[[
circos_prefix(a909_annot, name="a909_v4") cfgv4 <-
## This assumes you have a colors.conf in circos/colors/ and fonts.conf in circos/fonts/
## It also assumes you have conf/ideogram.conf, conf/ticks.conf, and conf/housekeeping.conf
## It will write circos/conf/a909_v4.conf with a reasonable first approximation config file.
## Wrote karyotype to circos/conf/ideograms/a909_v4.conf
## This should match the karyotype= line in a909_v4.conf
## Wrote ticks to circos/conf/ticks_a909_v4.conf
circos_karyotype(cfgv4, fasta="reference/sagalactiae_a909.fasta") karyv4 <-
## Wrote karyotype to circos/conf/karyotypes/a909_v4.conf
## This should match the karyotype= line in a909_v4.conf
circos_plus_minus(cfgv4, width=0.05, thickness=40) ring12v4 <-
## Writing data file: circos/data/a909_v4_plus_go.txt with the + strand GO data.
## Writing data file: circos/data/a909_v4_minus_go.txt with the - strand GO data.
## Wrote the +/- config files. Appending their inclusion to the master file.
## Returning the inner width: 0.9. Use it as the outer for the next ring.
circos_hist(cfgv4, combined_exprs, colname="control",
ring3v4 <-basename="control_rpkmv4", outer=ring12v4,
fill_color="vvdpblue", width=0.05)
## Writing data file: circos/data/a909_v4_control_rpkmv4control_hist.txt with the control_rpkmv4control column.
## Returning the inner width: 0.85. Use it as the outer for the next ring.
circos_tile(cfgv4, crazy_table, colname="previouscall",
ring4v4 <-colors=lindsey_colors, basename="previousv4", outer=ring3v4,
thickness=60, stroke_thickness=0.001, width=0.05)
## Writing data file: circos/data/a909_v4previouscall_tile.txt with the previousv4previouscall column.
## Returning the inner width: 0.8. Use it as the outer for the next ring.
circos_tile(cfgv4, crazy_table, colname="controlcall",
ring5v4 <-colors=lindsey_colors, basename="controlv4", outer=ring4v4,
thickness=60, stroke_thickness=0.001, width=0.05)
## Writing data file: circos/data/a909_v4controlcall_tile.txt with the controlv4controlcall column.
## Returning the inner width: 0.75. Use it as the outer for the next ring.
circos_tile(cfgv4, crazy_table, colname="lowcall",
ring6v4 <-colors=lindsey_colors, basename="lowv4", outer=ring5v4,
thickness=60, stroke_thickness=0.001, width=0.05)
## Writing data file: circos/data/a909_v4lowcall_tile.txt with the lowv4lowcall column.
## Returning the inner width: 0.7. Use it as the outer for the next ring.
circos_tile(cfgv4, crazy_table, colname="highcall",
ring7v4 <-colors=lindsey_colors, basename="highv4", outer=ring6v4,
thickness=60, stroke_thickness=0.001, width=0.05)
## Writing data file: circos/data/a909_v4highcall_tile.txt with the highv4highcall column.
## Returning the inner width: 0.65. Use it as the outer for the next ring.
circos_suffix(cfgv4)
suffixv4 <- circos_make(cfgv4, target="a909_v4") made <-
Here is today’s email:
“For the Bayesian representations, I was hoping we could display a total of 5 circles, four representing the data obtained for the four different conditions (i.e., THY -column BE-; RPMI -column BL-; RPMI+lowCal -column BS- and RPMI+highCal -column BZ-): for these 4 circles the non-essential genes would be shown in gree. Then, there would be one last circle, that is a summary of the results previously mentioned: this circle would integrate the data from column CI (also called Lindsey Call), and for this circle, I want the non-essential genes represented in yellow, so that we can see that this is a little different than the 4 other circles previously mentioned.”
If I read this correctly, the only thing I did wrong was to get some of the colors wrong?
c("e31212", "12e312", "444444", "0909AA", "090909")
main_colors <-names(main_colors) <- c("E", "NE", "S", "U", "NC")
c("e31212", "f0c505", "444444", "0909AA", "090909")
summary_colors <-names(summary_colors) <- c("E", "NE", "S", "U", "NC")
circos_prefix(a909_annot, name="a909_v5") cfgv5 <-
## This assumes you have a colors.conf in circos/colors/ and fonts.conf in circos/fonts/
## It also assumes you have conf/ideogram.conf, conf/ticks.conf, and conf/housekeeping.conf
## It will write circos/conf/a909_v5.conf with a reasonable first approximation config file.
## Wrote karyotype to circos/conf/ideograms/a909_v5.conf
## This should match the karyotype= line in a909_v5.conf
## Wrote ticks to circos/conf/ticks_a909_v5.conf
circos_karyotype(cfgv5, fasta="reference/sagalactiae_a909.fasta") karyv5 <-
## Wrote karyotype to circos/conf/karyotypes/a909_v5.conf
## This should match the karyotype= line in a909_v5.conf
circos_plus_minus(cfgv5, width=0.07, thickness=40) ring12v5 <-
## Writing data file: circos/data/a909_v5_plus_go.txt with the + strand GO data.
## Writing data file: circos/data/a909_v5_minus_go.txt with the - strand GO data.
## Wrote the +/- config files. Appending their inclusion to the master file.
## Returning the inner width: 0.86. Use it as the outer for the next ring.
circos_hist(cfgv5, combined_exprs, colname="control",
ring3v5 <-basename="control_rpkmv5", outer=ring12v5,
fill_color="vvdpblue", width=0.07)
## Writing data file: circos/data/a909_v5_control_rpkmv5control_hist.txt with the control_rpkmv5control column.
## Returning the inner width: 0.79. Use it as the outer for the next ring.
circos_tile(cfgv5, crazy_table, colname="previouscall",
ring4v5 <-colors=main_colors, basename="previousv5", outer=ring3v5,
thickness=70, stroke_thickness=0.001, width=0.07)
## Writing data file: circos/data/a909_v5previouscall_tile.txt with the previousv5previouscall column.
## Returning the inner width: 0.72. Use it as the outer for the next ring.
circos_tile(cfgv5, crazy_table, colname="controlcall",
ring5v5 <-colors=main_colors, basename="controlv5", outer=ring4v5,
thickness=70, stroke_thickness=0.001, width=0.07)
## Writing data file: circos/data/a909_v5controlcall_tile.txt with the controlv5controlcall column.
## Returning the inner width: 0.65. Use it as the outer for the next ring.
circos_tile(cfgv5, crazy_table, colname="lowcall",
ring6v5 <-colors=main_colors, basename="lowv5", outer=ring5v5,
thickness=70, stroke_thickness=0.001, width=0.07)
## Writing data file: circos/data/a909_v5lowcall_tile.txt with the lowv5lowcall column.
## Returning the inner width: 0.58. Use it as the outer for the next ring.
circos_tile(cfgv5, crazy_table, colname="highcall",
ring7v5 <-colors=main_colors, basename="highv5", outer=ring6v5,
thickness=70, stroke_thickness=0.001, width=0.07)
## Writing data file: circos/data/a909_v5highcall_tile.txt with the highv5highcall column.
## Returning the inner width: 0.51. Use it as the outer for the next ring.
circos_tile(cfgv5, crazy_table, colname="Final",
ring8v5 <-colors=summary_colors, basename="finalv5", outer=ring7v5,
thickness=70, stroke_thickness=0.001, width=0.07)
## Writing data file: circos/data/a909_v5Final_tile.txt with the finalv5Final column.
## Returning the inner width: 0.44. Use it as the outer for the next ring.
circos_suffix(cfgv5)
suffixv5 <- circos_make(cfgv5, target="a909_v5") made <-
create_expt(metadata="sample_sheets/sagalactiae_combined_with_previous_samples.xlsx",
a909_prev <-batch=FALSE, gene_info=a909_annot, file_column="file")
## Reading the sample metadata.
## The sample definitions comprises: 4 rows(samples) and 12 columns(metadata fields).
## Reading count tables.
## Reading count files with read.table().
## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/tnseq/sagalacticae_2019/preprocessing/combined_control/outputs/bowtie_sagalactiae_a909/trimmed_ca-v1m1.count.xz contains 2212 rows.
## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/tnseq/sagalacticae_2019/preprocessing/combined_low/outputs/bowtie_sagalactiae_a909/trimmed_ca-v1m1.count.xz contains 2212 rows and merges to 2212 rows.
## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/tnseq/sagalacticae_2019/preprocessing/combined_high/outputs/bowtie_sagalactiae_a909/trimmed_ca-v1m1.count.xz contains 2212 rows and merges to 2212 rows.
## preprocessing/combined_previous/outputs/bowtie_sagalactiae_a909/combined_ca_ta-v0M1.count_sagalactiae_a909_sno_gene_gene_id.count.xz contains 2212 rows and merges to 2212 rows.
## Finished reading count data.
## Matched 2048 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 2207 rows and 4 columns.
tnseq_saturation(
saturation_control <-"preprocessing/combined_control/outputs/essentiality_sagalactiae_a909/trimmed_ca-v1m1.wig")
$plot saturation_control
## Warning: Removed 16596 rows containing non-finite values (stat_bin).
## Warning: Removed 16596 rows containing non-finite values (stat_density).
## Warning: Removed 2 rows containing missing values (geom_bar).
plot_essentiality(
control_plts <-"preprocessing/combined_control/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v1m1_gene_tas_m2.csv")
"zbar"]] control_plts[[
tnseq_saturation(
saturation_low <-"preprocessing/combined_low/outputs/essentiality_sagalactiae_a909/trimmed_ca-v1m1.wig")
$plot saturation_control
## Warning: Removed 16596 rows containing non-finite values (stat_bin).
## Warning: Removed 16596 rows containing non-finite values (stat_density).
## Warning: Removed 2 rows containing missing values (geom_bar).
plot_essentiality(
low_plts <-"preprocessing/combined_low/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v1m1_gene_tas_m2.csv")
"zbar"]] low_plts[[
tnseq_saturation(
saturation_high <-"preprocessing/combined_high/outputs/essentiality_sagalactiae_a909/trimmed_ca-v1m1.wig")
$plot saturation_control
## Warning: Removed 16596 rows containing non-finite values (stat_bin).
## Warning: Removed 16596 rows containing non-finite values (stat_density).
## Warning: Removed 2 rows containing missing values (geom_bar).
plot_essentiality(
high_plts <-"preprocessing/combined_high/outputs/essentiality_sagalactiae_a909/mh_ess-trimmed_ca-v1m1_gene_tas_m2.csv")
"zbar"]] high_plts[[
tnseq_saturation(
saturation_prev <-"preprocessing/combined_previous/outputs/essentiality_sagalactiae_a909/combined_ca_ta-v0M1.wig")
$plot saturation_control
## Warning: Removed 16596 rows containing non-finite values (stat_bin).
## Warning: Removed 16596 rows containing non-finite values (stat_density).
## Warning: Removed 2 rows containing missing values (geom_bar).
plot_essentiality(
high_plts <-"preprocessing/combined_previous/outputs/essentiality_sagalactiae_a909/mh_ess-combined_ca_ta-v0M1_gene_tas_m2.csv")
"zbar"]] high_plts[[
tnseq_multi_saturation(meta=pData(a909_prev), meta_column="a909esswig")
plt <-$plot plt
$ggstats plt
## notch went outside hinges. Try setting notch=FALSE.
## notch went outside hinges. Try setting notch=FALSE.
## notch went outside hinges. Try setting notch=FALSE.
::pander(sessionInfo()) pander
R version 4.0.3 (2020-10-10)
Platform: x86_64-pc-linux-gnu (64-bit)
locale: LC_CTYPE=en_US.UTF-8, LC_NUMERIC=C, LC_TIME=en_US.UTF-8, LC_COLLATE=en_US.UTF-8, LC_MONETARY=en_US.UTF-8, LC_MESSAGES=en_US.UTF-8, LC_PAPER=en_US.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.UTF-8 and LC_IDENTIFICATION=C
attached base packages: splines, parallel, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: edgeR(v.3.30.3), ruv(v.0.9.7.1), lme4(v.1.1-23), Matrix(v.1.2-18), BiocParallel(v.1.22.0), variancePartition(v.1.18.3), hpgltools(v.1.0), testthat(v.2.3.2), R6(v.2.4.1), Biobase(v.2.48.0) and BiocGenerics(v.0.34.0)
loaded via a namespace (and not attached): corpcor(v.1.6.9), class(v.7.3-17), ps(v.1.4.0), Rsamtools(v.2.4.0), V8(v.3.2.0), lmtest(v.0.9-38), foreach(v.1.5.0), rprojroot(v.1.3-2), crayon(v.1.3.4), MASS(v.7.3-53), PMCMRplus(v.1.5.1), nlme(v.3.1-149), backports(v.1.1.10), metafor(v.2.4-0), ggcorrplot(v.0.1.3), sva(v.3.36.0), GOSemSim(v.2.14.2), rlang(v.0.4.7), XVector(v.0.28.0), readxl(v.1.3.1), performance(v.0.5.0), nloptr(v.1.2.2.2), callr(v.3.4.4), limma(v.3.44.3), bit64(v.4.0.5), loo(v.2.3.1), glue(v.1.4.2), rstan(v.2.21.2), pbkrtest(v.0.4-8.6), processx(v.3.4.4), AnnotationDbi(v.1.50.3), ggstatsplot(v.0.6.1), DOSE(v.3.14.0), haven(v.2.3.1), tidyselect(v.1.1.0), SummarizedExperiment(v.1.18.2), usethis(v.1.6.3), rio(v.0.5.16), XML(v.3.99-0.5), tidyr(v.1.1.2), zoo(v.1.8-8), SuppDists(v.1.1-9.5), GenomicAlignments(v.1.24.0), mc2d(v.0.1-18), xtable(v.1.8-4), MatrixModels(v.0.4-1), magrittr(v.1.5), evaluate(v.0.14), ggplot2(v.3.3.2), cli(v.2.0.2), zlibbioc(v.1.34.0), rstudioapi(v.0.11), miniUI(v.0.1.1.1), fastmatch(v.1.1-0), shiny(v.1.5.0), xfun(v.0.18), inline(v.0.3.16), askpass(v.1.1), parameters(v.0.8.6), pkgbuild(v.1.1.0), bridgesampling(v.1.0-0), caTools(v.1.18.0), tidygraph(v.1.2.0), WRS2(v.1.1-0), Vennerable(v.3.1.0.9000), Brobdingnag(v.1.2-6), expm(v.0.999-5), tibble(v.3.0.3), ggrepel(v.0.8.2), Biostrings(v.2.56.0), reshape(v.0.8.8), zeallot(v.0.1.0), ez(v.4.4-0), withr(v.2.3.0), rcompanion(v.2.3.25), slam(v.0.1-47), bitops(v.1.0-6), ggforce(v.0.3.2), RBGL(v.1.64.0), plyr(v.1.8.6), cellranger(v.1.1.0), e1071(v.1.7-3), coda(v.0.19-4), RcppParallel(v.5.0.2), pillar(v.1.4.6), gplots(v.3.1.0), GenomicFeatures(v.1.40.1), multcomp(v.1.4-14), Rmpfr(v.0.8-1), EBSeq(v.1.28.0), fs(v.1.5.0), europepmc(v.0.4), clusterProfiler(v.3.16.1), paletteer(v.1.2.0), vctrs(v.0.3.4), ellipsis(v.0.3.1), generics(v.0.0.2), nortest(v.1.0-4), urltools(v.1.7.3), devtools(v.2.3.2), tools(v.4.0.3), foreign(v.0.8-80), munsell(v.0.5.0), tweenr(v.1.0.1), fgsea(v.1.14.0), DelayedArray(v.0.14.1), fastmap(v.1.0.1), compiler(v.4.0.3), pkgload(v.1.1.0), abind(v.1.4-5), httpuv(v.1.5.4), rtracklayer(v.1.48.0), blockmodeling(v.1.0.0), sessioninfo(v.1.1.1), DescTools(v.0.99.38), ggExtra(v.0.9), GenomeInfoDbData(v.1.2.3), gridExtra(v.2.3), lattice(v.0.20-41), later(v.1.1.0.1), dplyr(v.1.0.2), prismatic(v.0.2.0), BiocFileCache(v.1.12.1), jsonlite(v.1.7.1), scales(v.1.1.1), graph(v.1.66.0), gld(v.2.6.2), pbapply(v.1.4-3), carData(v.3.0-4), genefilter(v.1.70.0), promises(v.1.1.1), tidyBF(v.0.3.0), car(v.3.0-10), BWStest(v.0.2.2), doParallel(v.1.0.15), R.utils(v.2.10.1), metaBMA(v.0.6.3), effectsize(v.0.3.3), pairwiseComparisons(v.3.0.0), sandwich(v.3.0-0), rmarkdown(v.2.4), openxlsx(v.4.2.2), cowplot(v.1.1.0), statmod(v.1.4.34), Rtsne(v.0.15), ipmisc(v.4.0.0), forcats(v.0.5.0), pander(v.0.6.3), downloader(v.0.4), selectr(v.0.4-2), logspline(v.2.1.16), igraph(v.1.2.6), survival(v.3.2-7), numDeriv(v.2016.8-1.1), yaml(v.2.2.1), metaplus(v.0.7-11), rstantools(v.2.1.1), htmltools(v.0.5.0), memoise(v.1.1.0), fastGHQuad(v.1.0), modeltools(v.0.2-23), locfit(v.1.5-9.4), graphlayouts(v.0.7.0), IRanges(v.2.22.2), quadprog(v.1.5-8), viridisLite(v.0.3.0), gmp(v.0.6-1), digest(v.0.6.25), assertthat(v.0.2.1), mime(v.0.9), rappdirs(v.0.3.1), bayestestR(v.0.7.2), RSQLite(v.2.2.1), LaplacesDemon(v.16.1.4), Exact(v.2.1), remotes(v.2.2.0), data.table(v.1.13.0), blob(v.1.2.1), R.oo(v.1.24.0), lpsymphony(v.1.16.0), S4Vectors(v.0.26.1), preprocessCore(v.1.50.0), labeling(v.0.3), rematch2(v.2.1.2), RCurl(v.1.98-1.2), broom(v.0.7.1), hms(v.0.5.3), colorspace(v.1.4-1), BiocManager(v.1.30.10), GenomicRanges(v.1.40.0), libcoin(v.1.0-6), broom.mixed(v.0.2.6), coin(v.1.3-1), Rcpp(v.1.0.5), mvtnorm(v.1.1-1), enrichplot(v.1.8.1), IHW(v.1.16.0), multcompView(v.0.1-8), fansi(v.0.4.1), grid(v.4.0.3), ggridges(v.0.5.2), lifecycle(v.0.2.0), EMT(v.1.1), StanHeaders(v.2.21.0-6), rootSolve(v.1.8.2.1), statsExpressions(v.0.5.1), zip(v.2.1.1), BayesFactor(v.0.9.12-4.2), curl(v.4.3), ggsignif(v.0.6.0), minqa(v.1.2.4), robustbase(v.0.93-6), broomExtra(v.4.0.6), fastcluster(v.1.1.25), DO.db(v.2.9), PROPER(v.1.20.0), qvalue(v.2.20.0), TH.data(v.1.0-10), desc(v.1.2.0), RColorBrewer(v.1.1-2), iterators(v.1.0.12), TMB(v.1.7.18), stringr(v.1.4.0), directlabels(v.2020.6.17), polyclip(v.1.10-0), triebeard(v.0.3.0), biomaRt(v.2.44.1), purrr(v.0.3.4), gridGraphics(v.0.5-0), rvest(v.0.3.6), mgcv(v.1.8-33), openssl(v.1.4.3), insight(v.0.9.6), lmom(v.2.8), bdsmatrix(v.1.3-4), codetools(v.0.2-16), matrixStats(v.0.57.0), GO.db(v.3.11.4), gtools(v.3.8.2), prettyunits(v.1.1.1), dbplyr(v.1.4.4), R.methodsS3(v.1.8.1), GenomeInfoDb(v.1.24.2), correlation(v.0.4.0), gtable(v.0.3.0), DBI(v.1.1.0), stats4(v.4.0.3), httr(v.1.4.2), KernSmooth(v.2.23-17), stringi(v.1.5.3), kSamples(v.1.2-9), progress(v.1.2.2), reshape2(v.1.4.4), farver(v.2.0.3), ggthemes(v.4.2.0), annotate(v.1.66.0), viridis(v.0.5.1), fdrtool(v.1.2.15), xml2(v.1.3.2), colorRamps(v.2.3), rvcheck(v.0.1.8), bbmle(v.1.0.23.1), boot(v.1.3-25), geneplotter(v.1.66.0), readr(v.1.4.0), ggplotify(v.0.0.5), DEoptimR(v.1.0-8), DESeq2(v.1.28.1), bit(v.4.0.4), scatterpie(v.0.1.5), ggraph(v.2.0.3), pkgconfig(v.2.0.3) and knitr(v.1.30)
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 4f025ebdf7b19ddfef0cf9ddaa9ebe2857477394
## This is hpgltools commit: Wed Aug 19 10:11:52 2020 -0400: 4f025ebdf7b19ddfef0cf9ddaa9ebe2857477394
paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
this_save <-message(paste0("Saving to ", this_save))
## Saving to index-v20191105.rda.xz
sm(saveme(filename=this_save)) tmp <-
loadme(filename=this_save)