## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>",
  dev = "png", dpi = 120, fig.retina = 1, out.width = "100%"
)
has_tx <- requireNamespace("genevintage", quietly = TRUE) &&
  "aggregate_transcripts" %in% getNamespaceExports("genevintage")

## ----library------------------------------------------------------------------
library(SexSeek)

## ----datasets-----------------------------------------------------------------
ExampleData()

## ----human-bulk---------------------------------------------------------------
hb <- ExampleData("human_bulk")
dim(hb$counts)

hb_calls <- EstimateSex(hb$counts, species = hb$species, group = colnames(hb$counts))
hb_cmp <- CompareSex(hb_calls, hb$recorded)
hb_cmp[, c("unit", "verdict", "confidence", "recorded", "agrees")]
table(recorded = hb_cmp$recorded_sex, call = hb_cmp$verdict)

## ----human-bulk-plot----------------------------------------------------------
PlotSexEvidence(hb_cmp, main = "GSE199208, human bulk")

## ----mouse-bulk---------------------------------------------------------------
mb <- ExampleData("mouse_bulk")
mb_cmp <- CompareSex(
  EstimateSex(mb$counts, species = mb$species, group = colnames(mb$counts)),
  mb$recorded
)
table(recorded = mb_cmp$recorded_sex, call = mb_cmp$verdict)

## ----mouse-bulk-plot----------------------------------------------------------
PlotSexEvidence(mb_cmp, main = "GSE85136, mouse bulk")

## ----tx-build, eval = FALSE---------------------------------------------------
# map <- genevintage::transcript2gene(rownames(tx), species = "Mus musculus")
# genes <- genevintage::aggregate_transcripts(tx, mapping = map, unmapped = "keep")

## ----tx, eval = has_tx--------------------------------------------------------
mt <- ExampleData("mouse_tx")
head(rownames(mt$counts))
head(mt$tx2gene)

genes <- genevintage::aggregate_transcripts(mt$counts,
  mapping = mt$tx2gene, label = "id", unmapped = "keep", quiet = TRUE
)
dim(genes)

mt_cmp <- CompareSex(
  EstimateSex(genes, species = mt$species, group = colnames(genes)),
  mt$recorded
)
mt_cmp[, c("unit", "verdict", "confidence", "recorded", "flags")]

## ----tx-notes, eval = has_tx--------------------------------------------------
unique(mt_cmp$notes[mt_cmp$verdict == "uncertain"])

## ----tx-plot, eval = has_tx---------------------------------------------------
PlotSexEvidence(mt_cmp, main = "GSE232600, mouse transcript level")

## ----human-sc-----------------------------------------------------------------
hs <- ExampleData("human_sc")
vapply(hs$counts, ncol, 1L) # cells per donor

hs_cmp <- CompareSex(EstimateSex(hs$counts, species = hs$species), hs$recorded)
table(recorded = hs_cmp$recorded_sex, call = hs_cmp$verdict)

## ----human-sc-plot------------------------------------------------------------
PlotSexEvidence(hs_cmp, main = "GSE196554, human single-cell (per donor)")

## ----human-sc-uncertain-------------------------------------------------------
hs_cmp[hs_cmp$verdict == "uncertain", c("unit", "inact_score", "inact_frac", "qc_score", "flags")]

## ----human-sc-cells, fig.height = 4-------------------------------------------
PlotCellFraction(hs$counts, hs$species, hs$recorded, main = "GSE196554, per donor")

## ----mouse-sc-----------------------------------------------------------------
ms <- ExampleData("mouse_sc")
ms_cmp <- CompareSex(EstimateSex(ms$counts, species = ms$species), ms$recorded)
table(recorded = ms_cmp$recorded_sex, call = ms_cmp$verdict)

## ----mouse-sc-plot------------------------------------------------------------
PlotSexEvidence(ms_cmp, main = "GSE132420, mouse single-cell (per mouse)")

## ----mouse-sc-cells, fig.height = 4-------------------------------------------
PlotCellFraction(ms$counts, ms$species, ms$recorded, main = "GSE132420, per mouse")

## ----species-status-----------------------------------------------------------
s <- SexSpecies()
s[
  s$scientific_name %in% c("Homo sapiens", "Mus musculus", "Rattus norvegicus", "Danio rerio"),
  c("scientific_name", "system", "inactivation_marker", "status")
]

## ----rat-bulk-----------------------------------------------------------------
rb <- ExampleData("rat_bulk")
rb_cmp <- CompareSex(
  EstimateSex(rb$counts, species = rb$species, group = colnames(rb$counts)),
  rb$recorded
)
rb_cmp[, c("unit", "verdict", "confidence", "recorded", "flags")]

## ----rat-bulk-plot, fig.height = 4--------------------------------------------
PlotSexEvidence(rb_cmp, main = "GSE100122, rat bulk")

## ----zebrafish----------------------------------------------------------------
zb <- ExampleData("zebrafish_bulk")
dim(zb$counts)
zb_cmp <- CompareSex(
  EstimateSex(zb$counts, species = zb$species, group = colnames(zb$counts)),
  zb$recorded
)
table(recorded = zb_cmp$recorded_sex, call = zb_cmp$verdict)
unique(zb_cmp$notes)

## ----summary------------------------------------------------------------------
cmps <- list("human bulk" = hb_cmp, "mouse bulk" = mb_cmp)
if (has_tx) cmps[["mouse transcripts"]] <- mt_cmp
cmps <- c(cmps, list(
  "human single-cell" = hs_cmp, "mouse single-cell" = ms_cmp,
  "rat bulk" = rb_cmp, "zebrafish bulk" = zb_cmp
))
SummariseComparison(cmps)