Examining various adipose ferritin transporter samples.

1 Notes 20241124

1. A discussion regarding the extant paper: I think a mouse must have been missing from one of the figures/analyses. It appears that mouse 3749 was deemed an outlier? In the PCA, there is a gray square directly next to a brown circle; that must be 3749. Najib is arguing to keep that mouse in all analyses. One might similarly argue that the IAD-CL sample may need to be removed.
2. Moving to the current dataset: a figure showing TfR1 ko reduces the mass of the iWAT (white adipose tissue); then showed scRNA of iWAT from a WT mouse which suggests a large population of T-cells as well as separate populations of adipocytes in the process of maturation. There is a large range of expression profiles among the preadipocytes which may contribute to variance observed in the previous/current bulk data.
3. Tamseel Isolated iWAT patch, took single cells from them iWAT patch. The names 1a/2b/etc are the paired sets of cells which were collected. Thus multiple combinations are possible: all controls vs all ko. Paired control/ko comparisons. I performed paired ko/control for the three pairs.
4. A contrast of 1+2ko/1+2control
5. Perhaps I can compare the scRNAseq data and see how well the data compares. Note that our individual cells all came from the same iWAT patch from one mouse. “A single-cell sequence analysis of mouse subcutaneous white adipose tissure reveals dynamic…” The assumption is that our wt samples should well match the preadipocyte cluster. The primary query would be to look at the importance of Fe genes.

1.1 TODO from Najib

1. Add distance+correlation pre/post normalized.
2. Add some text describing the intersection plot of the DESeq2 results.
3. Print out the venn/volcano/ma plots for some contrasts.

2 Introduction

Take a look at these samples. I will spend a few minutes (ok, having started the paper, hours) reading their previous document to fill this intro out.

3 Preprocessing

The following command was used to invoke a combination of trimomatic, fastp, fastqc, hisat2, samtools, htseq, and a couple gatk tools.

cd preprocessing
start=$(pwd)
module add cyoa
for i in $(/bin/ls -d *); do
    cd $i
    cyoa --method prnaseq --species mm39_112 --gff_type gene --gff_tag ID \
         --input $(/bin/ls unprocessed/*.fastq.gz | tr '\n' ':')
    cd $start
done

4 Annotation information

I used the ensembl mm39 release 112. That should work fine with the most recent biomart gene annotations (202411).

For the moment, the useast server seems to be down, I think I will just use the archive until that is fixed.

##mm_annot <- load_biomart_annotations(species = "mmusculus", archive = FALSE)
mm_annot <- load_biomart_annotations(species = "mmusculus")

## The biomart annotations file already exists, loading from it.

mm_genes <- mm_annot[["gene_annotations"]]

5 Extract metadata logs

I will probably disable the following block once I have the starting metadata in a completed state so that if someone wants to see all the numbers of mapped reads etc without having my working tree they will still be available.

sample_sheet <- glue("sample_sheets/TtR1_metadata_modified.xlsx")

meta <- gather_preprocessing_metadata("sample_sheets/all_samples.xlsx", species = "mm39_112",
                                      new_metadata = sample_sheet)

This results in a new column ‘hisat_count_table’, so let us use that, and the gene annotations, to create an expressionset.

6 Create Expressionset

mm_expt <- create_expt(sample_sheet,
                       gene_info = mm_genes,
                       file_column = "hisat_count_table") %>%
  set_expt_conditions(fact = "flox")

## Reading the sample metadata.

## Checking the state of the condition column.

## Checking the state of the batch column.

## Checking the condition factor.

## The sample definitions comprises: 24 rows(samples) and 27 columns(metadata fields).

## Matched 25409 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 25425 features and 24 samples.

## The numbers of samples by condition are:

## 
## control_1a control_2b control_s9 ko_1a_c25b ko_2b_54C1 ko_2b_54C2   ko_s9_s8  so_s9_r14 
##          3          3          3          3          3          3          3          3

mm_expt

## A modified expressionSet containing 25425 genes and 24 sample. There are 28 metadata columns and 13 annotation columns.
## The primary condition is comprised of:
## control_1a, control_2b, control_s9, ko_1a_c25b, ko_2b_54C1, ko_2b_54C2, ko_s9_s8, so_s9_r14.
## Its current state is: raw(data).

combined <- paste0(pData(mm_expt)[["state"]], "_",
                   pData(mm_expt)[["line"]])
pData(mm_expt)[["combined"]] <- combined
plot_meta_sankey(mm_expt, factors = c("state", "line", "batch"))

## A sankey plot describing the metadata of 24 samples,
## including 26 out of 0 nodes and traversing metadata factors:
## .

7 Initial visualization

Plot a few things to make sure there are not samples which are immeidately problematic.

plot_libsize(mm_expt)

## Library sizes of 24 samples, 
## ranging from 14,686,859 to 17,960,838.

plot_nonzero(mm_expt)

## The following samples have less than 16526.25 genes.

##  [1] "1_IATF-1A_" "2_IATF-1A_" "3_IATF-1A_" "4_IATF-2B_" "5_IATF-2B_" "6_IATF-2B_" "7_IATF-S9_"
##  [8] "8_IATF-S9_" "9_IATF-S9_" "10_IATF-1A" "11_IATF-1A" "12_IATF-1A" "13_IATF-2B" "14_IATF-2B"
## [15] "15_IATF-2B" "16_IATF-2B" "19_IATF-S9" "20_IATF-S9" "21_IATF-S9" "22_IATF-S9" "23_IATF-S9"
## [22] "24_IATF-S9"

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
## A non-zero genes plot of 24 samples.
## These samples have an average 16.65 CPM coverage and 15156 genes observed, ranging from 14576 to
## 16635.

plot_disheat(mm_expt)

## A heatmap of pairwise sample distances ranging from: 
## 16297.552607677 to 497679.486643764.

plot_corheat(mm_expt)

## A heatmap of pairwise sample correlations ranging from: 
## 0.573161898662224 to 0.999309770057417.

Nothing jumps out, except it is quite rare for raw data to cluster this tightly, that is pretty weird.

8 Normalized visualization

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

## Removing 13165 low-count genes (12260 remaining).

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

plot_disheat(mm_norm)

## A heatmap of pairwise sample distances ranging from: 
## 14.9259724949041 to 139.831122474879.

plot_corheat(mm_norm)

## A heatmap of pairwise sample correlations ranging from: 
## 0.836561417109657 to 0.998137561607948.

plot_pca(mm_norm)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by control_1a, control_2b, control_s9, ko_1a_c25b, ko_2b_54C1, ko_2b_54C2, ko_s9_s8, so_s9_r14
## Shapes are defined by b1, b2, b3.

Wow, these samples are clustering absurdly well.

9 Check out variance partition

mm_varpart <- simple_varpart(mm_expt, factors = c("state", "line", "batch"))

## Subsetting on features.

## remove_genes_expt(), before removal, there were 18527 genes, now there are 13824.

mm_varpart$partition_plot

mm_expt <- mm_varpart[["modified_expt"]]

10 Perform Differential expression

10.1 All ko/controls

Make the cell status the factor to compare; have the lineage as a secondary factor, perhaps add sva?

mm_state_line <- set_expt_conditions(mm_expt, fact = "state") %>%
  set_expt_batches(fact = "line")

## The numbers of samples by condition are:

## 
## control      ko 
##       9      15

## The number of samples by batch are:

## 
## l1a l2b ls9 
##   6   9   9

keepers <- list(
  "ko_vs_control" = c("ko", "control"))
mm_state_line_de <- all_pairwise(
  mm_state_line, model_batch = TRUE, keepers = keepers,
  keep_underscore = TRUE)

## 
## control      ko 
##       9      15 
## 
## l1a l2b ls9 
##   6   9   9

mm_state_line_de

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: batch in model/limma.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 k_vs_cntrl
## limma_vs_deseq      0.8935
## limma_vs_edger      0.8531
## limma_vs_basic      0.9365
## limma_vs_noiseq     0.7376
## deseq_vs_edger      0.9488
## deseq_vs_basic      0.9381
## deseq_vs_noiseq     0.6486
## edger_vs_basic      0.9392
## edger_vs_noiseq     0.5832
## basic_vs_noiseq     0.8687

mm_state_line_tables <- combine_de_tables(
  mm_state_line_de, label_column = "mgisymbol",
  excel = glue("excel/ko_vs_control_table-v{ver}.xlsx"))

## Deleting the file excel/ko_vs_control_table-v20241206.xlsx before writing the tables.

mm_state_line_tables

## A set of combined differential expression results.

##           table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup limma_sigdown
## 1 ko_vs_control         970          1064        1210          1065         889           703

## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

mm_state_line_sig <- extract_significant_genes(
  mm_state_line_tables, according_to = "deseq",
  excel = glue("excel/ko_vs_control_sig-v{ver}.xlsx"))

## Deleting the file excel/ko_vs_control_sig-v20241206.xlsx before writing the tables.

mm_state_line_sig

## A set of genes deemed significant according to deseq.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##               deseq_up deseq_down
## ko_vs_control      970       1064

10.2 Repeat with sva

all_de_sva <-   all_pairwise(
  mm_state_line, model_batch = "sva",
  keepers = keepers, keep_underscore = TRUE)

## 
## control      ko 
##       9      15 
## Error in density.default(x, adjust = adj) : 'x' contains missing values

## Warning in all_adjusters(input, estimate_type = model_type, surrogates = surrogates): It is highly likely
## that the underlying reason for this error is too many 0's in the dataset, please try doing a filtering of
## the data and retry.

## Warning in all_adjusters(input, estimate_type = model_type, surrogates = surrogates): sva will likely
## fail because some rowSums are 0.

## Error in density.default(x, adjust = adj): 'x' contains missing values

all_de_sva

## Error in eval(expr, envir, enclos): object 'all_de_sva' not found

all_tables_sva <- combine_de_tables(
  all_de_sva,
  excel = glue("excel/ko_vs_control_sva_table-v{ver}.xlsx"))

## Error in eval(expr, envir, enclos): object 'all_de_sva' not found

all_tables_sva

## Error in eval(expr, envir, enclos): object 'all_tables_sva' not found

all_sig_sva <- extract_significant_genes(
  all_tables_sva, excel = glue("excel/three_contrasts_sva_sig-v{ver}.xlsx"))

## Error in eval(expr, envir, enclos): object 'all_tables_sva' not found

all_sig_sva

## Error in eval(expr, envir, enclos): object 'all_sig_sva' not found

10.2.1 Some plots describing these results

mm_state_line_tables[["plots"]][["ko_vs_control"]][["deseq_scatter_plots"]][["scatter"]]

## Warning: Removed 6027 rows containing missing values or values outside the scale range (`geom_point()`).
## Removed 6027 rows containing missing values or values outside the scale range (`geom_point()`).

mm_state_line_tables[["plots"]][["ko_vs_control"]][["deseq_ma_plots"]]

mm_state_line_tables[["plots"]][["ko_vs_control"]][["deseq_vol_plots"]]

10.3 Semi-paired contrasts

The following compares the ko/control for each of the cell groups.

mm_combined <- set_expt_conditions(mm_expt, fact = "combined")

## The numbers of samples by condition are:

## 
## control_l1a control_l2b control_ls9      ko_l1a      ko_l2b      ko_ls9 
##           3           3           3           3           6           6

keepers <- list(
  "l1a" = c("ko_l1a", "control_l1a"),
  "l2b" = c("ko_l2b", "control_l2b"),
  "ls9" = c("ko_ls9", "control_ls9"))
mm_de <- all_pairwise(mm_combined, keepers = keepers, keep_underscore = TRUE)

## 
## control_l1a control_l2b control_ls9      ko_l1a      ko_l2b      ko_ls9 
##           3           3           3           3           6           6 
## 
## b1 b2 b3 
##  8  8  8

mm_de

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: batch in model/limma.

## The primary analysis performed 10 comparisons.

mm_tables <- combine_de_tables(mm_de, excel = "excel/three_contrasts.xlsx")
mm_tables

## A set of combined differential expression results.

##                   table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup limma_sigdown
## 1 ko_l1a_vs_control_l1a        1203          1651        1333          1714        1711          1457
## 2 ko_l2b_vs_control_l2b        1612          1557        1759          1684        1723          1655
## 3 ko_ls9_vs_control_ls9        1089           539        1081           581         752           735

## Plot describing unique/shared genes in a differential expression table.

mm_sig <- extract_significant_genes(mm_tables, excel = "excel/three_contrasts-sig.xlsx")

## Deleting the file excel/three_contrasts-sig.xlsx before writing the tables.

mm_sig

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##                       limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up basic_down
## ko_l1a_vs_control_l1a     1711       1457     1333       1714     1203       1651     1198       1079
## ko_l2b_vs_control_l2b     1723       1655     1759       1684     1612       1557     1261       1194
## ko_ls9_vs_control_ls9      752        735     1081        581     1089        539      536        588

11 Take a look with gprofiler

mm_gp <- all_gprofiler(mm_sig, species = "mmusculus", excel = "excel/all_gp.xlsx")
mm_gp

## Running gProfiler on every set of significant genes found:

##                             BP CORUM HP KEGG MIRNA MF REAC  TF WP
## ko_l1a_vs_control_l1a_up   330     0  0    0     0 63    0 150  0
## ko_l1a_vs_control_l1a_down 882     0  0    0     0 87    5 293  0
## ko_l2b_vs_control_l2b_up   181     0  0    0     0 69    0 106  0
## ko_l2b_vs_control_l2b_down 896     0  1    2     0 89   15 302  1
## ko_ls9_vs_control_ls9_up   517     0  0    4     0 61    5 254  0
## ko_ls9_vs_control_ls9_down 429     0  0    0     2 31    3  21  0

enrichplot::dotplot(mm_gp[["ko_l1a_vs_control_l1a_up"]][["BP_enrich"]])

enrichplot::dotplot(mm_gp[["ko_l1a_vs_control_l1a_down"]][["BP_enrich"]])

enrichplot::dotplot(mm_gp[["ko_l2b_vs_control_l2b_up"]][["BP_enrich"]])

enrichplot::dotplot(mm_gp[["ko_l2b_vs_control_l2b_down"]][["BP_enrich"]])

enrichplot::dotplot(mm_gp[["ko_ls9_vs_control_ls9_up"]][["BP_enrich"]])

enrichplot::dotplot(mm_gp[["ko_ls9_vs_control_ls9_down"]][["BP_enrich"]])

mm_cp <- all_cprofiler(mm_sig, mm_tables, orgdb = "org.Mm.eg.db", excel = "excel/all_cp.xlsx")

## preparing geneSet collections...

## GSEA analysis...

## leading edge analysis...

## done...

## Reading KEGG annotation online: "https://rest.kegg.jp/link/mmu/pathway"...

## Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/mmu"...

## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Deleting the file excel/all_cp.xlsx before writing the tables.
## 
## Writing a sheet containing the legend.
## 
## Writing the BP data.
## 
## Writing the MF data.
## 
## Writing the CC data.
## 
## Writing the KEGG data.
## 
## Finished writing excel file.
## 
## preparing geneSet collections...
## 
## GSEA analysis...
## 
## leading edge analysis...
## 
## done...
## 
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Deleting the file excel/all_cp.xlsx before writing the tables.
## 
## Writing a sheet containing the legend.
## 
## Writing the BP data.
## 
## Writing the MF data.
## 
## Writing the CC data.
## 
## Writing the KEGG data.
## 
## Finished writing excel file.
## 
## preparing geneSet collections...
## 
## GSEA analysis...
## 
## leading edge analysis...
## 
## done...
## 
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Deleting the file excel/all_cp.xlsx before writing the tables.
## 
## Writing a sheet containing the legend.
## 
## Writing the BP data.
## 
## Writing the MF data.
## 
## Writing the CC data.
## 
## Writing the KEGG data.
## 
## Finished writing excel file.
## 
## preparing geneSet collections...
## 
## GSEA analysis...
## 
## leading edge analysis...
## 
## done...
## 
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Deleting the file excel/all_cp.xlsx before writing the tables.
## 
## Writing a sheet containing the legend.
## 
## Writing the BP data.
## 
## Writing the MF data.
## 
## Writing the CC data.
## 
## Writing the KEGG data.
## 
## Finished writing excel file.
## 
## preparing geneSet collections...
## 
## GSEA analysis...
## 
## leading edge analysis...
## 
## done...
## 
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Deleting the file excel/all_cp.xlsx before writing the tables.
## 
## Writing a sheet containing the legend.
## 
## Writing the BP data.
## 
## Writing the MF data.
## 
## Writing the CC data.
## 
## Writing the KEGG data.
## 
## Finished writing excel file.
## 
## preparing geneSet collections...
## 
## GSEA analysis...
## 
## leading edge analysis...
## 
## done...
## 
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Deleting the file excel/all_cp.xlsx before writing the tables.
## 
## Writing a sheet containing the legend.
## 
## Writing the BP data.
## 
## Writing the MF data.
## 
## Writing the CC data.
## 
## Writing the KEGG data.
## 
## Finished writing excel file.

enrichplot::dotplot(mm_cp[["ko_l1a_vs_control_l1a_up"]][["enrich_objects"]][["BP_sig"]])

enrichplot::dotplot(mm_cp[["ko_l1a_vs_control_l1a_down"]][["enrich_objects"]][["BP_sig"]])

topn_gsea <- plot_topn_gsea(mm_cp[["ko_l1a_vs_control_l1a_up"]][["enrich_objects"]][["gse"]])
topn_gsea[[1]]