volcano plot now could plot with hyperbola curve (221b1504) · Commits · Chaos / chaos.tools

R/RNA-seq.R

+174 −120

Original line number	Diff line number	Diff line
		@@ -73,30 +73,68 @@ rnaseq123 <- function(
		.[,TSS_end := TSS_start + 1]
		}

		de.dt <- function(
		.data,
		fold_change = 2,
		p_value = 0.05,
		hyperbola = FALSE
		){
		if(isFALSE(hyperbola)){
		as.data.table(
		.data
		)[
		,regulation := fcase(
		adj.P.Val < p_value & log2FC > log2(fold_change), "up",
		adj.P.Val < p_value & log2FC < -log2(fold_change), "down",
		default = "not_sig"
		)
		][]
		}else{
		as.data.table(
		.data
		)[
		,hyperbola_log10P := fifelse(
		log2FC <0,
		-1/(log2FC + log2(fold_change)) - log10(p_value),
		1/(log2FC - log2(fold_change)) - log10(p_value)
		)
		][
		,regulation := fcase(
		log2FC > log2(fold_change) & log10P >= hyperbola_log10P,"up",
		log2FC < -log2(fold_change) & log10P >= hyperbola_log10P,"down",
		default = "not_sig"
		)
		][]
		}
		}

		de.gene <- function(
		.data,
		regulation,
		fold_change = 2,
		p_value = 0.05,
		keytype = "SYMBOL"
		keytype = "SYMBOL",
		hyperbola = FALSE
		){
		dt <-.data[adj.P.Val < p_value]

		if(tolower(regulation) == "up"){
		sub_dt <- dt[log2FC > log2(fold_change)]
		}else if(tolower(regulation) == "down"){
		sub_dt <- dt[log2FC < -log2(fold_change)]
		}else if(tolower(regulation) == "regulated"){
		sub_dt <- dt[log2FC > log2(fold_change) \| log2FC < -log2(fold_change)]
		}else if(tolower(regulation) == "not_sig"){
		sub_dt <- rbind(
		.data[adj.P.Val >= p_value],
		dt[log2FC <= log2(fold_change) & log2FC >= -log2(fold_change)]
		) \|>
		unique()
		}else if(tolower(regulation) == "all"){
		sub_dt <- .data
		re_str <- regulation %>%
		tolower() %>%
		gsub("-","_",.)

		dt <- de.dt(
		.data,
		fold_change,
		p_value,
		hyperbla
		)

		if(re_str %in% c("up","down","not_sig"))
		sub_dt <- dt[regulation == tolower(regulation)]
		else if(re_str == "regulated"){
		sub_dt <- dt[regulation %in% c("up","down")]
		}else{
		sub_dt <- dt
		}

		if(keytype == "SYMBOL"){
		sub_dt[,gene_name]
		}else if(keytype == "ENSEMBL"){
		@@ -105,6 +143,7 @@ de.gene <- function(
		separate_col(
		sub_dt,
		"gene_id",
		sep = "[.]",
		select = 1,
		into = "short_ENSEMBL",
		remove = F
		@@ -112,79 +151,6 @@ de.gene <- function(
		}
		}

		simple_ego <- function(
		.data,
		regulation,
		organism = "human",
		ont = "ALL",
		fold_change = 2,
		p_value = 0.05
		){
		gene <- de.gene(.data, regulation, fold_change, p_value)

		if(organism == "human")
		{
		library(org.Hs.eg.db)

		clusterProfiler::enrichGO(
		gene = gene,
		keyType = "SYMBOL",
		ont = ont,
		OrgDb = org.Hs.eg.db
		)
		}else if(organism == "mouse")
		{
		library(org.Mm.eg.db)

		clusterProfiler::enrichGO(
		gene = gene,
		keyType = "SYMBOL",
		ont = ont,
		OrgDb = org.Mm.eg.db
		)
		}
		}

		simple_kegg <- function(
		.data,
		regulation,
		organism = "human",
		fold_change = 2,
		p_value = 0.05
		){
		library(clusterProfiler)
		gene <- de.gene(.data, regulation, fold_change, p_value)

		if(organism == "human")
		{
		library(org.Hs.eg.db)
		organism_code_name <- "hsa"
		entrez_id <- bitr(
		geneID = gene,
		fromType = "SYMBOL",
		toType = "ENTREZID",
		OrgDb = org.Hs.eg.db
		)
		}else if(organism == "mouse")
		{
		library(org.Mm.eg.db)
		organism_code_name <- "mmu"
		entrez_id <- bitr(
		geneID = gene,
		fromType = "SYMBOL",
		toType = "ENTREZID",
		OrgDb = org.Mm.eg.db
		)
		}

		enrichKEGG(
		gene = entrez_id,
		keyType = "kegg",
		organism = organism_code_name
		)
		}


		volcano_plot <- function(
		.data,
		top_gene_number = 10,
		@@ -199,22 +165,23 @@ volcano_plot <- function(
		color_up = "#E64B35",
		color_down = "#3C5488",
		color_ns = "#BBBBBB",
		color_special = "#00a087"
		color_special = "#00a087",
		hyperbola = FALSE,
		p_line_type = "dashed",
		fc_line_type = "dashed",
		p_line_color = "black",
		left_fc_line_color = "black",
		right_fc_line_color = "black"
		){
		pacman::p_load(data.table,ggplot2,ggthemes,ggrepel)

		if(is.null(sub_gene)){sub_gene <- .data$gene_name}
		if(is.null(sub_id)){sub_id <- .data$gene_id}

		dt <- as.data.table(
		.data
		dt <- de.dt(
		.data,
		fold_change,
		p_value,
		hyperbola
		)[
		,regulation := fcase(
		adj.P.Val < p_value & log2FC > log2(fold_change), "up",
		adj.P.Val < p_value & log2FC < -log2(fold_change), "down",
		default = "not_sig"
		)
		][
		,sub_regulation := fcase(
		regulation == "up" & gene_name %in% sub_gene & gene_id %in% sub_id, "up",
		regulation == "down" & gene_name %in% sub_gene & gene_id %in% sub_id, "down",
		@@ -263,8 +230,7 @@ volcano_plot <- function(
		label_special := gene_name
		]

		if(!is.null(special.gene))
		{
		if(!is.null(special.gene)){
		special.gene <- special.gene %>% .[. %in% dt$gene_name]

		dt[
		@@ -273,7 +239,7 @@ volcano_plot <- function(
		]
		}

		p0 <- ggplot() +
		plot <- ggplot() +
		geom_point(
		data = dt[color == "not_sig"],
		aes(log2FC,log10P),
		@@ -301,15 +267,38 @@ volcano_plot <- function(
		color = color_special,
		size = 0.8,
		show.legend = F
		) +
		)

		if(isFALSE(hyperbola)){
		plot <- plot +
		geom_hline(
		yintercept = -log10(p_value),
		linetype = "dashed"
		linetype = p_line_type,
		color = p_line_color
		) +
		geom_vline(
		xintercept = c(-log2(fold_change),log2(fold_change)),
		linetype = "dashed"
		linetype = fc_line_type,
		color = right_fc_line_color
		)
		}else{
		plot <- plot +
		geom_line(
		data = dt[log2FC < 0],
		aes(log2FC,hyperbola_log10P),
		linetype = fc_line_type,
		color = left_fc_line_color
		) +
		geom_line(
		data = dt[log2FC > 0],
		aes(log2FC,hyperbola_log10P),
		linetype = fc_line_type,
		color = right_fc_line_color
		) +
		scale_y_continuous(limits = c(0,max(dt$log10P) + 1))
		}

		plot <- plot +
		theme_base() +
		geom_text_repel(
		data = dt,
		@@ -347,7 +336,7 @@ volcano_plot <- function(

		if(isTRUE(show_p))
		{
		p1 <- p0 +
		plot <- plot +
		annotate(
		geom = "text",
		x = min(dt$log2FC),
		@@ -356,14 +345,10 @@ volcano_plot <- function(
		vjust = 1,
		label = paste0("p-value = ",p_value)
		)
		}else
		{
		p1 <- p0
		}

		if(isTRUE(show_fc))
		{
		p2 <- p1 +
		if(isTRUE(show_fc)){
		plot <- plot +
		annotate(
		geom = "text",
		x = log2(fold_change) + 0.1,
		@@ -371,14 +356,11 @@ volcano_plot <- function(
		hjust = 0,
		label = paste0("Fold Change threshold = ",fold_change)
		)
		}else
		{
		p2 <- p1
		}

		if(isTRUE(show_gene_number))
		{
		p2 +
		plot <- plot +
		annotate(
		"text",
		x = c(-log2(fold_change),min(dt$log2FC)) \|> mean(),
		@@ -391,13 +373,85 @@ volcano_plot <- function(
		y = max(dt$log10P),
		label = stat[regulation == "up",label]
		)
		}else
		}

		plot
		}

		simple_ego <- function(
		.data,
		regulation,
		organism = "human",
		ont = "ALL",
		fold_change = 2,
		p_value = 0.05
		){
		gene <- de.gene(.data, regulation, fold_change, p_value)

		if(organism == "human")
		{
		library(org.Hs.eg.db)

		clusterProfiler::enrichGO(
		gene = gene,
		keyType = "SYMBOL",
		ont = ont,
		OrgDb = org.Hs.eg.db
		)
		}else if(organism == "mouse")
		{
		library(org.Mm.eg.db)

		clusterProfiler::enrichGO(
		gene = gene,
		keyType = "SYMBOL",
		ont = ont,
		OrgDb = org.Mm.eg.db
		)
		}
		}

		simple_kegg <- function(
		.data,
		regulation,
		organism = "human",
		fold_change = 2,
		p_value = 0.05
		){
		library(clusterProfiler)
		gene <- de.gene(.data, regulation, fold_change, p_value)

		if(organism == "human")
		{
		p2
		library(org.Hs.eg.db)
		organism_code_name <- "hsa"
		entrez_id <- bitr(
		geneID = gene,
		fromType = "SYMBOL",
		toType = "ENTREZID",
		OrgDb = org.Hs.eg.db
		)
		}else if(organism == "mouse")
		{
		library(org.Mm.eg.db)
		organism_code_name <- "mmu"
		entrez_id <- bitr(
		geneID = gene,
		fromType = "SYMBOL",
		toType = "ENTREZID",
		OrgDb = org.Mm.eg.db
		)
		}

		enrichKEGG(
		gene = entrez_id,
		keyType = "kegg",
		organism = organism_code_name
		)
		}



		sub_volcano_plot <- function(
		id = NULL,
		name = NULL,

Admin message