thyexpress <- set_expt_colors(thyexpress)
thyexpress_metrics <- sm(graph_metrics(thyexpress))thyexpress_norm <- sm(normalize_expt(thyexpress, transform="log2", convert="cpm", norm="quant", filter=TRUE))
thyexpress_norm_metrics <- sm(graph_metrics(thyexpress_norm))thyexpress_metrics$libsize
thyexpress_metrics$density
thyexpress_norm_metrics$corheat
thyexpress_norm_metrics$pcaplot
thyexpress_norm_metrics$smc
Holy ass crackers, the data is beautiful.
## First try a normal condition + batch model
thyexpress_pairwise <- all_pairwise(thyexpress)## Starting limma pairwise comparison.## libsize was not specified, this parameter has profound effects on limma's result.## Using the libsize from expt$libsize.## Choosing the intercept containing model.## Limma step 1/6: choosing model.## Limma step 2/6: running voom## The voom input was not cpm, converting now.## The voom input was not log2, transforming now.## Limma step 3/6: running lmFit## Limma step 4/6: making and fitting contrasts.## Limma step 5/6: Running eBayes and topTable.## Limma step 6/6: Writing limma outputs.## Limma step 6/6: 1/15: Printing table: rna_el.## Limma step 6/6: 2/15: Printing table: rna_ll.## Limma step 6/6: 3/15: Printing table: rna_st.## Limma step 6/6: 4/15: Printing table: rpoe_comp.## Limma step 6/6: 5/15: Printing table: rpoe_mut.## Limma step 6/6: 6/15: Printing table: rna_ll_vs_rna_el.## Limma step 6/6: 7/15: Printing table: rna_st_vs_rna_el.## Limma step 6/6: 8/15: Printing table: rpoe_comp_vs_rna_el.## Limma step 6/6: 9/15: Printing table: rpoe_mut_vs_rna_el.## Limma step 6/6: 10/15: Printing table: rna_st_vs_rna_ll.## Limma step 6/6: 11/15: Printing table: rpoe_comp_vs_rna_ll.## Limma step 6/6: 12/15: Printing table: rpoe_mut_vs_rna_ll.## Limma step 6/6: 13/15: Printing table: rpoe_comp_vs_rna_st.## Limma step 6/6: 14/15: Printing table: rpoe_mut_vs_rna_st.## Limma step 6/6: 15/15: Printing table: rpoe_mut_vs_rpoe_comp.## Starting DESeq2 pairwise comparisons.## The data should be suitable for EdgeR/DESeq.## If EdgeR/DESeq freaks out, check the state of the count table and ensure that it is in integer counts.## Choosing the intercept containing model.## DESeq2 step 1/5: Including batch and condition in the deseq model.## DESeq2 step 2/5: Estimate size factors.## DESeq2 step 3/5: Estimate dispersions.## DESeq2 step 4/5: nbinomWaldTest.## DESeq2 step 5/5: 1/10: Printing table: rna_ll_vs_rna_el## DESeq2 step 5/5: 2/10: Printing table: rna_st_vs_rna_el## DESeq2 step 5/5: 3/10: Printing table: rpoe_comp_vs_rna_el## DESeq2 step 5/5: 4/10: Printing table: rpoe_mut_vs_rna_el## DESeq2 step 5/5: 5/10: Printing table: rna_st_vs_rna_ll## DESeq2 step 5/5: 6/10: Printing table: rpoe_comp_vs_rna_ll## DESeq2 step 5/5: 7/10: Printing table: rpoe_mut_vs_rna_ll## DESeq2 step 5/5: 8/10: Printing table: rpoe_comp_vs_rna_st## DESeq2 step 5/5: 9/10: Printing table: rpoe_mut_vs_rna_st## DESeq2 step 5/5: 10/10: Printing table: rpoe_mut_vs_rpoe_comp## Collected coefficients for: rna_el## Collected coefficients for: rna_ll## Collected coefficients for: rna_st## Collected coefficients for: rpoe_comp## Collected coefficients for: rpoe_mut## Starting edgeR pairwise comparisons.## The data should be suitable for EdgeR/DESeq.## If EdgeR/DESeq freaks out, check the state of the count table and ensure that it is in integer counts.## Choosing the intercept containing model.## EdgeR step 1/9: normalizing data.## EdgeR step 9/9: 1/10: Printing table: rna_ll_vs_rna_el.## EdgeR step 9/9: 2/10: Printing table: rna_st_vs_rna_el.## EdgeR step 9/9: 3/10: Printing table: rpoe_comp_vs_rna_el.## EdgeR step 9/9: 4/10: Printing table: rpoe_mut_vs_rna_el.## EdgeR step 9/9: 5/10: Printing table: rna_st_vs_rna_ll.## EdgeR step 9/9: 6/10: Printing table: rpoe_comp_vs_rna_ll.## EdgeR step 9/9: 7/10: Printing table: rpoe_mut_vs_rna_ll.## EdgeR step 9/9: 8/10: Printing table: rpoe_comp_vs_rna_st.## EdgeR step 9/9: 9/10: Printing table: rpoe_mut_vs_rna_st.## EdgeR step 9/9: 10/10: Printing table: rpoe_mut_vs_rpoe_comp.## Starting basic pairwise comparison.## This basic pairwise function assumes log2, converted, normalized counts, normalizing now.## Basic step 1/3: Creating median and variance tables.## Basic step 2/3: Performing comparisons.## Basic step 2/3: 1/10: Performing log2 subtraction: rna_ll_vs_rna_el## Basic step 2/3: 2/10: Performing log2 subtraction: rna_st_vs_rna_el## Basic step 2/3: 3/10: Performing log2 subtraction: rpoe_comp_vs_rna_el## Basic step 2/3: 4/10: Performing log2 subtraction: rpoe_mut_vs_rna_el## Basic step 2/3: 5/10: Performing log2 subtraction: rna_st_vs_rna_ll## Basic step 2/3: 6/10: Performing log2 subtraction: rpoe_comp_vs_rna_ll## Basic step 2/3: 7/10: Performing log2 subtraction: rpoe_mut_vs_rna_ll## Basic step 2/3: 8/10: Performing log2 subtraction: rpoe_comp_vs_rna_st## Basic step 2/3: 9/10: Performing log2 subtraction: rpoe_mut_vs_rna_st## Basic step 2/3: 10/10: Performing log2 subtraction: rpoe_mut_vs_rpoe_comp## Basic step 3/3: Creating faux DE Tables.## Basic: Returning tables.## Comparing analyses 1/10: rna_ll_vs_rna_el## Comparing analyses 2/10: rna_st_vs_rna_el## Comparing analyses 3/10: rpoe_comp_vs_rna_el## Comparing analyses 4/10: rpoe_mut_vs_rna_el## Comparing analyses 5/10: rna_st_vs_rna_ll## Comparing analyses 6/10: rpoe_comp_vs_rna_ll## Comparing analyses 7/10: rpoe_mut_vs_rna_ll## Comparing analyses 8/10: rpoe_comp_vs_rna_st## Comparing analyses 9/10: rpoe_mut_vs_rna_st## Comparing analyses 10/10: rpoe_mut_vs_rpoe_comp
keepers <- list(
"ll_el" = c("rna_ll", "rna_el"),
"st_el" = c("rna_st","rna_el"),
"st_ll" = c("rna_st", "rna_ll"))
thyexpress_tables <- combine_de_tables(thyexpress_pairwise,
keepers=keepers,
extra_annot=Biobase::exprs(thyexpress_norm$expressionset),
excel=paste0("excel/thyexpress_tables-", ver, ".xlsx"))## Deleting the file excel/thyexpress_tables-01.xlsx before writing the tables.## Writing a legend of columns.## Working on 1/3: ll_el## Found table with rna_ll_vs_rna_el## Running create_combined_table with inverse=FALSE## The table is: rna_ll_vs_rna_el## Working on 2/3: st_el## Found table with rna_st_vs_rna_el## Running create_combined_table with inverse=FALSE## The table is: rna_st_vs_rna_el## Working on 3/3: st_ll## Found table with rna_st_vs_rna_ll## Running create_combined_table with inverse=FALSE## The table is: rna_st_vs_rna_ll## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## Attempting to add a coefficient plot for ll_el at column 79## Warning: Removed 85 rows containing non-finite values (stat_smooth).## Warning: Removed 85 rows containing missing values (geom_point).## Warning: Removed 45 rows containing missing values (geom_point).## Warning: Removed 22 rows containing missing values (geom_point).## Warning: Removed 85 rows containing missing values (geom_point).## Warning: Removed 7 rows containing non-finite values (stat_smooth).## Warning: Removed 7 rows containing missing values (geom_point).## Warning: Removed 1 rows containing missing values (geom_point).## Warning: Removed 2 rows containing missing values (geom_point).## Warning: Removed 7 rows containing missing values (geom_point).## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## Attempting to add a coefficient plot for st_el at column 79## Warning: Removed 93 rows containing non-finite values (stat_smooth).## Warning: Removed 93 rows containing missing values (geom_point).## Warning: Removed 10 rows containing missing values (geom_point).## Warning: Removed 17 rows containing missing values (geom_point).## Warning: Removed 93 rows containing missing values (geom_point).## Warning: Removed 8 rows containing non-finite values (stat_smooth).## Warning: Removed 8 rows containing missing values (geom_point).## Warning: Removed 1 rows containing missing values (geom_point).## Warning: Removed 8 rows containing missing values (geom_point).## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## Attempting to add a coefficient plot for st_ll at column 79## Warning: Removed 100 rows containing non-finite values (stat_smooth).## Warning: Removed 100 rows containing missing values (geom_point).## Warning: Removed 15 rows containing missing values (geom_point).## Warning: Removed 22 rows containing missing values (geom_point).## Warning: Removed 100 rows containing missing values (geom_point).## Warning: Removed 8 rows containing non-finite values (stat_smooth).## Warning: Removed 8 rows containing missing values (geom_point).## Warning: Removed 1 rows containing missing values (geom_point).## Warning: Removed 8 rows containing missing values (geom_point).## Writing summary information.## The sheet pairwise_summary already exists, it will get overwritten## Attempting to add the comparison plot to pairwise_summary at row: 21 and column: 1## Performing save of the workbook.thyexpress_sig <- extract_significant_genes(thyexpress_tables, according_to="all",
excel=paste0("excel/thyexpress_sig-", ver, ".xlsx"))## Writing excel data sheet 1/6## Assuming the fold changes are on the log scale and so taking >< 0## After (adj)p filter, the up genes table has 326 genes.## After (adj)p filter, the down genes table has 744 genes.## Assuming the fold changes are on the log scale and so taking -1 * fc## After fold change filter, the up genes table has 184 genes.## After fold change filter, the down genes table has 309 genes.## Writing excel data sheet 2/6## Assuming the fold changes are on the log scale and so taking >< 0## After (adj)p filter, the up genes table has 283 genes.## After (adj)p filter, the down genes table has 1107 genes.## Assuming the fold changes are on the log scale and so taking -1 * fc## After fold change filter, the up genes table has 206 genes.## After fold change filter, the down genes table has 829 genes.## Writing excel data sheet 3/6## Assuming the fold changes are on the log scale and so taking >< 0## After (adj)p filter, the up genes table has 232 genes.## After (adj)p filter, the down genes table has 726 genes.## Assuming the fold changes are on the log scale and so taking -1 * fc## After fold change filter, the up genes table has 178 genes.## After fold change filter, the down genes table has 556 genes.## Printing significant genes to the file: excel/thyexpress_sig-01.xlsx## Converted change_counts_up to characters.## Converted change_counts_down to characters.## 1/3: Writing excel data sheet up_limma_ll_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 1/3: Writing excel data sheet down_limma_ll_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 2/3: Writing excel data sheet up_limma_st_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 2/3: Writing excel data sheet down_limma_st_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 3/3: Writing excel data sheet up_limma_st_ll## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 3/3: Writing excel data sheet down_limma_st_ll## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## Writing changed genes summary on last sheet.## Writing excel data sheet 4/6## Assuming the fold changes are on the log scale and so taking >< 0## After (adj)p filter, the up genes table has 444 genes.## After (adj)p filter, the down genes table has 429 genes.## Assuming the fold changes are on the log scale and so taking -1 * fc## After fold change filter, the up genes table has 274 genes.## After fold change filter, the down genes table has 139 genes.## Writing excel data sheet 5/6## Assuming the fold changes are on the log scale and so taking >< 0## After (adj)p filter, the up genes table has 578 genes.## After (adj)p filter, the down genes table has 620 genes.## Assuming the fold changes are on the log scale and so taking -1 * fc## After fold change filter, the up genes table has 382 genes.## After fold change filter, the down genes table has 337 genes.## Writing excel data sheet 6/6## Assuming the fold changes are on the log scale and so taking >< 0## After (adj)p filter, the up genes table has 446 genes.## After (adj)p filter, the down genes table has 461 genes.## Assuming the fold changes are on the log scale and so taking -1 * fc## After fold change filter, the up genes table has 312 genes.## After fold change filter, the down genes table has 330 genes.## Printing significant genes to the file: excel/thyexpress_sig-01.xlsx## The sheet number_changed_genes already exists, it will get overwritten## Converted change_counts_up to characters.## Converted change_counts_down to characters.## 1/3: Writing excel data sheet up_edger_ll_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 1/3: Writing excel data sheet down_edger_ll_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## Warning in max(nchar(data[[col]]), na.rm = TRUE): no non-missing arguments to max;
## returning -Inf## 2/3: Writing excel data sheet up_edger_st_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 2/3: Writing excel data sheet down_edger_st_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 3/3: Writing excel data sheet up_edger_st_ll## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 3/3: Writing excel data sheet down_edger_st_ll## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## Writing changed genes summary on last sheet.## Writing excel data sheet 7/6## Assuming the fold changes are on the log scale and so taking >< 0## After (adj)p filter, the up genes table has 464 genes.## After (adj)p filter, the down genes table has 469 genes.## Assuming the fold changes are on the log scale and so taking -1 * fc## After fold change filter, the up genes table has 247 genes.## After fold change filter, the down genes table has 133 genes.## Writing excel data sheet 8/6## Assuming the fold changes are on the log scale and so taking >< 0## After (adj)p filter, the up genes table has 566 genes.## After (adj)p filter, the down genes table has 668 genes.## Assuming the fold changes are on the log scale and so taking -1 * fc## After fold change filter, the up genes table has 367 genes.## After fold change filter, the down genes table has 354 genes.## Writing excel data sheet 9/6## Assuming the fold changes are on the log scale and so taking >< 0## After (adj)p filter, the up genes table has 496 genes.## After (adj)p filter, the down genes table has 533 genes.## Assuming the fold changes are on the log scale and so taking -1 * fc## After fold change filter, the up genes table has 317 genes.## After fold change filter, the down genes table has 350 genes.## Printing significant genes to the file: excel/thyexpress_sig-01.xlsx## The sheet number_changed_genes already exists, it will get overwritten## Converted change_counts_up to characters.## Converted change_counts_down to characters.## 1/3: Writing excel data sheet up_deseq_ll_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 1/3: Writing excel data sheet down_deseq_ll_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## Warning in max(nchar(data[[col]]), na.rm = TRUE): no non-missing arguments to max;
## returning -Inf## 2/3: Writing excel data sheet up_deseq_st_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 2/3: Writing excel data sheet down_deseq_st_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 3/3: Writing excel data sheet up_deseq_st_ll## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 3/3: Writing excel data sheet down_deseq_st_ll## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## Writing changed genes summary on last sheet.## Writing excel data sheet 10/6## Assuming the fold changes are on the log scale and so taking >< 0## After (adj)p filter, the up genes table has 226 genes.## After (adj)p filter, the down genes table has 313 genes.## Assuming the fold changes are on the log scale and so taking -1 * fc## After fold change filter, the up genes table has 116 genes.## After fold change filter, the down genes table has 115 genes.## Writing excel data sheet 11/6## Assuming the fold changes are on the log scale and so taking >< 0## After (adj)p filter, the up genes table has 217 genes.## After (adj)p filter, the down genes table has 292 genes.## Assuming the fold changes are on the log scale and so taking -1 * fc## After fold change filter, the up genes table has 190 genes.## After fold change filter, the down genes table has 216 genes.## Writing excel data sheet 12/6## Assuming the fold changes are on the log scale and so taking >< 0## After (adj)p filter, the up genes table has 259 genes.## After (adj)p filter, the down genes table has 285 genes.## Assuming the fold changes are on the log scale and so taking -1 * fc## After fold change filter, the up genes table has 184 genes.## After fold change filter, the down genes table has 200 genes.## Printing significant genes to the file: excel/thyexpress_sig-01.xlsx## The sheet number_changed_genes already exists, it will get overwritten## Converted change_counts_up to characters.## Converted change_counts_down to characters.## 1/3: Writing excel data sheet up_basic_ll_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 1/3: Writing excel data sheet down_basic_ll_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 2/3: Writing excel data sheet up_basic_st_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 2/3: Writing excel data sheet down_basic_st_el## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 3/3: Writing excel data sheet up_basic_st_ll## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## 3/3: Writing excel data sheet down_basic_st_ll## Converted seqnames to characters.## Converted strand to characters.## Converted dbxref to characters.## Converted genesynonym to characters.## Writing changed genes summary on last sheet.