| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657 |
- # parse parameter ---------------------------------------------------------
- library(argparser, quietly=TRUE)
- # Create a parser
- p <- arg_parser("expression matrix statistics")
- # Add command line arguments
- p <- add_argument(p, "--dat", help="expression matrix", type="character")
- p <- add_argument(p, "--out", help="output prefix", type="character")
- # Parse the command line arguments
- argv <- parse_args(p)
- dat<-argv$dat
- outprefix<-argv$out
- library(ggplot2)
- library(tidyverse)
- library(formattable)
- ## read data and melt
- expr <- read_delim(dat,
- "\t", escape_double = FALSE, trim_ws = TRUE)
- expr <- rename(expr, featureId = X1)
- expr_tidy<-gather(expr, key = "sample", value = "expression", -1)
- # density
- f<-paste(outprefix, "density.pdf", sep = ".")
- p<-qplot(log10(expression), data=expr_tidy, geom="density", xlim=c(-2,6), colour=sample, fill=sample, alpha=I(1/10), main = "Density distribution of transcript expression") +
- theme_bw() +
- scale_y_continuous(expand = c(0,0))
- pdf(file=f)
- p
- dev.off()
- # boxplot
- p<-ggplot(expr_tidy, aes(x=sample, y=log10(expression), fill=sample)) + geom_boxplot() +
- theme_bw() +
- theme(axis.text.x=element_text(hjust=1,angle=60)) + scale_y_continuous(expand = c(0,0))
- f<-paste(outprefix, "boxplot.pdf", sep = ".")
- pdf(file=f)
- p
- dev.off()
- # expressed gene
- total_number = length(expr$featureId)
- expressed_tidy <- expr_tidy %>% filter(expression > 1) %>%
- group_by(sample) %>%
- summarise(expressed = n()) %>%
- mutate(percentage = percent(expressed/total_number))
- f<-paste(outprefix, "stat.txt", sep = ".")
- write.table(expressed_tidy, file = f, quote = F, row.names = F)
|
