58 lines
1.7 KiB
R
Executable File
58 lines
1.7 KiB
R
Executable File
#!/usr/bin/env Rscript
|
|
|
|
suppressMessages(library(dplyr))
|
|
suppressMessages(library(ggplot2))
|
|
suppressMessages(library(tidyr))
|
|
|
|
# Set up sample data frame:
|
|
coldata <- read.csv("de_analysis/coldata.tsv", row.names="alias", sep=",")
|
|
coldata$sample_id <- rownames(coldata)
|
|
coldata$condition <- factor(coldata$condition, levels=rev(levels(coldata$condition)))
|
|
coldata$type <-NULL
|
|
coldata$patient <-NULL
|
|
|
|
# Read stageR results:
|
|
stageR <- read.csv("de_analysis/results_dtu_stageR.tsv", sep="\t")
|
|
names(stageR) <- c("gene_id", "transcript_id", "p_gene", "p_transcript");
|
|
|
|
# Read filtered counts:
|
|
counts <- read.csv("merged/all_counts_filtered.tsv", sep="\t");
|
|
names(counts)[2]<-"transcript_id"
|
|
|
|
# Join counts and stageR results:
|
|
df <- counts %>% left_join(stageR, by = c("gene_id", "transcript_id"))
|
|
df <- df[order(df$p_gene),]
|
|
|
|
scols <- setdiff(names(df),c("gene_id", "transcript_id", "p_gene", "p_transcript"))
|
|
|
|
# Normalise counts:
|
|
for(sc in scols){
|
|
df[sc] <- df[sc] / sum(df[sc])
|
|
}
|
|
|
|
# Melt data frame:
|
|
tdf <- df %>% gather(key='sample', value='norm_count',-gene_id, -transcript_id, -p_gene, -p_transcript)
|
|
|
|
# Add sample group column:
|
|
sampleToGroup<-function(x){
|
|
return(coldata[x,]$condition)
|
|
}
|
|
|
|
tdf$group <- sampleToGroup(tdf$sample)
|
|
|
|
# Filter for significant genes:
|
|
sig_level <- 0.05
|
|
genes <- as.character(tdf[which(tdf$p_gene < sig_level),]$gene_id)
|
|
genes <- unique(genes)
|
|
|
|
pdf("de_analysis/dtu_plots.pdf")
|
|
|
|
for(gene in genes){
|
|
gdf<-tdf[which(tdf$gene_id==gene),]
|
|
p_gene <- unique(gdf$p_gene)
|
|
p <- ggplot(gdf, aes(x=transcript_id, y=norm_count)) + geom_bar(stat="identity", aes(fill=sample), position="dodge")
|
|
p <- p + facet_wrap(~ group) + coord_flip()
|
|
p <- p + ggtitle(paste(gene," : p_value=",p_gene,sep=""))
|
|
print(p)
|
|
}
|