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

R/Decay_curve.R

0 → 100644

+162 −0

Original line number	Diff line number	Diff line
		cvd <- function(
		hic_file = "/mnt/d/work/Hi-C/hic_files_hg19/DMSO_DHT.hic",
		resolution = 1e4,
		norm = "KR",
		chr = NA
		){
		BPPARAM <- hic_file \|> length() \|> MulticoreParam()
		bpstart(BPPARAM)

		c_count <- hic_file \|>
		bplapply(
		function(x){
		hic_interaction(
		hic_file = x[1],
		chr_list = chr,
		resolution = resolution,
		norm = norm,
		inter = "intra"
		)[
		,distance := chr2_bin - chr1_bin
		][
		,sum(counts),.(chr1,distance)
		][
		,.(sample = base_name(x[1]),chr = chr1,distance,counts = V1)
		]
		},
		BPPARAM = BPPARAM
		) \|>
		rbindlist()

		bpstop(BPPARAM)

		s_count <- c_count[
		,sum(counts),
		.(sample,distance)
		][
		,.(sample,chr = "all",distance,counts = V1)
		]

		rbind(c_count,s_count)
		}

		plot_cvd <- function(
		.data,
		sample = NA,
		chr = "all",
		slop_ratio = -1,
		slop_length = 2,
		slop_position = NA,
		min = 1e4,
		max = 1e7
		){
		s <- sample

		if(length(s) == 1)
		{
		if(is.na(s))
		{
		s <- unique(.data$sample)
		}
		}

		c <- chr

		# x_breaks <- data.table(
		# breaks = 4:log(max),
		# labels = c("10Kb","100Kb","1Mb","10Mb")
		# )
		# #as the default resolution is 10k, the closest distance is 10k except 0
		# #the contact that longer than 100Mb is rare, and we don't really need the plot the contact longer than 10Mb

		x_breaks <- log10(min):log10(max) \|>
		lapply(
		function(i)
		{
		if(i < 6)
		{
		data.table(breaks = i, labels = paste0(10^(i-3),"kb"))
		} else
		{
		data.table(breaks = i, labels = paste0(10^(i-6),"Mb"))
		}
		}
		) \|>
		rbindlist()

		dt <- .data[
		sample %in% s &
		chr %in% c &
		distance %between% c(min,max)
		]

		if(length(slop_position) == 1)
		{
		if(is.na(slop_position))
		{
		slop_position[1] <- log10(min)
		}

		slop_position[2] <- max(log10(dt$counts))
		}

		slop_position[3] <- slop_position[1] + slop_length/sqrt(1 + slop_ratio^2)
		slop_position[4] <- slop_position[2] - slop_ratio * (slop_position[1] - slop_position[3])

		p <- ggplot(
		dt,
		aes(
		log10(distance),
		log10(counts),
		color = sample
		)
		) +
		geom_line(alpha = 0.8) +
		theme_prism() +
		theme(legend.position = c(0.75,0.8)) +
		scale_x_continuous(
		breaks = x_breaks$breaks,
		minor_breaks = sapply(log10(min):(log10(max) - 1),function(x){log10(1:9) + x}),
		#1:9 代表大格之间的9个小格，x_breaks$breaks代表所有的大刻度
		labels = x_breaks$labels,
		guide = "prism_offset_minor"
		) +
		scale_y_continuous(
		limits = c(
		floor(log10(min(dt$counts))),
		ceiling(log10(max(dt$counts)))
		),
		guide = "prism_offset"
		) +
		labs(
		x = "Contact Distance",
		y = "Contact frequency"
		) +
		scale_color_npg() +
		annotate(
		geom = "text",
		x = mean(slop_position[c(1,3)]),
		y = mean(slop_position[c(2,4)]),
		parse = T,
		label = paste0("italic(s)^",slop_ratio)
		) +
		annotate(
		geom = "segment",
		x = slop_position[1],
		y = slop_position[2],
		xend = slop_position[3],
		yend = slop_position[4],
		linetype = 2,
		color = "#3c5488",
		alpha = 0.8
		)

		if(length(chr) == 1)
		{
		p
		}else
		{
		p + facet_grid(~chr)
		}
		}

R/anchor_overlap.R

0 → 100644

+154 −0

Original line number	Diff line number	Diff line
		reform <- function(
		.data,
		c = "V1",
		s = "V2",
		e = "V3",
		p = "V4",
		ref = FALSE
		){
		if(isFALSE(ref))
		{
		p_name <- "peak"
		}else
		{
		p_name <- "anchor"
		}

		.data %>%
		as.data.table() %>%
		setnames(
		c(c,s,e,p),
		c("chr","start","end",p_name)
		) %>%
		.[,`:=`("chr",chr_omit(chr))] %>%
		setkey(chr,start,end) %>%
		.[]
		}

		anchor_overlap <- function(
		d0 = anchor,
		d1 = tss,
		d2 = tre,
		d3 = cis_element,
		d4 = chip,
		d5 = SE
		){
		foverlaps(
		d0,
		rbind(d1,d2,d3,d4,d5) %>%
		setkey(chr,start,end),
		nomatch = NULL
		)[
		,.N,
		.(loop,loop_type,anchor,peak_type,description)
		][
		,description2 := fcase(
		description == "High-CTCF", "CTCF",
		description == "High-H3K27ac", "AC",
		description == "High-H3K27ac,High-CTCF", "both",
		description == "not_sig", "ns"
		)
		][
		is.na(description2),
		description2 := description
		][
		,peak_type := gsub("cis-element","CE",peak_type)
		][
		,.(loop,anchor,peak_type,description2,N)
		]
		}

		anchor_overlap_stat <- function(
		.data
		){

		sub <- function(
		.data = t,
		l1 = "TRE",
		l1_status = "acquired",
		l2 = "TSS",
		l2_status = "up"
		){
		left_l1 <- paste("left",l1,l1_status,sep = "_")
		right_l1 <- paste("right",l1,l1_status,sep = "_")
		left_l2 <- paste("left",l2,l2_status,sep = "_")
		right_l2 <- paste("right",l2,l2_status,sep = "_")

		.data[
		peak_type %in% c(l1,l2)
		] %>%
		dcast(
		loop ~ anchor + peak_type + description2,
		value.var = "N",
		fill = 0
		) %>%
		setnames(
		c(left_l1,right_l1,left_l2,right_l2),
		c("left_l1","right_l1","left_l2","right_l2")
		) %>%
		.[
		(left_l1 > 0 & right_l2 > 0) \|
		(left_l2 > 0 & right_l1 > 0),
		.(loop,left_l1,right_l1,left_l2,right_l2)
		] %>%
		.[
		,status := paste(l1,l1_status,l2,l2_status,sep = "_")
		]
		}

		.data[,.(peak_type,description2)] %>%
		unique() %>%
		.[,l1 := paste(peak_type,description2,sep = "_")] %>%
		.[,l1] %>%
		dt_comb(rep = 2,inter = "inter")%>%
		separate_col("V1",into = c("l1","l1_status")) %>%
		separate_col("V2",into = c("l2","l2_status")) %>%
		.[
		(l1 != l2) &
		!((l1 == "AR" & l2 == "ERG") \| (l1 == "ERG" & l2 == "AR"))
		] %>%
		apply(
		1,
		function(x)
		{
		sub(
		.data = .data,
		l1 = x[1],
		l1_status = x[2],
		l2 = x[3],
		l2_status = x[4]
		)
		}
		) %>%
		rbindlist() %>%
		separate_col(
		"loop",
		select = 1,
		into = "loop_type",
		remove = F
		)
		}

		anchor_overlap_stat_plot <- function(
		.data,
		order = "acquired",
		remove = "stable"
		){
		all_loop_type <- c("acquired","strengthened","stable","weakened","lost")

		.data[,.N,.(status,loop_type)] %>%
		complete_dt(
		c("status","loop_type"),
		fill = 0
		) %>%
		.[loop_type != remove] %>%
		.[!status %in% status[str_detect(status,"AR_0h")]] %>%
		.[,p := N/sum(N),.(status)] %>%
		.[,status := factor(status,.[loop_type == order][order(p),status])] %>%
		.[,loop_type := factor(loop_type,all_loop_type[!all_loop_type %in% remove])] %>%
		ggplot(aes(status,N,fill = loop_type)) +
		geom_bar(stat = "identity",position = "fill") +
		theme_bw() +
		theme(axis.text.x = element_text(angle = 45,hjust = 1)) +
		scale_fill_npg()
		}

R/apa_plot.R

0 → 100644

+157 −0

Original line number	Diff line number	Diff line

		apa_plot <- function(
		.data,
		corner_size = 6,
		min = NA,
		min_color = "#3c5488",
		max = NA,
		max_color = "#e64b35",
		smooth = FALSE,
		border_color = "#000000",
		number_size = 5,
		number_color = "#000000",
		digit = 3,
		legend_title = "value"
		){
		cs <- corner_size

		shift <- 0.53

		dt <- data.table(.data)

		l <- nrow(dt)

		a <- l - cs + 1

		central_pixel <- paste0(
		"V",
		ceiling(l/2),
		ceiling(l/2)
		)

		all <- unlist(dt)

		p2m <- (all[central_pixel]/mean(all[names(all) != central_pixel])) %>% round(digit)
		p2ul <- (all[central_pixel]/mean(unlist(dt[1:cs,1:cs]))) %>% round(digit)
		p2ur <- (all[central_pixel]/mean(unlist(dt[1:cs,a:l]))) %>% round(digit)
		p2ll <- (all[central_pixel]/mean(unlist(dt[a:l,1:cs]))) %>% round(digit)
		p2lr <- (all[central_pixel]/mean(unlist(dt[a:l,a:l]))) %>% round(digit)

		if(is.na(min)){min <- min(all)}
		if(is.na(max)){max <- max(all)}

		p_base <- dt[,rn := factor(1:.N,levels = .N:1)] %>%
		melt("rn") %>%
		ggplot(aes(variable, rn, fill = value))

		if(smooth)
		{
		p_heatmap <- p_base + geom_raster(interpolate = T)
		}else
		{
		p_heatmap <- p_base + geom_tile()
		}

		p_heatmap +
		annotate(
		"rect",
		xmin= c(0,0,l - cs, l - cs) + shift,
		xmax = c(cs +1,cs + 1, l + 1, l + 1) - shift,
		ymin = c(0, l - cs, 0, l - cs) + shift,
		ymax =c(cs + 1, l + 1, cs + 1, l + 1) - shift,
		color = border_color,
		alpha = 0
		) +
		annotate(
		"text",
		x= (cs/2) + shift,
		y = (cs/2) + shift,
		label = p2ll,
		size = number_size,
		color = number_color
		) +
		annotate(
		"text",
		x= (l - cs/2) + shift,
		y = (cs/2) + shift,
		label = p2lr,
		size = number_size,
		color = number_color
		) +
		annotate(
		"text",
		x= (cs/2) + shift,
		y = (l - cs/2) + shift,
		label = p2ul,
		size = number_size,
		color = number_color
		) +
		annotate(
		"text",
		x= (l - cs/2) + shift,
		y = (l - cs/2) + shift,
		label = p2ur,
		size = number_size,
		color = number_color
		) +
		theme_void() +
		scale_fill_gradient2(
		limits = c(min,max),
		low = min_color,
		mid = "white",
		high = max_color,
		midpoint = mean(c(min,max))
		) +
		labs(fill = legend_title)
		}


		apa_compare_plot <- function(
		dt1,
		dt2,
		corner_size = 6,
		min = NA,
		min_color = "#3c5488",
		max = NA,
		max_color = "#e64b35",
		smooth = FALSE,
		border_color = "#000000",
		number_size = 5,
		number_color = "#000000",
		digit = 3,
		legend_title = "value"
		){

		if(is.na(min)){min <- floor(min(dt1,dt2))}
		if(is.na(max)){max <- ceiling(max(dt1,dt2))}

		apa_plot(
		data.table(dt1),
		corner_size = corner_size,
		min = min,
		max = max,
		min_color = min_color,
		max_color = max_color,
		smooth = smooth,
		border_color = border_color,
		number_size = number_size,
		number_color = number_color,
		digit = digit,
		legend_title = legend_title
		) +
		apa_plot(
		data.table(dt2),
		corner_size = corner_size,
		min = min,
		max = max,
		min_color = min_color,
		max_color = max_color,
		smooth = smooth,
		border_color = border_color,
		number_size = number_size,
		number_color = number_color,
		digit = digit,
		legend_title = legend_title
		) +
		plot_layout(guides = "collect")
		}

R/bin.R

0 → 100644

+31 −0

Original line number	Diff line number	Diff line
		bin <- function(
		.data,
		column = NA,
		x = 0,
		positive = 1,
		negative = 0
		){
		dt <- data.table(.data)

		if(length(column) == 1)
		{
		if(is.na(column)){
		BPPARAM <- colnames(.data) \|>
		length() \|>
		MulticoreParam()
		bpstart(BPPARAM)
		column <- bplapply(
		colnames(.data),
		function(i){if(is.numeric(dt[[i]])){i}},
		BPPARAM = BPPARAM
		) \|>
		unlist()
		bpstop(BPPARAM)
		}
		}

		numeric_dt <- dt[, ..column]
		numeric_dt[numeric_dt > x] <- positive
		numeric_dt[numeric_dt <= x] <- negative
		dt[,(column) := numeric_dt][]
		}

R/chaos_color.R

0 → 100644

+31 −0

Original line number	Diff line number	Diff line
		chaos_color <- function(
		){
		c(
		"#E64B35", ## "Cinnabar"
		"#4DBBD5", ## "Shakespeare"
		"#00A087", ## "PersianGreen"
		"#3C5488", ## "Chambray"
		"#F39B7F", ## "Apricot"
		"#8491B4", ## "WildBlueYonder"
		"#91D1C2", ## "MonteCarlo"
		"#DC0000", ## "Monza"
		"#7E6148", ## "RomanCoffee"
		"#B09C85" ## "Sandrift"
		)
		}

		scale_fill_chao <- function(
		n
		){
		scale_fill_manual(
		values = colorRampPalette(chaos_color())(n)
		)
		}

		scale_color_chao <- function(
		n
		){
		scale_color_manual(
		values = colorRampPalette(chaos_color())(n)
		)
		}
		No newline at end of file

Admin message