Skip to content

Commit

Permalink
Merge pull request #86 from BIGslu/dev_msc
Browse files Browse the repository at this point in the history
Replace "." with "_" in parameter names
  • Loading branch information
kdillmcfarland authored May 31, 2023
2 parents eb0cda9 + 8e0ff5b commit f3381fd
Show file tree
Hide file tree
Showing 21 changed files with 533 additions and 342 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
Package: kimma
Type: Package
Title: Kinship In Mixed Model Analysis of RNA-seq
Version: 1.4.4
Version: 1.5.0
Author: Kim Dill-McFarland
Maintainer: Kim Dill-McFarland <kadm@uw.edu>
Description: Linear mixed effects models with pairwise kinship for RNA-seq data.
Expand Down
8 changes: 4 additions & 4 deletions R/extract_lmFit.R
Original file line number Diff line number Diff line change
Expand Up @@ -36,10 +36,10 @@
extract_lmFit <- function(design, fit, contrast_mat=NULL,
dat_genes=NULL, name_genes="geneName",
contrast.mat=NULL, dat.genes=NULL, name.genes=NULL){
#Back compatability
contrast_mat <- contrast.mat
dat_genes <- dat.genes
name_genes <- name.genes
#Back compatibility
if(!is.null(contrast.mat)){contrast_mat <- contrast.mat}
if(!is.null(dat.genes)){dat_genes <- dat.genes}
if(!is.null(name.genes)){name_genes <- name.genes}

#Empty df to hold results
pval.result <- data.frame()
Expand Down
32 changes: 16 additions & 16 deletions R/kimma_cleaning.R
Original file line number Diff line number Diff line change
Expand Up @@ -13,20 +13,20 @@
#' @param weights Matrix or data frame of gene specific weights
#'
#' Subset data (optional)
#' @param subset.var Character list of variable name(s) to filter data by.
#' @param subset.lvl Character list of variable value(s) or level(s) to filter data to. Must match order of subset.var
#' @param subset.genes Character vector of genes to include in models.
#' @param model.lm Character vector of simple linear model version of model provided
#' @param genotype.name Character string. Used internally for kmFit_eQTL
#' @param subset_var Character list of variable name(s) to filter data by.
#' @param subset_lvl Character list of variable value(s) or level(s) to filter data to. Must match order of subset_var
#' @param subset_genes Character vector of genes to include in models.
#' @param model_lm Character vector of simple linear model version of model provided
#' @param genotype_name Character string. Used internally for kmFit_eQTL
#'
#' @return Data frame formatted for use in kmFit
#' @keywords internal
#' @importFrom tidyselect vars_select_helpers

kimma_cleaning <- function(dat=NULL, kin=NULL, patientID="ptID", libraryID="libID",
counts=NULL, meta=NULL, genes=NULL, weights=NULL,
subset.var = NULL, subset.lvl = NULL, subset.genes = NULL,
model.lm = NULL, genotype.name = NULL){
subset_var = NULL, subset_lvl = NULL, subset_genes = NULL,
model_lm = NULL, genotype_name = NULL){
i <- rowname <- NULL
#If data are NOT a voom EList, create a mock version
if(is.null(dat)) {
Expand Down Expand Up @@ -135,10 +135,10 @@ kimma_cleaning <- function(dat=NULL, kin=NULL, patientID="ptID", libraryID="libI
dat.subset <- dat.format

#Subset samples
if(!is.null(subset.var)){
for(i in 1:length(subset.var)) {
if(!is.null(subset_var)){
for(i in 1:length(subset_var)) {
dat.subset$targets <- dplyr::filter(dat.subset$targets,
get(subset.var[[i]]) %in% subset.lvl[[i]])
get(subset_var[[i]]) %in% subset_lvl[[i]])

dat.subset$E <- dplyr::select(as.data.frame(dat.subset$E),
rowname,
Expand All @@ -153,12 +153,12 @@ kimma_cleaning <- function(dat=NULL, kin=NULL, patientID="ptID", libraryID="libI
}

#Subset genes
if(!is.null(subset.genes)){
if(!is.null(subset_genes)){
dat.subset$E <- dplyr::filter(as.data.frame(dat.subset$E),
rowname %in% subset.genes)
rowname %in% subset_genes)
if(!is.null(dat.subset$weights)){
dat.subset$weights <- dplyr::filter(as.data.frame(dat.subset$weights),
rowname %in% subset.genes)
rowname %in% subset_genes)
}
}

Expand Down Expand Up @@ -239,7 +239,7 @@ kimma_cleaning <- function(dat=NULL, kin=NULL, patientID="ptID", libraryID="libI
}

##Missing other variables
all_vars <- gsub("expression~", "", model.lm)
all_vars <- gsub("expression~", "", model_lm)
all_vars <- unique(strsplit(all_vars, "\\+|\\*|:")[[1]])

complete <- to.model %>%
Expand Down Expand Up @@ -286,10 +286,10 @@ kimma_cleaning <- function(dat=NULL, kin=NULL, patientID="ptID", libraryID="libI
}

## Missing data
all_vars <- gsub("expression~", "", model.lm)
all_vars <- gsub("expression~", "", model_lm)
all_vars <- unique(strsplit(all_vars, "\\+|\\*|:")[[1]])

if(!is.null(genotype.name)){ all_vars <- c(all_vars, genotype.name)}
if(!is.null(genotype_name)){ all_vars <- c(all_vars, genotype_name)}

complete <- to.model %>%
dplyr::select(tidyselect::all_of(c(libraryID, all_vars))) %>%
Expand Down
16 changes: 8 additions & 8 deletions R/kimma_lm.R
Original file line number Diff line number Diff line change
@@ -1,28 +1,28 @@
#' Run kimma linear model
#'
#' @param model.lm Character model created in kmFit
#' @param to.model.gene Data frame formatted in kmFit, subset to gene of interest
#' @param model_lm Character model created in kmFit
#' @param to_model_gene Data frame formatted in kmFit, subset to gene of interest
#' @param gene Character of gene to model
#' @param use.weights Logical if gene specific weights should be used in model. Default is FALSE
#' @param use_weights Logical if gene specific weights should be used in model. Default is FALSE
#' @param metrics Logical if should calculate model fit metrics such as AIC, BIC, R-squared. Default is FALSE
#'
#' @return Linear model results data frame for 1 gene
#' @keywords internal

kimma_lm <- function(model.lm, to.model.gene, gene, use.weights, metrics){
kimma_lm <- function(model_lm, to_model_gene, gene, use_weights, metrics){
#Place holder LM results
p.lm <- NaN
sigma.lm <- 0
results.lm <- NULL
.GlobalEnv$to.model.gene <- to.model.gene
.GlobalEnv$to_model_gene <- to_model_gene

#Fit model
if(use.weights){
if(use_weights){
set.seed(42)
fit.lm <- stats::lm(model.lm, data=to.model.gene, weights=to.model.gene$gene_weight)
fit.lm <- stats::lm(model_lm, data=to_model_gene, weights=to_model_gene$gene_weight)
} else{
set.seed(42)
fit.lm <- stats::lm(model.lm, data=to.model.gene, weights=NULL)
fit.lm <- stats::lm(model_lm, data=to_model_gene, weights=NULL)
}

p.lm <- broom::tidy(fit.lm)
Expand Down
26 changes: 13 additions & 13 deletions R/kimma_lme.R
Original file line number Diff line number Diff line change
@@ -1,30 +1,30 @@
#' Run kimma linear mixed effects model
#'
#' @param model.lme Character model created in kmFit
#' @param to.model.gene Data frame formatted in kmFit, subset to gene of interest
#' @param model_lme Character model created in kmFit
#' @param to_model_gene Data frame formatted in kmFit, subset to gene of interest
#' @param gene Character of gene to model
#' @param use.weights Logical if gene specific weights should be used in model. Default is FALSE
#' @param use_weights Logical if gene specific weights should be used in model. Default is FALSE
#' @param metrics Logical if should calculate model fit metrics such as AIC, BIC, R-squared. Default is FALSE
#'
#' @return Linear mixed effect results data frame for 1 gene
#' @keywords internal

kimma_lme <- function(model.lme, to.model.gene, gene, use.weights, metrics){
kimma_lme <- function(model_lme, to_model_gene, gene, use_weights, metrics){
rowname <- NULL
#Place holder LME results
p.lme <- NaN
sigma.lme <- 0
results.lme <- NULL
.GlobalEnv$to.model.gene <- to.model.gene
.GlobalEnv$to_model_gene <- to_model_gene

#Fit LME model
if(use.weights){
if(use_weights){
set.seed(42)
fit.lme <- lme4::lmer(model.lme, data=to.model.gene,
weights=to.model.gene$gene_weight)
fit.lme <- lme4::lmer(model_lme, data=to_model_gene,
weights=to_model_gene$gene_weight)
} else{
set.seed(42)
fit.lme <- lme4::lmer(model.lme, data=to.model.gene, weights=NULL)
fit.lme <- lme4::lmer(model_lme, data=to_model_gene, weights=NULL)
}

#Estimate P-value
Expand Down Expand Up @@ -62,11 +62,11 @@ kimma_lme <- function(model.lme, to.model.gene, gene, use.weights, metrics){

if(metrics){
#Calculate R-squared
if(use.weights){
null <- stats::glm(formula = expression ~ 1, data = to.model.gene,
weights = to.model.gene$gene_weight)
if(use_weights){
null <- stats::glm(formula = expression ~ 1, data = to_model_gene,
weights = to_model_gene$gene_weight)
} else{
null <- stats::glm(formula = expression ~ 1, data = to.model.gene)
null <- stats::glm(formula = expression ~ 1, data = to_model_gene)
}

L0 <- as.vector(stats::logLik(null))
Expand Down
30 changes: 15 additions & 15 deletions R/kimma_lmerel.R
Original file line number Diff line number Diff line change
@@ -1,37 +1,37 @@
#' Run kimma linear mixed effects model with kinship
#'
#' @param model.lme Character model created in kmFit
#' @param to.model.gene Data frame formatted in kmFit, subset to gene of interest
#' @param model_lme Character model created in kmFit
#' @param to_model_gene Data frame formatted in kmFit, subset to gene of interest
#' @param gene Character of gene to model
#' @param kin.subset Pairwise kinship matrix
#' @param use.weights Logical if gene specific weights should be used in model. Default is FALSE
#' @param kin_subset Pairwise kinship matrix
#' @param use_weights Logical if gene specific weights should be used in model. Default is FALSE
#' @param patientID Character of variable name to match dat$targets to kinship row and column names.
#' @param metrics Logical if should calculate model fit metrics such as AIC, BIC, R-squared. Default is FALSE
#'
#' @return Linear mixed effect with kinship results data frame for 1 gene
#' @keywords internal

kimma_lmerel <- function(model.lme, to.model.gene, gene, kin.subset, use.weights,
kimma_lmerel <- function(model_lme, to_model_gene, gene, kin_subset, use_weights,
patientID, metrics){
rowname <- NULL
#Place holder LME results
p.kin <- NaN
sigma.kin <- 0
results.kin <- NULL
.GlobalEnv$to.model.gene <- to.model.gene
.GlobalEnv$to_model_gene <- to_model_gene

#Format kinship matrix to list
kin.ls <- list()
kin.ls[[patientID]] <- as.matrix(kin.subset)
kin.ls[[patientID]] <- as.matrix(kin_subset)
#Fit kin model
if(use.weights){
if(use_weights){
set.seed(42)
fit.kin <- lme4qtl::relmatLmer(model.lme, data=to.model.gene,
fit.kin <- lme4qtl::relmatLmer(model_lme, data=to_model_gene,
relmat = kin.ls,
weights=to.model.gene$gene_weight)
weights=to_model_gene$gene_weight)
} else{
set.seed(42)
fit.kin <- lme4qtl::relmatLmer(model.lme, data=to.model.gene,
fit.kin <- lme4qtl::relmatLmer(model_lme, data=to_model_gene,
relmat = kin.ls)
}

Expand Down Expand Up @@ -70,11 +70,11 @@ kimma_lmerel <- function(model.lme, to.model.gene, gene, kin.subset, use.weights

if(metrics){
#Calculate R-squared
if(use.weights){
null <- stats::glm(formula = expression ~ 1, data = to.model.gene,
weights = to.model.gene$gene_weight)
if(use_weights){
null <- stats::glm(formula = expression ~ 1, data = to_model_gene,
weights = to_model_gene$gene_weight)
} else{
null <- stats::glm(formula = expression ~ 1, data = to.model.gene)
null <- stats::glm(formula = expression ~ 1, data = to_model_gene)
}

L0 <- as.vector(stats::logLik(null))
Expand Down
Loading

0 comments on commit f3381fd

Please sign in to comment.