modified RNA-seq.R in R (2f82902d) · Commits · Chaos / chaos.tools

R/RNA-seq.R

+29 −42

Original line number	Diff line number	Diff line
		@@ -4,22 +4,24 @@ rnaseq123 <- function(
		organism = "human",
		count_type = "star"
		){
		if(tolower(count_type) == "star")
		{
		x <- readDGE(
		c(ctl_count_file,obs_count_file),
		header = F,
		skip = 4
		)

		if(tolower(count_type) == "star"){
		cols <- c(1, 2)
		header_str <- TRUE
		skip_rows <- 4
		}else if(tolower(count_type) == "featurecounts")
		{
		cols <- c(1, 7)
		header_str <- FALSE
		skip_rows <- 1
		}

		x <- readDGE(
		c(ctl_count_file,obs_count_file),
		columns = c(1,7),
		sep = "\t",
		skip = 1
		columns = cols,
		header = header_str,
		skip = skip_rows
		)
		}

		group <- c(
		rep("ctl",length(ctl_count_file)),
		@@ -28,26 +30,6 @@ rnaseq123 <- function(

		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)] \|>
		@@ -480,28 +462,25 @@ rnaseqqc <- function(
		name = NA,
		group = NA,
		count_type = "star",
		plot_type = "mds"
		out_type = "mds"
		){
		g.number <- length(group)

		if(g.number != length(count_file)){
		stop()
		}else{
		g.color <- as.factor(group)
		if(g.number <= 10){
		levels(g.color) <- sample(chaos_color(),g.number)
		}else{
		levels(g.color) <- colorRampPalette(chaos_color())(n)
		levels(g.color) <- colorRampPalette(chaos_color())(g.number)
		}
		g.color <- as.character(g.color)
		}


		if(length(name) == 1){
		if(is.na(name)){
		name <- base_name(count_file,"ReadsPerGene.out.tab")
		}
		}


		if(tolower(count_type) == "star"){
		cols <- c(1, 2)
		header_str <- TRUE
		@@ -526,10 +505,18 @@ rnaseqqc <- function(
		calcNormFactors(method = "TMM") \|>
		cpm(log = T)

		if(plot_type == "mds"){
		limma::plotMDS(x, labels = sample,col = g.color)
		}else if(plot_type == "glmds"){
		Glimma::glMDSPlot(x, labels = sample, launch = F,groups = group)
		colnames(x) <- name

		if(out_type == "pca"){
		limma::plotMDS(x, labels = name,col = g.color)
		}else if(out_type == "pca_data"){
		limma::plotMDS(x, plot = F)
		}else if(out_type == "pca_html"){
		Glimma::glMDSPlot(x, labels = name, launch = F,groups = group)
		print("figure saved in ./glimma-plots")
		}else if(out_type == "cor_data"){
		cor(x)
		}else if(out_type == "cor"){
		cor(x) \|> pheatmap::pheatmap(display_numbers = T)
		}
		}
		No newline at end of file

Admin message