upload some function (c9fa20e0) · Commits · Chaos / chaos.tools

R/RNA-seq.R

0 → 100755

+477 −0

Original line number	Diff line number	Diff line
		rnaseq123 <- function(
		ctl_count_file,
		obs_count_file,
		organism = "human",
		count_type = "star"
		){
		if(tolower(count_type) == "star")
		{
		x <- readDGE(
		c(ctl_count_file,obs_count_file),
		header = F,
		skip = 4
		)
		}else if(tolower(count_type) == "featurecounts")
		{
		x <- readDGE(
		c(ctl_count_file,obs_count_file),
		columns = c(1,7),
		sep = "\t",
		skip = 1
		)
		}

		group <- c(
		rep("ctl",length(ctl_count_file)),
		rep("obs",length(obs_count_file))
		)

		x$samples$group <- group

		# gtf_dt_file <- gsub("gtf.*","gtf.rdata",gtf_file)
		#
		# if(!isTRUE(str_detect(gtf_file,"rdata")))
		# {
		# if(! file.exists(gtf_dt_file)))
		# {
		# gtf_dt <- gtf_file %>%
		# import() %>%
		# as.data.table()
		#
		# save(gtf_dt,gtf_dt_file)
		# }else
		# {
		# gtf_dt <- load(gtf_dt_file)
		# }
		# }else
		# {
		# gtf_dt <- load(gtf_file)
		# }

		gtf_info <- gtf(organism,"gene") %>%
		as.data.table() %>%
		.[,.(gene_id,gene_name,type,chr = seqnames,start,end,strand)] \|>
		setkey(gene_id,gene_name)

		## it will take few minutes to load the gtf_file
		## filter the unneeded rows

		x$genes <- gtf_info[,.(gene_id,gene_name)] %>%
		unique()

		#get the correspondece between gene_id (Ensembl ID) and gene_name

		x <- x[
		filterByExpr(x, group = group),
		keep.lib.sizes = F
		] %>%
		calcNormFactors(method = "TMM")

		design <- model.matrix(~0+group)
		colnames(design) <- gsub("group", "", colnames(design))

		contr.matrix <- makeContrasts(
		contrasts = "obs-ctl",
		levels = colnames(design)
		)

		x %>%
		voomWithQualityWeights(design, plot = F) %>%
		lmFit(design) %>%
		contrasts.fit(contrasts = contr.matrix) %>%
		eBayes() %>%
		topTable(n = Inf) %>%
		as.data.table() %>%
		.[order(-logFC)] %>%
		.[,log10P := -log10(adj.P.Val)] %>%
		setnames(old = "logFC", new = "log2FC") %>%
		setkey(gene_id,gene_name) %>%
		.[gtf_info, nomatch = F] %>%
		.[,TSS_start := fifelse(strand == "+",start,end - 1)] %>%
		.[,TSS_end := TSS_start + 1]
		}

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

		if(tolower(regulation) == "up")
		{
		dt[log2FC > log2(fold_change), gene_name]
		}else if(tolower(regulation) == "down")
		{
		dt[log2FC < -log2(fold_change), gene_name]
		}else if(tolower(regulation) == "regulated")
		{
		dt[log2FC > log2(fold_change) \| log2FC < -log2(fold_change), gene_name]
		}else if(tolower(regulation) == "not_sig")
		{
		c(
		.data[adj.P.Val >= p_value, gene_name],
		dt[log2FC <= log2(fold_change) & log2FC >= -log2(fold_change), gene_name]
		)
		}else if(tolower(regulation) == "all")
		{
		.data[, gene_name]
		}
		}

		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,
		fold_change = 1.5,
		p_value = 0.05,
		special.gene = NULL,
		show_p = TRUE,
		show_fc = TRUE,
		show_gene_number = TRUE,
		sub_gene = NULL,
		sub_id = NULL,
		color_up = "#E64B35",
		color_down = "#3C5488",
		color_ns = "#BBBBBB",
		color_special = "#00a087"
		){
		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
		)[
		,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",
		default = "not_sig"
		)
		]

		stat <- merge(
		dt[,.(main = .N),regulation],
		dt[,.(sub = .N),sub_regulation] \|>
		setnames("sub_regulation","regulation")
		)[
		,label := fifelse(
		main == sub,
		paste0("N = ",main),
		paste0("N = ",sub,"/",main)
		)
		]

		top.genes <- sapply(
		c("up","down"),
		\(x)
		{
		dt[
		sub_regulation == x
		][
		order(-abs(log2FC))
		][
		0:top_gene_number,
		gene_name
		]
		}
		) \|>
		as.data.table()

		dt[
		,color := sub_regulation
		][
		gene_name %in% top.genes$up,
		label_up := gene_name
		][
		gene_name %in% top.genes$down,
		label_down := gene_name
		][
		gene_name %in% special.gene,
		label_special := gene_name
		]

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

		dt[
		label_special %in% special.gene,
		color := "special"
		]
		}

		p0 <- ggplot() +
		geom_point(
		data = dt[color == "not_sig"],
		aes(log2FC,log10P),
		color = color_ns,
		size = 0.8,
		show.legend = F
		) +
		geom_point(
		data = dt[color == "up"],
		aes(log2FC,log10P),
		color = color_up,
		size = 0.8,
		show.legend = F
		) +
		geom_point(
		data = dt[color == "down"],
		aes(log2FC,log10P),
		color = color_down,
		size = 0.8,
		show.legend = F
		) +
		geom_point(
		data = dt[color == "special"],
		aes(log2FC,log10P),
		color = color_special,
		size = 0.8,
		show.legend = F
		) +
		geom_hline(
		yintercept = -log10(p_value),
		linetype = "dashed"
		) +
		geom_vline(
		xintercept = c(-log2(fold_change),log2(fold_change)),
		linetype = "dashed"
		) +
		theme_base() +
		geom_text_repel(
		data = dt,
		aes(x = log2FC, y = log10P),
		label = dt$label_up,
		max.overlaps = Inf,
		box.padding = 1,
		na.rm = T,
		color = color_up
		) +
		geom_text_repel(
		data = dt,
		aes(x = log2FC, y = log10P),
		label = dt$label_down,
		max.overlaps = Inf,
		box.padding = 1,
		na.rm = T,
		color = color_down
		) +
		geom_text_repel(
		data = dt,
		aes(x = log2FC, y = log10P),
		label = dt$label_special,
		max.overlaps = Inf,
		box.padding = 1,
		na.rm = T,
		color = color_special
		) +
		labs(
		x = expression(log[2]("Fold Change")),
		y = expression(-log[10](adj.p-Value))
		) +
		guides(color = "none") +
		theme(plot.margin = unit(rep(1,4),'lines'))

		if(isTRUE(show_p))
		{
		p1 <- p0 +
		annotate(
		geom = "text",
		x = min(dt$log2FC),
		y = -log10(p_value),
		hjust = 0,
		vjust = 1,
		label = paste0("p-value = ",p_value)
		)
		}else
		{
		p1 <- p0
		}

		if(isTRUE(show_fc))
		{
		p2 <- p1 +
		annotate(
		geom = "text",
		x = log2(fold_change) + 0.1,
		y = max(dt$log10P) + 1,
		hjust = 0,
		label = paste0("Fold Change threshold = ",fold_change)
		)
		}else
		{
		p2 <- p1
		}

		if(isTRUE(show_gene_number))
		{
		p2 +
		annotate(
		"text",
		x = c(-log2(fold_change),min(dt$log2FC)) \|> mean(),
		y = max(dt$log10P),
		label = stat[regulation == "down",label]
		) +
		annotate(
		"text",
		x = c(log2(fold_change),max(dt$log2FC)) \|> mean(),
		y = max(dt$log10P),
		label = stat[regulation == "up",label]
		)
		}else
		{
		p2
		}
		}


		sub_volcano_plot <- function(
		id = NULL,
		name = NULL,
		fold_change = NA,
		p_value = NA,
		ref_data = rna_seq_result$DE,
		plot = TRUE,
		label_position = c(-4,12,4,12)
		){
		total_stat <- sapply(
		c("up","down"),
		function(x){
		ref_data[,.N,.(regulation)][regulation == x,N]
		}
		)

		if(is.na(fold_change))
		{
		dt <- ref_data
		}else
		{
		dt <- ref_data[log2FC > log2(fold_change) \| log2FC < -log2(fold_change)]
		}

		if(!is.na(p_value))
		{
		dt <- dt[adj.P.Val < p_value]
		}


		if(is.NULL(id))
		{
		id <- ref_data$gene_id
		}


		if(is.NULL(name))
		{
		name <- ref_data$gene_name
		}


		dt <- dt[gene_id %in% id & gene_name %in% name]

		dt_stat <- sapply(
		c("up","down"),
		function(x){
		dt[,.N,.(regulation)][regulation == x,N]
		}
		)

		label_up <- paste0("N = ",dt_stat[1]," / ",total_stat[1],"\n",round(dt_stat[1]/total_stat[1] * 100,2),"%")
		label_down <- paste0("N = ",dt_stat[2]," / ",total_stat[2],"\n",round(dt_stat[2]/total_stat[2] * 100,2),"%")

		if(isTRUE(plot))
		{
		volcano_plot(
		dt,
		top_gene_number = 0
		) +
		theme(plot.title = element_text(hjust = 0.5)) +
		annotate("text",x = label_position[1], y = label_position[2], label = label_down) +
		annotate("text",x = label_position[3], y = label_position[4], label = label_up)
		} else
		{
		dt
		}
		}

R/hello.R

deleted100755 → 0

+0 −18

Original line number	Diff line number	Diff line
		# Hello, world!
		#
		# This is an example function named 'hello'
		# which prints 'Hello, world!'.
		#
		# You can learn more about package authoring with RStudio at:
		#
		# http://r-pkgs.had.co.nz/
		#
		# Some useful keyboard shortcuts for package authoring:
		#
		# Install Package: 'Ctrl + Shift + B'
		# Check Package: 'Ctrl + Shift + E'
		# Test Package: 'Ctrl + Shift + T'

		hello <- function() {
		print("Hello, world!")
		}

R/separate_col.R

0 → 100755

+82 −0

Original line number	Diff line number	Diff line
		#' Fast separate column
		#' @examples
		#' data.frame(x=paste0(1:5,"_",letters[1:5],".",LETTERS[1:5])) %>% separate_col("x")
		#' data.frame(x=paste0(1:5,"_",letters[1:5],".",LETTERS[1:5])) %>% separate_col("x",sep = "_")
		#' data.frame(x=paste0(1:5,"_",letters[1:5],".",LETTERS[1:5])) %>% separate_col("x",select = 1)
		#' data.frame(x=paste0(1:5,"_",letters[1:5],".",LETTERS[1:5])) %>% separate_col("x",select = 2,into = "y")
		#' data.frame(x=paste0(1:5,"_",letters[1:5],".",LETTERS[1:5])) %>% separate_col("x",select = 3,into = "y",remove = F)
		#' data.frame(x=paste0(1:5,"_",letters[1:5],".",LETTERS[1:5])) %>% separate_col("x",select = 1,remove = F)
		#' data.frame(x=paste0(1:5,"_",letters[1:5],".",LETTERS[1:5])) %>% separate_col("x",select = c(1,3),remove = F)
		#' data.frame(x=paste0(1:5,"_",letters[1:5],".",LETTERS[1:5])) %>% separate_col("x",select = c(1,3),into = c("y","z"),remove = F)
		#' @export
		separate_col <- function(
		.data,
		column,
		into = NA,
		select = NA,
		sep = "[^[:alnum:]]+",
		remove = TRUE
		){
		split_columns <- as.data.table(.data)%>%
		.[[column]] %>%
		tstrsplit(split = sep) %>%
		setDT()

		default_name <- paste(
		"splited",
		column,
		1:ncol(split_columns),
		sep = "_"
		)

		if(length(select) == 1)
		{
		if(is.na(select))
		{
		select = 1:ncol(split_columns)
		}
		}

		selected_columns <- paste("splited",column,select,sep = "_")

		if(length(into) == 1)
		{
		if(is.na(into))
		{
		into <- default_name[select]

		##a simplified way when there is only 1 column were selected
		##and the original column was not need
		if(length(select) == 1)
		{
		if(!is.na(select))
		{
		if(isTRUE(remove))
		{
		into <- column
		}
		}
		}
		}
		}

		renamed_split_columns <- split_columns %>%
		setnames(default_name) %>%
		.[,..selected_columns] %>%
		setnames(into)

		if(isTRUE(remove))
		{
		cbind(
		.data[,(column) := NULL],
		renamed_split_columns
		)
		} else
		{
		cbind(
		.data,
		renamed_split_columns
		)
		}

		}

R/shinyapp_test.R

0 → 100755

+151 −0

Original line number	Diff line number	Diff line
		cvd_ui <- function(cvd_data)
		{
		cvd_server <- function(
		input,
		output
		){
		output$cvd_caption <- renderText(
		{
		paste0("Chromosome ", input$chr)
		})

		output$cvd_plot <- renderPlot(
		{
		plot_cvd(
		cvd_data,
		chr = input$chr,
		min = input$interaction_distance[1],
		max = input$interaction_distance[2],
		slop_position = c(input$slop_position_x,input$slop_position_y),
		slop_length = input$slop_length

		)
		})
		output$cvd_table <- renderDataTable({
		cvd_data[chr == input$chr]
		})
		}

		cvd_ui <- fluidPage(
		titlePanel("Decaycurve"),
		sidebarLayout(
		sidebarPanel(
		selectInput(
		inputId = "chr",
		label = "chr",
		unique(cvd_data$chr)
		),
		sliderInput(
		"interaction_distance",
		"interaction distance",
		value = c(1e4,1e7),
		min = 1e4,
		max = 1e8,
		animate = T
		),
		sliderInput(
		"slop_position_x",
		"slop position x",
		value = 4,
		min = 4,
		max = 7,
		step = 0.1,
		animate = T
		),
		sliderInput(
		"slop_position_y",
		"slop position y",
		value = 6.8,
		min = 0,
		max = 10,
		step = 0.1,
		animate = T
		),
		sliderInput(
		"slop_length",
		"slop length",
		value = 2,
		min = 0,
		max = 5,
		step = 0.1,
		animate = T
		)
		),
		mainPanel(
		h3(textOutput("cvd_caption")),
		plotOutput("cvd_plot"),
		h3("Data"),
		dataTableOutput("cvd_table")
		)
		)
		)
		shinyApp(cvd_ui,cvd_server)
		}

		rna_ui <- function(.data)
		{
		rna_server <- function(
		input,
		output
		){
		output$rna_caption <- renderText(
		{
		"Volcano plot"
		})

		output$rna_plot <- renderPlot(
		{
		volcano_plot(
		.data,
		top_gene_number = input$top_gene_number,
		fold_change = input$fold_change,
		p_value = input$p_value
		)
		})
		output$rna_table <- renderDataTable({
		.data
		})
		}

		rna_ui <- fluidPage(
		titlePanel("Transcriptome"),
		sidebarLayout(
		sidebarPanel(
		sliderInput(
		"top_gene_number",
		"Top gene number",
		value = 10,
		min = 0,
		max = 100,
		step = 1,
		animate = T
		),
		sliderInput(
		"fold_change",
		"fold change",
		value = 1.5,
		min = 0,
		max = 5,
		step = 0.1,
		animate = T
		),
		sliderInput(
		"p_value",
		"P value",
		value = 0.05,
		min = 0.001,
		max = 0.5,
		step = 0.001,
		animate = T
		)
		),
		mainPanel(
		h3(textOutput("rna_caption")),
		plotOutput("rna_plot"),
		h3("Data"),
		dataTableOutput("rna_table")
		)
		)
		)
		shinyApp(rna_ui,rna_server)
		}

R/shuffle.R

0 → 100755

+68 −0

Original line number	Diff line number	Diff line
		shuffle <- function(
		.data,
		ref = "hg19"
		){
		dt <- as.data.table(.data)


		chr_info <- chr_size(
		ref = ref,
		extra = T,
		mit = T
		) %>%
		setkey()

		first_3_names <- colnames(dt[,1:3])

		dt_first_3 <- dt[,1:3] %>%
		setnames(c("chr","start","end")) %>%
		setkey()

		if(ncol(dt) >=4)
		{
		dt_rest <- dt[,4:ncol(dt)]
		}else
		{
		dt_rest <- NULL
		}

		chr_character <- dt_first_3$chr %>%
		grepl(pattern = "chr") %>%
		unique()

		if(length(chr_character) == 1)
		{
		if(isFALSE(chr_character))
		{
		chr_info[
		,chr := gsub(chr,pattern = "chr",replacement = "")
		] %>%
		setkey()
		}
		}else
		{
		dt_first_3[
		,chr := gsub(chr,pattern = "chr",replacement = "")
		][
		,chr := paste0("chr",chr)
		] %>%
		setkey()
		}

		merge(
		dt_first_3,
		chr_info
		)[
		,peak_length := end - start
		][
		,range := length - peak_length
		][
		,new_start := sapply(range,function(x){sample(1:x,1)})
		][
		,new_end := new_start + peak_length
		][
		,.(chr,new_start,new_end)
		] %>%
		setnames(first_3_names) %>%
		cbind(dt_rest)
		}
		No newline at end of file

Admin message