Introduction to SexSeek
Here we use SexSeek on a few public GEO datasets where
the sex information is available per sample from the authors and fetched
with seqout.
ExampleData()
#> dataset study species assay
#> 1 human_bulk GSE199208 Homo sapiens bulk
#> 2 mouse_bulk GSE85136 Mus musculus bulk
#> 3 mouse_tx GSE232600 Mus musculus bulk, transcript level
#> 4 human_sc GSE196554 Homo sapiens single-cell
#> 5 mouse_sc GSE132420 Mus musculus single-cell
#> 6 rat_bulk GSE100122 Rattus norvegicus bulk
#> 7 zebrafish_bulk GSE123439 Danio rerio bulk
#> tissue
#> 1 submandibular gland
#> 2 nucleus accumbens
#> 3 parathyroid gland
#> 4 skeletal muscle stem cells
#> 5 exorbital lacrimal gland
#> 6 dorsal root ganglion
#> 7 gonadBulk RNA-seq
Human, gene symbols
GSE199208 contains bulk counts from human submandibular gland. The rows of the matrix are gene symbols.
hb <- ExampleData("human_bulk")
dim(hb$counts)
#> [1] 14 15
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")]
#> unit verdict confidence recorded agrees
#> 1 GSM5966687 male high M TRUE
#> 2 GSM5966688 male high M TRUE
#> 3 GSM5966689 male high M TRUE
#> 4 GSM5966690 male high M TRUE
#> 5 GSM5966691 male high M TRUE
#> 6 GSM5966692 male high M TRUE
#> 7 GSM5966693 male high M TRUE
#> 8 GSM5966694 female high F TRUE
#> 9 GSM5966695 female high F TRUE
#> 10 GSM5966696 female high F TRUE
#> 11 GSM5966697 female high F TRUE
#> 12 GSM5966698 male high M TRUE
#> 13 GSM5966699 male high M TRUE
#> 14 GSM5966700 male high M TRUE
#> 15 GSM5966701 female high F TRUE
table(recorded = hb_cmp$recorded_sex, call = hb_cmp$verdict)
#> call
#> recorded female male
#> female 5 0
#> male 0 10
Mouse, Ensembl gene IDs
GSE85136 contains bulk counts from mouse nucleus accumbens. Each group contains male and female mice. The rows of the matrix are Ensembl gene IDs.
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)
#> call
#> recorded female male
#> female 12 0
#> male 0 12
Mouse, transcript IDs
GSE232600 contains
counts from mouse parathyroid. Salmon calculated these counts for each
transcript ID. Before you use EstimateSex(), use genevintage to add
the transcript counts of each gene.
map <- genevintage::transcript2gene(rownames(tx), species = "Mus musculus")
genes <- genevintage::aggregate_transcripts(tx, mapping = map, unmapped = "keep")mt <- ExampleData("mouse_tx")
head(rownames(mt$counts))
#> [1] "ENSMUST00000189069.7" "ENSMUST00000055032.14" "ENSMUST00000188910.7"
#> [4] "ENSMUST00000186696.7" "ENSMUST00000186726.7" "ENSMUST00000189626.7"
head(mt$tx2gene)
#> input transcript_id gene_id gene_name
#> 1 ENSMUST00000189069.7 ENSMUST00000189069 ENSMUSG00000056673 Kdm5d
#> 2 ENSMUST00000055032.14 ENSMUST00000055032 ENSMUSG00000056673 Kdm5d
#> 3 ENSMUST00000188910.7 ENSMUST00000188910 ENSMUSG00000056673 Kdm5d
#> 4 ENSMUST00000186696.7 ENSMUST00000186696 ENSMUSG00000056673 Kdm5d
#> 5 ENSMUST00000186726.7 ENSMUST00000186726 ENSMUSG00000056673 Kdm5d
#> 6 ENSMUST00000189626.7 ENSMUST00000189626 ENSMUSG00000056673 Kdm5d
genes <- genevintage::aggregate_transcripts(mt$counts,
mapping = mt$tx2gene, label = "id", unmapped = "keep", quiet = TRUE
)
dim(genes)
#> [1] 9 8
mt_cmp <- CompareSex(
EstimateSex(genes, species = mt$species, group = colnames(genes)),
mt$recorded
)
mt_cmp[, c("unit", "verdict", "confidence", "recorded", "flags")]
#> unit verdict confidence recorded
#> 1 GSM7359815 male medium male
#> 2 GSM7359816 male medium male
#> 3 GSM7359817 male medium male
#> 4 GSM7359818 male medium male
#> 5 GSM7359843 uncertain low female
#> 6 GSM7359844 uncertain low female
#> 7 GSM7359845 uncertain low female
#> 8 GSM7359846 uncertain low female
#> flags
#> 1 normalised_input
#> 2 normalised_input
#> 3 normalised_input
#> 4 normalised_input
#> 5 inactivation_marker_missing_from_input;normalised_input
#> 6 inactivation_marker_missing_from_input;normalised_input
#> 7 inactivation_marker_missing_from_input;normalised_input
#> 8 inactivation_marker_missing_from_input;normalised_inputThe data has no Xist transcript. Thus, all females get
the verdict uncertain:
unique(mt_cmp$notes[mt_cmp$verdict == "uncertain"])
#> [1] "No sex-specific signal, and the inactivation marker row is absent from this matrix, so female cannot be called. Re-run with an unfiltered gene set to recover it."
Single-cell RNA-seq
Human
GSE196554 contains
single-cell counts from human skeletal muscle stem cells. Each library
contains the cells of one donor. EstimateSex() gives one
verdict for each library.
hs <- ExampleData("human_sc")
vapply(hs$counts, ncol, 1L) # cells per donor
#> GSM5889068 GSM5889069 GSM5889070 GSM5889071 GSM5889072 GSM5889073 GSM5889074
#> 5206 2394 6897 19993 9847 20868 16881
#> GSM5889075 GSM5889076 GSM5889077 GSM5889078 GSM5889079
#> 2369 3913 4412 4833 882
hs_cmp <- CompareSex(EstimateSex(hs$counts, species = hs$species), hs$recorded)
table(recorded = hs_cmp$recorded_sex, call = hs_cmp$verdict)
#> call
#> recorded female male uncertain
#> female 4 0 2
#> male 0 6 0
If a female has an inact_frac less than 0.05 and no Y
signal, the verdict is uncertain.
hs_cmp[hs_cmp$verdict == "uncertain", c("unit", "inact_score", "inact_frac", "qc_score", "flags")]
#> unit inact_score inact_frac qc_score flags
#> 7 GSM5889074 164.1085 0.02770882 115972026 no_signal
#> 12 GSM5889079 205.3167 0.03257616 6838217 no_signalThis plot shows the fraction of cells that have one or more Y-marker counts:
Mouse
GSE132420 contains single-cell counts from mouse lacrimal gland. The mice are from wild-type strains and from Sjögren’s model strains.
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)
#> call
#> recorded female male
#> female 6 0
#> male 0 6
Species without XIST
s <- SexSpecies()
s[
s$scientific_name %in% c("Homo sapiens", "Mus musculus", "Rattus norvegicus", "Danio rerio"),
c("scientific_name", "system", "inactivation_marker", "status")
]
#> scientific_name system inactivation_marker status
#> 1 Homo sapiens XY XIST validated
#> 2 Mus musculus XY Xist validated
#> 13 Rattus norvegicus XY NA experimental
#> 28 Danio rerio none NA unsupportedRat: males only
GSE100122 contains bulk
counts from rat dorsal root ganglion. The rat panel has no
Xist marker. Thus, EstimateSex() can find
males, but all females get the verdict uncertain.
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")]
#> unit verdict confidence recorded
#> 1 GSM2671558 male medium Male
#> 2 GSM2671559 male medium Male
#> 3 GSM2671560 male medium Male
#> 4 GSM2671561 male medium Male
#> 5 GSM2671562 uncertain low Female
#> 6 GSM2671563 uncertain low Female
#> 7 GSM2671564 uncertain low Female
#> 8 GSM2671565 uncertain low Female
#> flags
#> 1 experimental_panel
#> 2 experimental_panel
#> 3 experimental_panel
#> 4 experimental_panel
#> 5 no_inactivation_marker;experimental_panel
#> 6 no_inactivation_marker;experimental_panel
#> 7 no_inactivation_marker;experimental_panel
#> 8 no_inactivation_marker;experimental_panel
Zebrafish: no panel
GSE123439 contains bulk
counts from zebrafish gonads. SexSeek has no panel for zebrafish. The
fish are from strain AB. This strain does not have a differentiated sex
chromosome. Thus, all verdicts are unknown.
zb <- ExampleData("zebrafish_bulk")
dim(zb$counts)
#> [1] 1 16
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)
#> call
#> recorded unknown
#> female 8
#> male 8
unique(zb_cmp$notes)
#> [1] "Common laboratory strains (AB, TU) have lost the differentiated sex chromosome and determine sex polygenically with a strong environmental component; wild and natural strains DO carry WZ/ZZ sex linkage on chromosome 4. The correct statement is strain-dependent, not that zebrafish lack sex chromosomes - and a counts matrix does not reveal the strain."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)
#> comparison samples agree disagree abstain
#> 1 human bulk 15 15 0 0
#> 2 mouse bulk 24 24 0 0
#> 3 mouse transcripts 8 4 0 4
#> 4 human single-cell 12 10 0 2
#> 5 mouse single-cell 12 12 0 0
#> 6 rat bulk 8 4 0 4
#> 7 zebrafish bulk 16 0 0 16