-
Notifications
You must be signed in to change notification settings - Fork 5
/
stats.R
30 lines (26 loc) · 1.33 KB
/
stats.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
#######
####### get some stats and some plots
####### takes the output directory as sole input
####### outputs two files: 1. a .txt file with ids and associated demultiplexed reads and 2. pdf file with barplot of demultiplexed datasets
args = commandArgs(trailingOnly = T)
if (length(args) != 1){
stop('Please provide fastQ output directory as solitary input')
} else if (length(args) == 1) output.dir = args
output.dir = '~/Barcode_Partitioning/test.cases/'
setwd(output.dir)
output.fastqs = list.files(path = output.dir, pattern = 'fastq')
output.table = data.frame(sample.ids = character(length(output.fastqs)), read.counts = numeric(length(output.fastqs)), stringsAsFactors = F)
for (i in seq(1, length(output.fastqs))){
command = paste0('wc -l ', output.fastqs[i])
test = system(command, intern = T)
results = unlist(strsplit(x = test, split = ' '))
ind = which(results == output.fastqs[i])
id = results[ind]
counts = as.numeric(results[ind-1])/4
output.table$sample.ids[i] = id
output.table$read.counts[i] = counts
}
write.table(x = output.table, file = 'demultiplexingResults.txt', append = F, quote = F, sep = '\t', row.names = F, col.names = T)
pdf('graphicalDemultiplexingResults.pdf')
barplot(height = output.table$read.counts, names.arg = output.table$sample.ids, ylab = '# reads/file', col = 'red', cex.names = 0.3)
dev.off()