#' Relevant probes for gain, loss, and change in methylation status {VALIDATED}
#'
#' This function returns relevant probe lists for each node that exhibit a
#' gain or loss in methylation status. It also returns three sets of branch
#' labels: number of methylation changes, number of methylation gains, and
#' number of methylation losses.
#'
#' @param tree ape::phylo tree
#' @param probe_node_matrix matrix storing the methylation status across all
#' probes and nodes on a tree
#'
#' @return a structure containing a list of probes with differentially
#' methylated probes for each node, a list of branch labels for the number of
#' total changes in methylation, a list of branch labels for the number of
#' gains in methylation, a list of branch labels for the number of losses in
#' methylation
#' @examples
#' probe_information <- generate_branch_changes(tree, probe_node_matrix)
#' relevant_probes <- probe_information$relevant_probes
#' total_changes_label <- probe_information$total_changes_label
#' meth_gains_label <- probe_information$meth_gains_label
#' meth_losses_label <- probe_information$meth_losses_label
#' ... some visualization
#' @export
generate_branch_changes <- function(tree, probe_node_matrix) {
branch_changes = data.frame(total_changes=integer(), meth_gains=integer(), meth_losses=integer())
relevant_probes = c()
for (i in 1:nnodes(tree)) {
row = as.matrix(tree$edge[i,])
i_ancestor = toString(row[1])
i_descendant = toString(row[2])
probe_node_matrixcol = rownames(probe_node_matrix)
branch_changes[i, 'total_changes'] = length(probe_node_matrixcol[probe_node_matrix[[i_ancestor]] != probe_node_matrix[[i_descendant]]])
meth_gains_probes = probe_node_matrixcol[probe_node_matrix[[i_ancestor]] != probe_node_matrix[[i_descendant]] & probe_node_matrix[[i_descendant]] == 1 & probe_node_matrix[[i_ancestor]] == 0]
branch_changes[i, 'meth_gains'] = length(meth_gains_probes)
meth_losses_probes = probe_node_matrixcol[probe_node_matrix[[i_ancestor]] != probe_node_matrix[[i_descendant]] & probe_node_matrix[[i_descendant]] == 0 & probe_node_matrix[[i_ancestor]] == 1]
branch_changes[i, 'meth_losses'] = length(meth_losses_probes)
relevant_probes[[paste(toString(i), "_gain", sep = "")]] = meth_gains_probes
relevant_probes[[paste(toString(i), "_loss", sep = "")]] = meth_losses_probes
}
# Define the branch labels based on the total number of changes
total_changes_label = as.numeric(branch_changes$total_changes)
# Define branch labels based on the number of methyl gains
meth_gains_label = as.numeric(branch_changes$meth_gains)
# Define branch labels based on the number of methyl losses
meth_losses_label = as.numeric(branch_changes$meth_losses)
return(list('relevant_probes' = relevant_probes, 'total_changes_label' = total_changes_label, 'meth_gains_label' = meth_gains_label, 'meth_losses_label' = meth_losses_label))
}
#' Relevant genes for gain and loss of methylation {VALIDATED}
#'
#' This function returns relevant gene lists for each node that exhibit a
#' gain or loss in methylation status.
#'
#' @param tree ape::phylo tree
#' @param relevant_probes list of differentially methylated probes for each
#' node
#' @param reference reference file
#'
#' @return list of differentially methylated genes for each node
#' @examples
#' relevant_genes <- get_relevant_genes(tree, relevant_probes, reference)
#' ... some visualization
#' @export
get_relevant_genes <- function(tree, relevant_probes, reference) {
relevant_genes = c()
edges = paste(tree$edge[,1], tree$edge[,2], sep="-")
gain_name = paste(edges, "_gain", sep="")
loss_name = paste(edges, "_loss", sep="")
edge_names = append(gain_name, loss_name)
for (i in 1:nnodes(tree)) {
i_gain = paste(toString(i), "_gain", sep = "")
i_loss = paste(toString(i), "_loss", sep = "")
# Retrieve genes from specified probes
relevant_genes[[i_gain]] = reference[rownames(reference) %in% relevant_probes[[i_gain]], "gene"]
relevant_genes[[i_loss]] = reference[rownames(reference) %in% relevant_probes[[i_loss]], "gene"]
# Remove NA genes
relevant_genes[[i_gain]] = relevant_genes[[i_gain]][!is.na(relevant_genes[[i_gain]])]
relevant_genes[[i_loss]] = relevant_genes[[i_loss]][!is.na(relevant_genes[[i_loss]])]
i_descendant = as.numeric(tree$edge[i, 1])
i_ancestor = as.numeric(tree$edge[i, 2])
names(relevant_genes)[2 * i - 1] = edge_names[grepl(i_descendant, edge_names) & grepl(i_ancestor, edge_names) & grepl("gain", edge_names)]
names(relevant_genes)[2 * i] = edge_names[grepl(i_descendant, edge_names) & grepl(i_ancestor, edge_names) & grepl("loss", edge_names)]
}
return(relevant_genes)
}
#' Probes along a given gene {VALIDATED}
#'
#' This function returns a list of probes along a particular gene
#'
#' @param reference reference file
#' @param gene gene name
#'
#' @return list of probes along a particular gene
#' @examples
#' gene_probes <- get_probes_of_gene(reference, gene)
#' ... some visualization
#' @export
get_probes_of_gene <- function(reference, gene) {
gene_probes = rownames(reference[reference$gene == gene & !is.na(reference$gene), ])
return(gene_probes)
}
#' Proportion of probes methylated, unmethylated, ambiguous at each node
#' across a tree {VALIDATED}
#'
#' This function returns the proportion of each probes that are methylated,
#' unmethylated, ambiguous at each node across a tree.
#'
#' @param tree ape::phylo tree
#' @param probe_node_matrix matrix storing the methylation status across all
#' probes and nodes on a tree
#' @param reference reference file
#' @param gene gene name
#'
#' @return a structure composed of two lists that contain methylation state
#' proportions for internal nodes and methylation state proportions for leaves
#' @examples
#' gene_information <- track_gene(tree, probe_node_matrix, reference, 'CD81')
#' thermo_prop_internal_nodes = gene_information$thermo_prop_internal_nodes
#' thermo_prop_leaves = gene_information$thermo_prop_leaves
#' ... some visualization
#' @export
track_gene <- function(tree, probe_node_matrix, reference, gene) {
gene_probes <- get_probes_of_gene(reference, gene)
probe_node_matrix_states = probe_node_matrix[rownames(probe_node_matrix) %in% gene_probes,]
probe_node_matrix_states = data.frame(t(probe_node_matrix_states))
probe_node_matrix_states[probe_node_matrix_states == "1.0"] = 1
probe_node_matrix_states[probe_node_matrix_states == "0.0"] = 0
nprobes = ncol(probe_node_matrix_states)
thermo_prop = c()
thermo_prop$"1" = rowSums(probe_node_matrix_states[,1:nprobes] == 1) / nprobes
thermo_prop$"0" = rowSums(probe_node_matrix_states[,1:nprobes] == 0) / nprobes
thermo_prop$"0.5" = rowSums(probe_node_matrix_states[,1:nprobes] == 0.5) / nprobes
thermo_prop = data.frame(thermo_prop, check.names=FALSE, row.names = NULL)
rownames(thermo_prop) = rownames(probe_node_matrix_states)
# Use leaves(tree) function
thermo_prop_internal_nodes <- thermo_prop[1:5,]
thermo_prop_leaves <- thermo_prop[6:11,]
return(list('thermo_prop_internal_nodes' = thermo_prop_internal_nodes, 'thermo_prop_leaves' = thermo_prop_leaves))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.