Skip to content

Commit

Permalink
Merge fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
kgori committed Jun 1, 2021
1 parent 42851d1 commit 7e4f593
Show file tree
Hide file tree
Showing 5 changed files with 1 addition and 238 deletions.
28 changes: 0 additions & 28 deletions R/sigfit_utility.R
Original file line number Diff line number Diff line change
Expand Up @@ -28,13 +28,6 @@ build_opps_matrix <- function(nsamples, ncat, opps) {
opps <- matrix(rep(human_trinuc_freqs(opps, strand), nsamples),
nrow = nsamples, byrow = TRUE)
}
<<<<<<< HEAD
else {
opps <- as.matrix(opps)
}
# Normalise to sum to nsamples (to avoid very large/small values)
opps / sum(opps) * nsamples
=======
else if (is.numeric(opps) & length(opps) == ncat) {
opps <- matrix(rep(opps, nsamples), nrow = nsamples, byrow = TRUE)
}
Expand All @@ -45,7 +38,6 @@ build_opps_matrix <- function(nsamples, ncat, opps) {
#opps / sum(opps) * nsamples
# Normalise each row to sum to 1
opps / rowSums(opps)
>>>>>>> dev
}
}

Expand Down Expand Up @@ -482,23 +474,12 @@ build_catalogues <- function(variants) {
#' @param signatures Either a numeric matrix of mutational signatures, with one row per signature
#' and one column per mutation type, or a list of matrices generated via \code{\link{retrieve_pars}}.
#' @param opportunities_from Mutational opportunities that are currently imposed on the signatures.
<<<<<<< HEAD
#' This can be a numeric matrix of mutational opportunities, with one row per signature and
#' one column per mutation type (if the signatures were inferred together, then every row should
#' contain the same opportunities). Character values \code{"human-genome"} or \code{"human-exome"}
#' are also admitted, in which case the mutational opportunities of the reference human
#' genome/exome will be used for every signature. If this argument is provided, the signatures will
#' be normalised (divided) by the opportunities. Otherwise (if \code{opportunities_from = NULL}),
#' the signatures will be interpreted as being already normalised (i.e. not relative to any
#' opportunities).
=======
#' This can be a numeric vector of mutational opportunities, with one element per mutation type.
#' Character values \code{"human-genome"} or \code{"human-exome"} are also admitted, in which case
#' the mutational opportunities of the reference human genome/exome will be used.
#' If this argument is provided, the signatures will be normalised (divided) by the opportunities.
#' Otherwise (if \code{opportunities_from = NULL}), the signatures will be interpreted as being
#' already normalised (i.e. not relative to any opportunities).
>>>>>>> dev
#' @param opportunities_to Mutational opportunities that are to be imposed on the signatures.
#' Admits the same values as \code{opportunities_from}. If this argument is provided, the
#' signatures will be multiplied by the opportunities. Otherwise (if \code{opportunities_to = NULL}),
Expand Down Expand Up @@ -536,14 +517,6 @@ convert_signatures <- function(signatures, opportunities_from = NULL, opportunit
stop("Either 'opportunities_from' or 'opportunities_to' must be provided.")
}
signatures <- to_matrix(signatures)
<<<<<<< HEAD
opportunities_from <- build_opps_matrix(nrow(signatures), ncol(signatures), opportunities_from)
opportunities_to <- build_opps_matrix(nrow(signatures), ncol(signatures), opportunities_to)

conv_sigs <- signatures
for (i in 1:nrow(conv_sigs)) {
conv_sigs[i, ] <- conv_sigs[i, ] / opportunities_from[i, ] * opportunities_to[i, ]
=======
opportunities_from <- build_opps_matrix(1, ncol(signatures), opportunities_from)
opportunities_to <- build_opps_matrix(1, ncol(signatures), opportunities_to)
if (nrow(opportunities_from) > 1)
Expand All @@ -554,7 +527,6 @@ convert_signatures <- function(signatures, opportunities_from = NULL, opportunit
conv_sigs <- signatures
for (i in 1:nrow(conv_sigs)) {
conv_sigs[i, ] <- conv_sigs[i, ] / opportunities_from * opportunities_to
>>>>>>> dev
conv_sigs[i, ] <- conv_sigs[i, ] / sum(conv_sigs[i, ])
}
conv_sigs
Expand Down
3 changes: 0 additions & 3 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,6 @@ sigfit is an R package. As it is in early development it is not yet on CRAN, but

devtools::install_github("kgori/sigfit", build_vignettes = TRUE,
build_opts = c("--no-resave-data", "--no-manual"))
<<<<<<< HEAD
=======

The arguments `build_vignettes` and `build_opts` are necessary for the package vignette to be built.

Expand All @@ -41,7 +39,6 @@ For solutions to some of the problems that may arise during installation, see th
Problem:

Error: 'rstan_config' is not an exported object from 'namespace:rstantools'
>>>>>>> dev

The arguments `build_vignettes` and `build_opts` are necessary for the package vignette to be built.

Expand Down
113 changes: 0 additions & 113 deletions doc/sigfit_vignette.R
Original file line number Diff line number Diff line change
@@ -1,22 +1,8 @@
<<<<<<< HEAD
## ----setup, include=FALSE------------------------------------------------
=======
## ----setup, include=FALSE-----------------------------------------------------
>>>>>>> dev
knitr::opts_chunk$set(echo = TRUE)
suppressPackageStartupMessages(library(sigfit))
par(mar = c(6, 4, 6, 4))

<<<<<<< HEAD
## ----devtools_instructions, eval=FALSE-----------------------------------
# devtools::install_github("kgori/sigfit", build_opts = c("--no-resave-data", "--no-manual"))

## ----fetch---------------------------------------------------------------
library(sigfit)
data("cosmic_signatures_v2")

## ----sim-----------------------------------------------------------------
=======
## ----devtools_instructions, eval=FALSE----------------------------------------
# devtools::install_github("kgori/sigfit", build_opts = c("--no-resave-data", "--no-manual"))

Expand All @@ -25,37 +11,24 @@ library(sigfit)
data("cosmic_signatures_v2")

## ----sim----------------------------------------------------------------------
>>>>>>> dev
set.seed(1)
probs <- c(0.4, 0.3, 0.2, 0.1) %*% as.matrix(cosmic_signatures_v2[c(1, 3, 7, 11), ])
mutations <- matrix(rmultinom(1, 20000, probs), nrow = 1)
colnames(mutations) <- colnames(cosmic_signatures_v2)

<<<<<<< HEAD
## ----plotsim, fig.width=20, fig.height=7, out.width="100%", echo=-1------
par(mar = c(4.5,5.5,6.5,1))
plot_spectrum(mutations, name = "Simulated counts")

## ----fitting, warning=FALSE----------------------------------------------
=======
## ----plotsim, fig.width=20, fig.height=7, out.width="100%", echo=-1-----------
par(mar = c(4.5,5.5,6.5,1))
plot_spectrum(mutations, name = "Simulated counts")

## ----fitting, warning=FALSE---------------------------------------------------
>>>>>>> dev
mcmc_samples_fit <- fit_signatures(counts = mutations,
signatures = cosmic_signatures_v2,
iter = 2000,
warmup = 1000,
chains = 1,
seed = 1756)

<<<<<<< HEAD
## ----retrieve_exp--------------------------------------------------------
=======
## ----retrieve_exp-------------------------------------------------------------
>>>>>>> dev
exposures <- retrieve_pars(mcmc_samples_fit,
par = "exposures",
hpd_prob = 0.90)
Expand All @@ -66,11 +39,7 @@ exposures$mean
par(mar=c(8,5,3.5,0))
plot_exposures(mcmc_samples = mcmc_samples_fit)

<<<<<<< HEAD
## ----refitting, eval=FALSE-----------------------------------------------
=======
## ----refitting, eval=FALSE----------------------------------------------------
>>>>>>> dev
# selected_signatures <- c(1, 3, 7, 11)
# mcmc_samples_fit_2 <- fit_signatures(mutations,
# cosmic_signatures_v2[selected_signatures, ],
Expand All @@ -83,35 +52,20 @@ plot_exposures(mcmc_samples = mcmc_samples_fit)
par(mar=c(5,6,6.5,1))
plot_reconstruction(mcmc_samples = mcmc_samples_fit)

<<<<<<< HEAD
## ----plot_all, eval=FALSE------------------------------------------------
=======
## ----plot_all, eval=FALSE-----------------------------------------------------
>>>>>>> dev
# plot_all(mcmc_samples = mcmc_samples_fit,
# out_path = "your/output/dir/here",
# prefix = "Fitting")

<<<<<<< HEAD
## ----load_mutations------------------------------------------------------
=======
## ----load_mutations-----------------------------------------------------------
>>>>>>> dev
library(sigfit)
data("variants_21breast")
head(variants_21breast)

<<<<<<< HEAD
## ----show_samples--------------------------------------------------------
unique(variants_21breast[, 1])

## ----build_catalogues----------------------------------------------------
=======
## ----show_samples-------------------------------------------------------------
unique(variants_21breast[, 1])

## ----build_catalogues---------------------------------------------------------
>>>>>>> dev
counts_21breast <- build_catalogues(variants_21breast)
dim(counts_21breast)
counts_21breast[1:5, 1:8]
Expand All @@ -121,19 +75,12 @@ par(mar = c(5,6,7,2))
par(mfrow = c(7, 3))
plot_spectrum(counts_21breast)

<<<<<<< HEAD
## ----extraction, eval=FALSE----------------------------------------------
=======
## ----extraction, eval=FALSE---------------------------------------------------
>>>>>>> dev
# mcmc_samples_extr <- extract_signatures(counts = counts_21breast,
# nsignatures = 2:7,
# iter = 1000,
# seed = 1756)
<<<<<<< HEAD
=======
# plot_gof(mcmc_samples_extr)
>>>>>>> dev

## ----plot_gof_silent, echo=FALSE, fig.width=9, fig.height=6, out.width="100%"----
## Plot precalculated GOF in order to avoid running the model
Expand All @@ -145,16 +92,6 @@ plot(nS, gof, type = "o", lty = 3, pch = 16, col = "dodgerblue4",
points(nS[best], gof[best], pch = 16, col = "orangered", cex = 1.1)
cat("Estimated best number of signatures:", nS[best], "\n")

<<<<<<< HEAD
## ----extr_names, eval=FALSE----------------------------------------------
# names(mcmc_samples_extr)

## ----extr_names_silent, echo=FALSE---------------------------------------
print(c("nsignatures=1", "nsignatures=2", "nsignatures=3", "nsignatures=4", "nsignatures=5", "nsignatures=6", "nsignatures=7", "best"))
signatures <- extr_signatures

## ----retrieve_sigs, eval=FALSE-------------------------------------------
=======
## ----extr_names, eval=FALSE---------------------------------------------------
# names(mcmc_samples_extr)

Expand All @@ -163,68 +100,43 @@ print(c("nsignatures=1", "nsignatures=2", "nsignatures=3", "nsignatures=4", "nsi
signatures <- extr_signatures

## ----retrieve_sigs, eval=FALSE------------------------------------------------
>>>>>>> dev
# ## Note: mcmc_samples_extr[[N]] contains the extraction results for N signatures
# signatures <- retrieve_pars(mcmc_samples_extr[[4]],
# par = "signatures")

<<<<<<< HEAD
## ----show_signames-------------------------------------------------------
=======
## ----show_signames------------------------------------------------------------
>>>>>>> dev
rownames(signatures$mean)

## ----plot_sigs, warning=FALSE, fig.width=25, fig.height=10, out.width='100%', fig.align="center", echo=-1----
par(mar = c(6,7,6,1))
par(mfrow = c(2, 2))
plot_spectrum(signatures)

<<<<<<< HEAD
## ----match_sigs----------------------------------------------------------
data("cosmic_signatures_v2")
match_signatures(signatures, cosmic_signatures_v2)

## ----refitting2, eval=FALSE----------------------------------------------
=======
## ----match_sigs---------------------------------------------------------------
data("cosmic_signatures_v2")
match_signatures(signatures, cosmic_signatures_v2)

## ----refitting2, eval=FALSE---------------------------------------------------
>>>>>>> dev
# signatures <- retrieve_pars(mcmc_samples_extr_2, "signatures")
# mcmc_samples_refit <- fit_signatures(counts = counts_21breast,
# signatures = signatures,
# iter = 2000,
# warmup = 1000)
# exposures <- retrieve_pars(mcmc_samples_refit, "exposures")

<<<<<<< HEAD
## ----strandwise, fig.width=22, fig.height=7, out.width="100%", echo=-1----
=======
## ----strandwise, fig.width=22, fig.height=7, out.width="100%", echo=-1--------
>>>>>>> dev
par(mar = c(4.5,5.5,6.5,1))
# Load and plot a strand-wise catalogue
data("counts_88liver_strand")
plot_spectrum(counts_88liver_strand[1, ], name = rownames(counts_88liver_strand)[1])

<<<<<<< HEAD
## ----generic, fig.width=22, fig.height=8, out.width="100%", echo=-1------
=======
## ----generic, fig.width=22, fig.height=8, out.width="100%", echo=-1-----------
>>>>>>> dev
par(mar = c(6.5,5.5,6.5,1)); set.seed(0xC0FFEE)
# Create and plot an arbitrary catalogue
counts <- as.numeric(rmultinom(1, 5000, runif(100)))
plot_spectrum(counts, name = "Arbitrary catalogue")

<<<<<<< HEAD
## ----convert_sigs, eval=F------------------------------------------------
=======
## ----convert_sigs, eval=F-----------------------------------------------------
>>>>>>> dev
# # (Following from the signature extraction example presented above)
# # Retrieve genome-derived signatures
# genome_signatures <- retrieve_pars(mcmc_samples_extr[[4]],
Expand All @@ -247,11 +159,7 @@ par(mfrow = c(2, 1), mar = c(5,5.5,6.5,1))
plot_spectrum(signatures$mean[4,], name="Signature D, Genome-relative probabilities")
plot_spectrum(exome_signatures[4,], name="Signature D, Exome-relative probabilities")

<<<<<<< HEAD
## ----multichain, eval=F--------------------------------------------------
=======
## ----multichain, eval=F-------------------------------------------------------
>>>>>>> dev
# # Load mutation counts
# data("counts_21breast")
#
Expand All @@ -276,18 +184,6 @@ plot_spectrum(exome_signatures[4,], name="Signature D, Exome-relative probabilit
# chains = ncores,
# init = inits)

<<<<<<< HEAD
## ------------------------------------------------------------------------
data("cosmic_signatures_v2")
data("cosmic_signatures_v3")
data("cosmic_signatures_v3_strand")

## ------------------------------------------------------------------------
data("variants_21breast")
data("counts_21breast")

## ------------------------------------------------------------------------
=======
## -----------------------------------------------------------------------------
data("cosmic_signatures_v2")
data("cosmic_signatures_v3")
Expand All @@ -299,27 +195,18 @@ plot_spectrum(exome_signatures[4,], name="Signature D, Exome-relative probabilit
data("counts_21breast")

## -----------------------------------------------------------------------------
>>>>>>> dev
data("variants_119breast")
data("variants_119breast_strand")
data("counts_119breast")
data("counts_119breast_strand")

<<<<<<< HEAD
## ------------------------------------------------------------------------
=======
## -----------------------------------------------------------------------------
>>>>>>> dev
data("variants_88liver")
data("variants_88liver_strand")
data("counts_88liver")
data("counts_88liver_strand")

<<<<<<< HEAD
## ------------------------------------------------------------------------
=======
## -----------------------------------------------------------------------------
>>>>>>> dev
data("methylation_50breast")
data("methylation_27normal")

Loading

0 comments on commit 7e4f593

Please sign in to comment.