tRAX/readlengthhistogram.R at master · UCSC-LoweLab/tRAX · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78


library(ggplot2)
library(reshape2)
library(scales)

args <- commandArgs(trailingOnly = TRUE)

#args <- c("hg19-genomictrailertable.txt","hg19-genomicbarplot.png")
#args =c("repfragtypenormcounts.txt","barplot.png")


#Rscript trailerbarplot.R hg19-trailertable.txt hg19-barplot.png
readlengths <- read.table(args[1], header = TRUE)
sampledata <- read.table(args[2])


#temp = cbind(selectcounts, seq = factor(rownames(selectcounts),rev(rownames(selectcounts)), ordered = TRUE))

#levels(temp$seq) <- rev(rownames(selectcounts))

#countsmelt = melt(temp, id.vars = c('seq'))


if(FALSE){
ggplot(readlengths,aes(x = variable, y = value,fill = seq, stat="identity")) +
	geom_bar(position = "fill",stat="identity") +
    geom_bar(position = "fill",stat="identity",color="black",show_guide=FALSE) +
    scale_y_continuous(labels = percent_format()) +
    xlab("Sample") +
    ylab("Percentage of Reads") +
    theme(axis.title.x = element_text(face="bold", size=15), axis.text.x = element_text(face="bold", size=9,angle = 90, vjust = .5))
    }


binw <- 2
maxlen <- min(c(200, max(readlengths[,1])))
readlengths$bin <- floor(readlengths$Length / binw)*binw

#head(readlengths)
binnedtrnareadlengths <- aggregate(trnas ~ bin * Sample, readlengths, sum)

#binnedreadlengths <- aggregate(other ~ bin, readlengths, sum)
binnedreadlengths <- aggregate(other ~ bin * Sample ,readlengths, sum)
#head(binnedreadlengths)
#unique(binnedreadlengths$Sample)
temp = merge(binnedtrnareadlengths,binnedreadlengths,by=c("bin", "Sample"))
#head(temp)
#temp
#temp = binnedtrnareadlengths
binnedlengths = melt(temp, id.vars = c('bin', 'Sample'))
#binnedlengths = temp
#head(binnedlengths)
binnedlengths[is.na(binnedlengths)] <- 0
#unique(binnedlengths$variable)
binnedlengths <- binnedlengths[order(binnedlengths$Sample,binnedlengths$bin),]

#binnedlengths <- binnedlengths[binnedlengths$variable =='trnas',]

#binnedlengths
#binnedlengths$variable <- factor(binnedlengths$variable, levels = c("trnas", "other"))
binnedlengths$variable <- factor(binnedlengths$variable, levels = c("other", "trnas"))
binnedlengths$Sample <- factor(binnedlengths$Sample, levels = sampledata[,1])

#levels(binnedreadlengths$Sample)
#levels(sampledata[,1])

ggplot(data=binnedlengths, aes(x = bin, y = value, fill = variable))+geom_bar(stat="identity")+facet_wrap( ~ Sample, scales="free",ncol = 3)+ xlim(0,maxlen) + ylab("Count") + xlab("Read Length") + scale_fill_discrete(name="Gene Type")  + scale_x_continuous(labels = comma, limits=c(0, maxlen)) + theme_bw()
#length(unique(binnedlengths$Sample))
ggsave(filename=args[3],limitsize=FALSE,width = 8, height = .5 * length(unique(binnedlengths$Sample)))#, width = 3 * length(unique(binnedlengths$Sample)))