/
StatQualFiles.R
90 lines (80 loc) · 2.52 KB
/
StatQualFiles.R
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
79
80
81
82
83
84
85
86
87
88
89
########################################################################################
## This file is a part of YAP package of scripts. https://github.com/shpakoo/YAP
## Distributed under the MIT license: http://www.opensource.org/licenses/mit-license.php
## Copyright (c) 2011-2013 Sebastian Szpakowski
########################################################################################
makePlots = function(tmp, dataset, f)
{
for ( stat in c("length", "meanq", "medq", "R50.meanq", "cov"))
{
print (stat)
if ( max(dataset[,stat])>10 )
{
br = seq(0-1.5, max(dataset[,stat])+1.5, by=1)
}
else
{
br = 100
}
#print (br)
hist (tmp[,stat], br = br, main=paste(f, stat, sep="\n"), xlim=c(0, max(dataset[,stat]) ) )
}
x=sample(1:nrow(tmp), 0.1* nrow(tmp) )
plot (tmp$length[x], tmp$meanq[x], pch=19, col="gray", cex=0.1, xlim=c(0, max(dataset[,"length"]) ), ylim=c(0, max(dataset[,"meanq"]) ))
points(lowess(tmp$length, tmp$meanq, f=0.05) , type="l", lwd=6, col="gray50")
abline (v= mean(tmp[,"length"]), col="green", lwd=3)
abline (h= mean(tmp[,"meanq"]), col="blue", lwd=3)
legend ("bottom", legend=c(round(mean(tmp[,"length"]),0), round(mean(tmp[,"meanq"]),0), nrow(tmp)), fill=c("green", "blue", "black"), cex=1.5 , ncol=3)
}
processFiles = function(filelist, prefix, parentchild=list(), raster=TRUE)
{
dataset=data.frame()
for (filename in filelist)
{
cat("Loading ", filename,"...", sep="")
tmp = read.table(filename, header=TRUE, as.is=TRUE)
if (nrow(dataset)==0)
{
dataset=tmp
}
else
{
dataset=rbind(dataset, tmp)
}
cat("\n")
}
dataset$cov = dataset$sdq / dataset$meanq
files = unique(dataset$file)
panels = length(files)+length(parentchild)
if (raster)
{
png(paste(prefix, ".stats.png", sep=""), width=26, height=5*(panels), unit="in", res=400)
}
else
{
pdf(paste(prefix, ".stats.pdf", sep=""), width=24, height=5*(panels), paper="special")
}
par(mfrow= c(panels,6))
for (f in files)
{
tmp= dataset[dataset$file==f, ]
makePlots(tmp, dataset, f)
for (k in names(parentchild))
{
if (k == f)
{
cat("Loading ", parentchild[[k]], " for ", k, "...", sep="")
subsets = read.table(parentchild[[k]], header=TRUE, as.is=TRUE)
tmp = tmp[tmp$read %in% subsets$QueryName,]
makePlots(tmp, dataset, parentchild[[k]])
cat("\n")
}
}
}
dev.off()
invisible(dataset)
}
# lista=list()
# lista["T1.trim.qual"] = "T1.trim.unique.fragclust.align.report"
# filelist=c("T1.trim.qual.stats.tab.txt")
# processFiles(filelist, "test", lista)