Commit 4a88e8fe authored by Chaos's avatar Chaos
Browse files

add rnaseq123_multi

parent cbc6b45e
Loading
Loading
Loading
Loading
+82 −1
Original line number Diff line number Diff line
@@ -601,7 +601,88 @@ rnaseqqc <- function(
	}
}

rnaseq123_multi <- function(
				dge_list = NULL,
				count_files,
				group,
				compare,
				organism = "human",
				count_type = "star"
){
	if(is.null(dge_list)){
		if(tolower(count_type) == "star"){
			cols <- c(1, 2)
			header_str <- FALSE
			skip_rows <- 4
		}else if(tolower(count_type) == "featurecounts"){
			cols <- c(1, 7)
			header_str <- FALSE
			skip_rows <- 1
		}else if(tolower(count_type) == "salmon"){
			cols <- c(1,5)
			header_str <- TRUE
			skip_rows <- 0
		}

		x <-	readDGE(
					count_files,
					columns = cols,
					header = header_str,
					skip = skip_rows
				)
	}else if(class(dge_list) == "DGEList"){
		x <- dge_list
	}

	x$samples$group <- group

	if(count_type %in% c("star","featurecounts")){
		gtf_info <- gtf(organism,"gene") %>% 
					as.data.table() %>% 
					.[,.(gene_id,gene_name,type = gene_type,chr = seqnames,start,end,strand)] |>
					setkey(gene_id,gene_name)
	}else if(count_type == "salmon"){
		gtf_info <- gtf(organism,"transcript") %>% 
					as.data.table() %>% 
					.[,.(gene_id = transcript_id,gene_name = transcript_name,type = transcript_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 = compare,
						levels = colnames(design)
					)

	x %>%
	voomWithQualityWeights(design, plot = F) %>%
	lmFit(design) %>%
	contrasts.fit(contrasts = contr.matrix) %>%
	eBayes() %>%
	topTable(n = Inf) %>%
	as.data.table() %>%
	.[,log10P := -log10(adj.P.Val)] %>%
	merge(gtf_info[,.(gene_id,gene_name,type)])
}

## the input data for venn_plot should come from rnaseq123_multi()
venn_plot <-	function(
					.data,
					regulation = "up",
@@ -611,7 +692,7 @@ venn_plot <- function(
					list = FALSE,
					group = NULL
){
	other_cols <-  c("gene_id","gene_name","AveExpr","F","P.Value","adj.P.Val")
	other_cols <-  c("gene_id","gene_name","AveExpr","F","P.Value","adj.P.Val","log10P","type")
	samples <-	colnames(.data) %>% 
				.[! . %in% other_cols]