update (35a48526) · Commits · Chaos / chaos.tools

R/RNA-seq.R

+33 −25

Original line number	Diff line number	Diff line
		@@ -394,7 +394,7 @@ volcano_plot <- function(
		) +
		labs(
		x = expression(log[2]("Fold Change")),
		y = expression(-log[10]("BH.adj.p.Value"))
		y = expression(-log[10]("adj.p.Value"))
		) +
		theme(
		plot.margin = unit(rep(1,4),'lines'),
		@@ -666,18 +666,20 @@ rnaseq123_multi <- function(
		compare,
		organism = "human",
		gtf_dt = NULL,
		count_type = "star"
		count_type = "star",
		result = "DE"
		){
		if(is.null(dge_list)){
		if(tolower(count_type) == "star"){
		ct <- tolower(count_type)
		if(ct == "star"){
		cols <- c(1, 2)
		header_str <- FALSE
		skip_rows <- 4
		}else if(tolower(count_type) == "featurecounts"){
		}else if(ct == "featurecounts"){
		cols <- c(1, 7)
		header_str <- FALSE
		skip_rows <- 1
		}else if(tolower(count_type) == "salmon"){
		}else if(ct == "salmon"){
		cols <- c(1,5)
		header_str <- TRUE
		skip_rows <- 0
		@@ -697,58 +699,64 @@ rnaseq123_multi <- function(

		if(is.null(gtf_dt)){
		if(count_type %in% c("star","featurecounts")){
		gtf_info <- gtf(organism,"gene") %>%
		gtf_dt <- 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") %>%
		gtf_dt <- 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)
		}
		}else{
		gtf_info <- gtf_dt %>%
		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()
		x$genes <- gtf_dt

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

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

		re <- toupper(result)
		if(re == "DE"){
		design <- model.matrix(~0+group)
		colnames(design) <- gsub("group", "", colnames(design))

		contr.matrix <- makeContrasts(
		DE <- x \|>
		voomWithQualityWeights(design, plot = F) \|>
		lmFit(design) \|>
		contrasts.fit(
		contrasts = makeContrasts(
		contrasts = compare,
		levels = colnames(design)
		)

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

		if("logFC" %in% colnames(dt)){
		setnames(dt,"logFC","log2FC")
		setnames(DE,"logFC","log2FC")
		}

		dt
		DE
		}else if(re == "PCA"){
		x \|>
		cpm(log = T) \|>
		limma::plotMDS()
		}else if(re == "PCA_DATA"){
		x \|>
		cpm(log = T) \|>
		limma::plotMDS(plot = F)
		}
		}

		## the input data for venn_plot should come from rnaseq123_multi()

Admin message