Commit f7c29c93 authored by Chaos's avatar Chaos
Browse files

modified gtf.R added dturtle.R

parent 0c4e14c9
Loading
Loading
Loading
Loading

R/dturtle.R

0 → 100644
+100 −0
Original line number Diff line number Diff line
tx2gene <- function(organism){
	tx2gene_file <- file.path(
						chaos_info("data_dir"),
						paste0("tx2gene_for_dturtle",organism,"rds",sep = ".")
					)

	if(!file.exists(tx2_gene_file)){
		tx2gene <-	gtf("mouse") |>
					as.data.table() |> 
					move_columns_to_front(
						columns = c("transcript_name", "gene_name")
					)

		tx2gene$gene_name <- one_to_one_mapping(
								name = tx2gene$gene_name, 
								id = tx2gene$gene_id
							) 
		tx2gene$transcript_name <- one_to_one_mapping(
										name = tx2gene$transcript_name, 
										id = tx2gene$transcript_id
									)
		saveRDS(tx2gene,tx2gene_file,compress = "xz")
	}else{
		tx2gene <- readRDS(tx2gene_file)
	}
	tx2gene
}

dtu <-	function(
			ctl_count_file,
			ctl_name,
			obs_count_file,
			obs_name,
			result_dir,
			organism = "mouse"
){
	tx2gene <- tx2gene("mouse")

	names(ctl_count_file) <-	gsub("/quant.sf","",ctl_count_file) |>
								basename()

	names(obs_count_files) <-	gsub("/quant.sf","",ctl_count_file) |>
								basename()

	cts <-	c(ctl_count_files,obs_count_files) |> 
			import_counts(
				type = "salmon",
				tx2gene = tx2gene[,c("transcript_id", "gene_name")]
			)

	rownames(cts) <- tx2gene$transcript_name[match(rownames(cts), tx2gene$transcript_id)]

	pd <-	data.frame(
				id = colnames(cts),
				group = sub_dt[,samples],
				stringsAsFactors = FALSE
			)

	dturtle <-	cts |>
				run_drimseq(
					tx2gene = tx2gene, 
					pd = pd, 
					id_col = "id",
					cond_col = "group",
					cond_levels = c(obs_name,ctl_name),
					filtering_strategy = "bulk", 
					BPPARAM = biocpar
				) |> 
				posthoc_and_stager(
					ofdr = 0.05,
					posthoc = 0.1
				)

	dturtle$dge_analysis <- files |> 
							import_dge_counts(
								type = "salmon",
								tx2gene = tx2gene[,c("transcript_id", "gene_name")]
							) |> 
							run_deseq2(
								pd = pd,
								id_col = "id",
								cond_col = "group",
								cond_levels = c(obs_name,ctl_name),
								lfc_threshold = 0.5,
								sig_threshold = 0.01,
								dge_calling_strategy = "bulk",
								BPPARAM = biocpar
							)

	dtu_dge_genes <-	intersect(
							dturtle$sig_gene,
							dturtle$dge_analysis$results_sig$gene
						)
		
	dturtle <-	create_dtu_table(
					dturtle = dturtle,
					add_gene_metadata = list("chromosome"="seqnames"), 
					add_tx_metadata = list("tx_expr_in_max" = c("exp_in", max))
				)
}
 No newline at end of file
+0 −24
Original line number Diff line number Diff line
@@ -91,27 +91,3 @@ txdb <- function(
		txdb
	}
}

tx2gene <- function(organism){
	tx2gene_file <- file.path(dir,paste0("tx2gene_for_dturtle",organism,"rds",sep = "."))
	if(!file.exists(tx2_gene_file)){
		tx2gene <-	gtf("mouse") |>
					as.data.table() |> 
					move_columns_to_front(
						columns = c("transcript_name", "gene_name")
					)

		tx2gene$gene_name <- one_to_one_mapping(
								name = tx2gene$gene_name, 
								id = tx2gene$gene_id
							) 
		tx2gene$transcript_name <- one_to_one_mapping(
										name = tx2gene$transcript_name, 
										id = tx2gene$transcript_id
									)
		saveRDS(tx2gene,tx2gene_file,compress = "xz")
	}else{
		tx2gene <- readRDS(tx2gene_file)
	}
	tx2gene
}
 No newline at end of file