Skip to content

Commit fc6f91e

Browse files
blex-maxblex-max
committed
Merge branch 'hotfix/3.5.5'
2 parents f2a469a + 9026421 commit fc6f91e

File tree

3 files changed

+175
-168
lines changed

3 files changed

+175
-168
lines changed

‎.gitignore‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1,2 @@
1-
bugfix-scratch/
1+
bugfix-scratch/
2+
bug-data/

‎R/nanoseq_results_plotter.R‎

Lines changed: 167 additions & 161 deletions
Original file line numberDiff line numberDiff line change
@@ -148,6 +148,11 @@ n_variants <- burdens[ismasked == 0][isvariant == 1]$count
148148
n_reference <- burdens[ismasked == 0][isvariant == 0]$count
149149
n_unique <- nrow(unique_variants[ismasked == 0])
150150

151+
if(length(n_reference) == 0) {
152+
message("0 reference calls. Exiting…\n")
153+
quit(save = "no", status = 1)
154+
}
155+
151156
if (length(n_variants) == 0) { n_variants = 0; }
152157
metrics <- data.frame("Metric" = c("total variants", "unique variants",
153158
"total variant + reference", "uncorrected burden", "physical coverage"),
@@ -568,181 +573,182 @@ write.table(tosave,file=paste(out_name, ".obs_burdens.pre_vs_post_masking.tsv",
568573
# Mismatches
569574
# To understand error profile - damage during library preparation
570575
mismatches = read.table(paste(dirname, 'mismatches.csv', sep = "/"), sep = ",", header = T)
571-
mismatches$sub = paste(unlist(sapply(mismatches$mismatch, function(x) unlist(strsplit(x, ""))[2])), ">", unlist(sapply(mismatches$mismatch, function(x) unlist(strsplit(x, ""))[5])), sep = "")
572-
mismatches = mismatches[which(mismatches$ismasked == 0),]
573-
mismatches = mismatches[which(mismatches$mismatch != "UNKNOWN"),]
574-
total = nrow(mismatches)
575576

576-
pdf(width = 5, height = 5, file = paste(out_name, ".mismatches.subst_asym.pdf", sep = ""))
577-
578-
subs = c("C>A", "G>T", "C>G", "G>C", "C>T", "G>A", "T>A", "A>T", "T>C", "A>G", "T>G", "A>C")
579-
colors = rep(c(rgb(34 / 255, 159 / 255, 198 / 255), rgb(26 / 255, 26 / 255, 26 / 255),
580-
rgb(201 / 255, 93 / 255, 94 / 255), rgb(178 / 255, 182 / 255, 180 / 255),
581-
rgb(153 / 255, 208 / 255, 62 / 255), rgb(217 / 255, 190 / 255, 217 / 255)), each = 2)
582-
par(mar = c(5, 5, 5, 5))
583-
counts = table(mismatches$sub)
584-
tmp_ = names(counts)
585-
counts = as.vector(counts)
586-
names(counts) = tmp_
587-
counts[setdiff(subs, names(counts))] = 0
588-
counts = counts[intersect(subs, names(counts))]
589-
counts = counts[subs]
590-
write.table(counts, file = paste(out_name, ".mismatches.subst_asym.tsv", sep = ""), quote = F, sep = "\t", row.names = F, col.names = F)
591-
592-
tick_10 = 10 * total / 100
593-
bar = barplot(counts, las = 2, col = colors,
594-
xlab = "Substitution (read strand)", ylab = "Number of mutations", border = "NA",
595-
ylim = c(0, max(counts) + 0.2 * max(counts)))
596-
axis(side = 4, at = c(0, tick_10, tick_10 * 2, tick_10 * 3), labels = c("0", "10", "20", "30"))
597-
mtext("% of mutations", side = 4, line = 2)
598-
# Add significance indicators
599-
pvalues = vector()
600-
for (i in c(1, 3, 5, 7, 9, 11)) {
601-
j = i + 1
602-
if (counts[i] + counts[j] > 0) {
603-
pvalues[subs[i]] = binom.test(counts[i], counts[i] + counts[j])$p.value
604-
} else {
605-
pvalues[subs[i]] = 1
606-
}
607-
}
608-
write.table(pvalues, file = paste(out_name, ".mismatches.subst_asym.pvals", sep = ""), quote = F, sep = "\t", row.names = T, col.names = F)
609-
qvalues = p.adjust(pvalues, method = "BH")
610-
for (i in c(1, 3, 5, 7, 9, 11)) {
611-
indicator = ""
612-
if (qvalues[subs[i]] < 0.01) {
613-
indicator = "**"
614-
} else if (qvalues[subs[i]] < 0.1) {
615-
indicator = "*"
577+
if(nrow(mismatches) > 0) {
578+
mismatches$sub = paste(unlist(sapply(mismatches$mismatch, function(x) unlist(strsplit(x, ""))[2])), ">", unlist(sapply(mismatches$mismatch, function(x) unlist(strsplit(x, ""))[5])), sep = "")
579+
mismatches = mismatches[which(mismatches$ismasked == 0),]
580+
mismatches = mismatches[which(mismatches$mismatch != "UNKNOWN"),]
581+
total = nrow(mismatches)
582+
583+
pdf(width = 5, height = 5, file = paste(out_name, ".mismatches.subst_asym.pdf", sep = ""))
584+
585+
subs = c("C>A", "G>T", "C>G", "G>C", "C>T", "G>A", "T>A", "A>T", "T>C", "A>G", "T>G", "A>C")
586+
colors = rep(c(rgb(34 / 255, 159 / 255, 198 / 255), rgb(26 / 255, 26 / 255, 26 / 255),
587+
rgb(201 / 255, 93 / 255, 94 / 255), rgb(178 / 255, 182 / 255, 180 / 255),
588+
rgb(153 / 255, 208 / 255, 62 / 255), rgb(217 / 255, 190 / 255, 217 / 255)), each = 2)
589+
par(mar = c(5, 5, 5, 5))
590+
counts = table(mismatches$sub)
591+
tmp_ = names(counts)
592+
counts = as.vector(counts)
593+
names(counts) = tmp_
594+
counts[setdiff(subs, names(counts))] = 0
595+
counts = counts[intersect(subs, names(counts))]
596+
counts = counts[subs]
597+
write.table(counts, file = paste(out_name, ".mismatches.subst_asym.tsv", sep = ""), quote = F, sep = "\t", row.names = F, col.names = F)
598+
599+
tick_10 = 10 * total / 100
600+
bar = barplot(counts, las = 2, col = colors,
601+
xlab = "Substitution (read strand)", ylab = "Number of mutations", border = "NA",
602+
ylim = c(0, max(counts) + 0.2 * max(counts)))
603+
axis(side = 4, at = c(0, tick_10, tick_10 * 2, tick_10 * 3), labels = c("0", "10", "20", "30"))
604+
mtext("% of mutations", side = 4, line = 2)
605+
# Add significance indicators
606+
pvalues = vector()
607+
for (i in c(1, 3, 5, 7, 9, 11)) {
608+
j = i + 1
609+
if (counts[i] + counts[j] > 0) {
610+
pvalues[subs[i]] = binom.test(counts[i], counts[i] + counts[j])$p.value
611+
} else {
612+
pvalues[subs[i]] = 1
613+
}
616614
}
617-
text((bar[i] + bar[i + 1]) / 2, max(counts[i], counts[i + 1]), indicator, pos = 3)
618-
}
619-
dev.off()
620-
621-
#########
622-
# Now the trinuc profiles
623-
pdf(width = 15, height = 8, file = paste(out_name, ".mismatches.trinuc-profile.pdf", sep = ""))
624-
par(mfrow = c(2, 1))
625-
par(mar = c(4, 6, 2, 2))
626-
627-
complement_subs2pyr = function(sub) {
628-
first = unlist(strsplit(sub, ""))[1]
629-
from = unlist(strsplit(sub, ""))[2]
630-
third = unlist(strsplit(sub, ""))[3]
631-
to = unlist(strsplit(sub, ""))[5]
632-
complement = vector()
633-
complement[c("A", "C", "G", "T")] = c("T", "G", "C", "A")
634-
pyr = c("C", "T")
635-
if (from %in% pyr) {
636-
return(sub)
637-
} else {
638-
sub = paste(complement[third], complement[from], complement[first], ">", complement[to], sep = "")
639-
return(sub)
615+
write.table(pvalues, file = paste(out_name, ".mismatches.subst_asym.pvals", sep = ""), quote = F, sep = "\t", row.names = T, col.names = F)
616+
qvalues = p.adjust(pvalues, method = "BH")
617+
for (i in c(1, 3, 5, 7, 9, 11)) {
618+
indicator = ""
619+
if (qvalues[subs[i]] < 0.01) {
620+
indicator = "**"
621+
} else if (qvalues[subs[i]] < 0.1) {
622+
indicator = "*"
623+
}
624+
text((bar[i] + bar[i + 1]) / 2, max(counts[i], counts[i + 1]), indicator, pos = 3)
640625
}
641-
}
626+
dev.off()
642627

643-
mismatches$mismatch_pyr = sapply(mismatches$mismatch, function(x) complement_subs2pyr(x))
628+
#########
629+
# Now the trinuc profiles
630+
pdf(width = 15, height = 8, file = paste(out_name, ".mismatches.trinuc-profile.pdf", sep = ""))
631+
par(mfrow = c(2, 1))
632+
par(mar = c(4, 6, 2, 2))
644633

645-
colours = rep(c("deepskyblue", "black", "firebrick2", "gray", "darkolivegreen3", "rosybrown2"), each = 16)
646-
sub_vec = c("C>A", "C>G", "C>T", "T>A", "T>C", "T>G")
647-
ctx_vec = paste(rep(c("A", "C", "G", "T"), each = 4), rep(c("A", "C", "G", "T"), times = 4), sep = "-")
648-
full_vec = paste(rep(sub_vec, each = 16), rep(ctx_vec, times = 6), sep = ",")
649-
xstr = paste(substr(full_vec, 5, 5), substr(full_vec, 1, 1), substr(full_vec, 7, 7), sep = "")
650-
ordered_names = paste(xstr, ">", rep(c("A", "G", "T", "A", "C", "G"), each = 16), sep = "")
651-
tmp_ = table(mismatches[which(mismatches$mismatch == mismatches$mismatch_pyr), "mismatch"])
652-
tri_obs = as.vector(tmp_)
653-
names(tri_obs) = names(tmp_)
654-
tri_obs[setdiff(ordered_names, names(tri_obs))] = 0
655-
tri_obs = tri_obs[ordered_names]
634+
complement_subs2pyr = function(sub) {
635+
first = unlist(strsplit(sub, ""))[1]
636+
from = unlist(strsplit(sub, ""))[2]
637+
third = unlist(strsplit(sub, ""))[3]
638+
to = unlist(strsplit(sub, ""))[5]
639+
complement = vector()
640+
complement[c("A", "C", "G", "T")] = c("T", "G", "C", "A")
641+
pyr = c("C", "T")
642+
if (from %in% pyr) {
643+
return(sub)
644+
} else {
645+
sub = paste(complement[third], complement[from], complement[first], ">", complement[to], sep = "")
646+
return(sub)
647+
}
648+
}
656649

657-
# Observed counts
658-
y = tri_obs;
659-
maxy = max(y)
660-
h = barplot(y, las = 2, col = colours, border = NA, ylim = c(0, maxy * 1.5), space = 1, cex.names = 0.6, names.arg = xstr, ylab = "Observed Pyrimidine mismatches counts", main = "Observed Pyirimidine mismatches counts")
661-
for (j in 1:length(sub_vec)) {
662-
xpos = h[c((j - 1) * 16 + 1, j * 16)]
663-
rect(xpos[1] - 0.5, maxy * 1.2, xpos[2] + 0.5, maxy * 1.3, border = NA, col = colours[j * 16])
664-
text(x = mean(xpos), y = maxy * 1.3, pos = 3, label = sub_vec[j])
665-
}
666-
write.table(tri_obs, file = paste(out_name, ".SSC-mismatches-Pyrimidine.triprofiles.tsv", sep = ""), quote = F, sep = "\t", row.names = T, col.names = F)
667-
tri_obs_pyr = tri_obs
650+
mismatches$mismatch_pyr = sapply(mismatches$mismatch, function(x) complement_subs2pyr(x))
668651

669-
tmp_ = table(mismatches[which(mismatches$mismatch != mismatches$mismatch_pyr), "mismatch_pyr"])
670-
tri_obs = as.vector(tmp_)
671-
names(tri_obs) = names(tmp_)
672-
tri_obs[setdiff(ordered_names, names(tri_obs))] = 0
673-
tri_obs = tri_obs[ordered_names]
652+
colours = rep(c("deepskyblue", "black", "firebrick2", "gray", "darkolivegreen3", "rosybrown2"), each = 16)
653+
sub_vec = c("C>A", "C>G", "C>T", "T>A", "T>C", "T>G")
654+
ctx_vec = paste(rep(c("A", "C", "G", "T"), each = 4), rep(c("A", "C", "G", "T"), times = 4), sep = "-")
655+
full_vec = paste(rep(sub_vec, each = 16), rep(ctx_vec, times = 6), sep = ",")
656+
xstr = paste(substr(full_vec, 5, 5), substr(full_vec, 1, 1), substr(full_vec, 7, 7), sep = "")
657+
ordered_names = paste(xstr, ">", rep(c("A", "G", "T", "A", "C", "G"), each = 16), sep = "")
658+
tmp_ = table(mismatches[which(mismatches$mismatch == mismatches$mismatch_pyr), "mismatch"])
659+
tri_obs = as.vector(tmp_)
660+
names(tri_obs) = names(tmp_)
661+
tri_obs[setdiff(ordered_names, names(tri_obs))] = 0
662+
tri_obs = tri_obs[ordered_names]
674663

675-
# Observed counts
676-
y = tri_obs;
677-
maxy = max(y)
678-
h = barplot(y, las = 2, col = colours, border = NA, ylim = c(0, maxy * 1.5), space = 1, cex.names = 0.6, names.arg = xstr, ylab = "Observed Purine mismatches counts", main = "Observed Purine mismatches counts")
679-
for (j in 1:length(sub_vec)) {
680-
xpos = h[c((j - 1) * 16 + 1, j * 16)]
681-
rect(xpos[1] - 0.5, maxy * 1.2, xpos[2] + 0.5, maxy * 1.3, border = NA, col = colours[j * 16])
682-
text(x = mean(xpos), y = maxy * 1.3, pos = 3, label = sub_vec[j])
683-
}
684-
dev.off()
685-
write.table(tri_obs, file = paste(out_name, ".SSC-mismatches-Purine.triprofiles.tsv", sep = ""), quote = F, sep = "\t", row.names = T, col.names = F)
686-
tri_obs = tri_obs_pyr + tri_obs
687-
write.table(tri_obs, file = paste(out_name, ".SSC-mismatches-Both.triprofiles.tsv", sep = ""), quote = F, sep = "\t", row.names = T, col.names = F)
664+
# Observed counts
665+
y = tri_obs;
666+
maxy = max(y)
667+
h = barplot(y, las = 2, col = colours, border = NA, ylim = c(0, maxy * 1.5), space = 1, cex.names = 0.6, names.arg = xstr, ylab = "Observed Pyrimidine mismatches counts", main = "Observed Pyirimidine mismatches counts")
668+
for (j in 1:length(sub_vec)) {
669+
xpos = h[c((j - 1) * 16 + 1, j * 16)]
670+
rect(xpos[1] - 0.5, maxy * 1.2, xpos[2] + 0.5, maxy * 1.3, border = NA, col = colours[j * 16])
671+
text(x = mean(xpos), y = maxy * 1.3, pos = 3, label = sub_vec[j])
672+
}
673+
write.table(tri_obs, file = paste(out_name, ".SSC-mismatches-Pyrimidine.triprofiles.tsv", sep = ""), quote = F, sep = "\t", row.names = T, col.names = F)
674+
tri_obs_pyr = tri_obs
688675

676+
tmp_ = table(mismatches[which(mismatches$mismatch != mismatches$mismatch_pyr), "mismatch_pyr"])
677+
tri_obs = as.vector(tmp_)
678+
names(tri_obs) = names(tmp_)
679+
tri_obs[setdiff(ordered_names, names(tri_obs))] = 0
680+
tri_obs = tri_obs[ordered_names]
689681

690-
##########################################################################################
691-
# Now error rates based on SSC, by each of 96 trinuc-sub channels and global
692-
tmp_ = table(mismatches[which(mismatches$mismatch == mismatches$mismatch_pyr), "mismatch"])
693-
tri_obs = as.vector(tmp_)
694-
names(tri_obs) = names(tmp_)
695-
tri_obs[setdiff(ordered_names, names(tri_obs))] = 0
696-
tri_obs = tri_obs[ordered_names]
697-
tri_obs_pyr = tri_obs
698-
tmp_ = table(mismatches[which(mismatches$mismatch != mismatches$mismatch_pyr), "mismatch_pyr"])
699-
tri_obs = as.vector(tmp_)
700-
names(tri_obs) = names(tmp_)
701-
tri_obs[setdiff(ordered_names, names(tri_obs))] = 0
702-
tri_obs = tri_obs[ordered_names]
703-
tri_obs_pur = tri_obs
682+
# Observed counts
683+
y = tri_obs;
684+
maxy = max(y)
685+
h = barplot(y, las = 2, col = colours, border = NA, ylim = c(0, maxy * 1.5), space = 1, cex.names = 0.6, names.arg = xstr, ylab = "Observed Purine mismatches counts", main = "Observed Purine mismatches counts")
686+
for (j in 1:length(sub_vec)) {
687+
xpos = h[c((j - 1) * 16 + 1, j * 16)]
688+
rect(xpos[1] - 0.5, maxy * 1.2, xpos[2] + 0.5, maxy * 1.3, border = NA, col = colours[j * 16])
689+
text(x = mean(xpos), y = maxy * 1.3, pos = 3, label = sub_vec[j])
690+
}
691+
dev.off()
692+
write.table(tri_obs, file = paste(out_name, ".SSC-mismatches-Purine.triprofiles.tsv", sep = ""), quote = F, sep = "\t", row.names = T, col.names = F)
693+
tri_obs = tri_obs_pyr + tri_obs
694+
write.table(tri_obs, file = paste(out_name, ".SSC-mismatches-Both.triprofiles.tsv", sep = ""), quote = F, sep = "\t", row.names = T, col.names = F)
704695

705-
pdf(width = 15, height = 4, file = paste(out_name, ".DSC_estimated_error_rates.pdf", sep = ""))
706-
par(mar = c(4, 6, 2, 2))
707-
tris = sapply(names(tri_obs), function(x) unlist(strsplit(x, ">"))[1])
708-
y = (tri_obs_pyr / (tri_bg[tris] / 2)) * (tri_obs_pur / (tri_bg[tris] / 2));
709-
maxy = max(y)
710-
##y = (tri_obs_pyr/(tri_bg[tris])) * (tri_obs_pur/(tri_bg[tris])); maxy = max(y)
711-
h = barplot(y, las = 2, col = colours, border = NA, ylim = c(0, maxy * 1.5), space = 1, cex.names = 0.6, names.arg = xstr, ylab = "Error rates", main = "Error rates")
712-
for (j in 1:length(sub_vec)) {
713-
xpos = h[c((j - 1) * 16 + 1, j * 16)]
714-
rect(xpos[1] - 0.5, maxy * 1.2, xpos[2] + 0.5, maxy * 1.3, border = NA, col = colours[j * 16])
715-
text(x = mean(xpos), y = maxy * 1.3, pos = 3, label = sub_vec[j])
716-
}
717-
dev.off()
718696

697+
##########################################################################################
698+
# Now error rates based on SSC, by each of 96 trinuc-sub channels and global
699+
tmp_ = table(mismatches[which(mismatches$mismatch == mismatches$mismatch_pyr), "mismatch"])
700+
tri_obs = as.vector(tmp_)
701+
names(tri_obs) = names(tmp_)
702+
tri_obs[setdiff(ordered_names, names(tri_obs))] = 0
703+
tri_obs = tri_obs[ordered_names]
704+
tri_obs_pyr = tri_obs
705+
tmp_ = table(mismatches[which(mismatches$mismatch != mismatches$mismatch_pyr), "mismatch_pyr"])
706+
tri_obs = as.vector(tmp_)
707+
names(tri_obs) = names(tmp_)
708+
tri_obs[setdiff(ordered_names, names(tri_obs))] = 0
709+
tri_obs = tri_obs[ordered_names]
710+
tri_obs_pur = tri_obs
719711

720-
# Global:
721-
# Number of errors per channel:
722-
errors_per_channel = y * tri_bg[tris]
723-
total_error_rate = sum(errors_per_channel) / sum(tri_bg)
712+
pdf(width = 15, height = 4, file = paste(out_name, ".DSC_estimated_error_rates.pdf", sep = ""))
713+
par(mar = c(4, 6, 2, 2))
714+
tris = sapply(names(tri_obs), function(x) unlist(strsplit(x, ">"))[1])
715+
y = (tri_obs_pyr / (tri_bg[tris] / 2)) * (tri_obs_pur / (tri_bg[tris] / 2));
716+
maxy = max(y)
717+
##y = (tri_obs_pyr/(tri_bg[tris])) * (tri_obs_pur/(tri_bg[tris])); maxy = max(y)
718+
h = barplot(y, las = 2, col = colours, border = NA, ylim = c(0, maxy * 1.5), space = 1, cex.names = 0.6, names.arg = xstr, ylab = "Error rates", main = "Error rates")
719+
for (j in 1:length(sub_vec)) {
720+
xpos = h[c((j - 1) * 16 + 1, j * 16)]
721+
rect(xpos[1] - 0.5, maxy * 1.2, xpos[2] + 0.5, maxy * 1.3, border = NA, col = colours[j * 16])
722+
text(x = mean(xpos), y = maxy * 1.3, pos = 3, label = sub_vec[j])
723+
}
724+
dev.off()
724725

725-
pdf(width = 15, height = 4, file = paste(out_name, ".DSC_errors_per_channel.pdf", sep = ""))
726-
par(mar = c(4, 6, 2, 2))
727-
tris = sapply(names(tri_obs), function(x) unlist(strsplit(x, ">"))[1])
728-
y = errors_per_channel;
729-
maxy = max(y)
730-
h = barplot(y, las = 2, col = colours, border = NA, ylim = c(0, maxy * 1.5), space = 1, cex.names = 0.6, names.arg = xstr, ylab = "Predicted errors", main = "Predicted errors")
731-
for (j in 1:length(sub_vec)) {
732-
xpos = h[c((j - 1) * 16 + 1, j * 16)]
733-
rect(xpos[1] - 0.5, maxy * 1.2, xpos[2] + 0.5, maxy * 1.3, border = NA, col = colours[j * 16])
734-
text(x = mean(xpos), y = maxy * 1.3, pos = 3, label = sub_vec[j])
735-
}
736-
dev.off()
737726

738-
output = list()
739-
output[["errors_per_channel"]] = errors_per_channel
740-
output[["error_rates_per_channel"]] = y
741-
output[["total_error_rate"]] = total_error_rate
742-
output[["total_errors"]] = sum(errors_per_channel)
743-
kkk = vector()
744-
kkk["total_error_rate"] = total_error_rate
745-
kkk["total_errors"] = output[["total_errors"]]
746-
write.table(kkk, file = paste(out_name, ".estimated_error_rates.tsv", sep = ""), quote = F, sep = "\t", row.names = T, col.names = F)
727+
# Global:
728+
# Number of errors per channel:
729+
errors_per_channel = y * tri_bg[tris]
730+
total_error_rate = sum(errors_per_channel) / sum(tri_bg)
747731

732+
pdf(width = 15, height = 4, file = paste(out_name, ".DSC_errors_per_channel.pdf", sep = ""))
733+
par(mar = c(4, 6, 2, 2))
734+
tris = sapply(names(tri_obs), function(x) unlist(strsplit(x, ">"))[1])
735+
y = errors_per_channel;
736+
maxy = max(y)
737+
h = barplot(y, las = 2, col = colours, border = NA, ylim = c(0, maxy * 1.5), space = 1, cex.names = 0.6, names.arg = xstr, ylab = "Predicted errors", main = "Predicted errors")
738+
for (j in 1:length(sub_vec)) {
739+
xpos = h[c((j - 1) * 16 + 1, j * 16)]
740+
rect(xpos[1] - 0.5, maxy * 1.2, xpos[2] + 0.5, maxy * 1.3, border = NA, col = colours[j * 16])
741+
text(x = mean(xpos), y = maxy * 1.3, pos = 3, label = sub_vec[j])
742+
}
743+
dev.off()
748744

745+
output = list()
746+
output[["errors_per_channel"]] = errors_per_channel
747+
output[["error_rates_per_channel"]] = y
748+
output[["total_error_rate"]] = total_error_rate
749+
output[["total_errors"]] = sum(errors_per_channel)
750+
kkk = vector()
751+
kkk["total_error_rate"] = total_error_rate
752+
kkk["total_errors"] = output[["total_errors"]]
753+
write.table(kkk, file = paste(out_name, ".estimated_error_rates.tsv", sep = ""), quote = F, sep = "\t", row.names = T, col.names = F)
754+
}

‎python/runNanoSeq.py‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -423,9 +423,9 @@ def runCommand(command):
423423
print("\nExecuting: %s\n" % ijob)
424424
p = subprocess.Popen(
425425
ijob, shell=True, stdout=subprocess.DEVNULL, stderr=subprocess.PIPE)
426-
p.wait()
426+
std: tuple = p.communicate()
427427
if (p.returncode != 0):
428-
error = p.stderr.read().decode()
428+
error = std[1].decode()
429429
sys.stderr.write("\n!Error processing: %s\n" % ijob)
430430
raise ValueError(error)
431431
return
@@ -1205,10 +1205,10 @@ def vcfHeader(args):
12051205
ifile = tmpDir+"/indel/%s.indel.filtered.vcf.gz" % (i+1)
12061206
if ( os.stat(ifile).st_size == 0 ) : continue
12071207
vcf2Merge.append(ifile)
1208-
cmd = "cp %s merged.vcf.gz ;"% vcf2Merge.pop(0)
1209-
for ifile in vcf2Merge[1:]:
1210-
cmd += " bcftools concat --no-version -Oz -o tmp.xxx.vcf.gz merged.vcf.gz %s; mv tmp.xxx.vcf.gz merged.vcf.gz; "%ifile
1211-
cmd += "bcftools sort -Oz -o %s/post/%s.indel.vcf.gz merged.vcf.gz;" % (tmpDir, args.name)
1208+
cmd = "cp %s %s/indel/merged.vcf.gz ;" % (vcf2Merge.pop(0), tmpDir)
1209+
for ifile in vcf2Merge[0:]:
1210+
cmd += "bcftools concat --no-version -Oz -o %s/indel/tmp.xxx.vcf.gz %s/indel/merged.vcf.gz %s; mv %s/indel/tmp.xxx.vcf.gz %s/indel/merged.vcf.gz; " % (tmpDir, tmpDir, ifile, tmpDir, tmpDir)
1211+
cmd += "bcftools sort -Oz -o %s/post/%s.indel.vcf.gz %s/indel/merged.vcf.gz;" % (tmpDir, args.name, tmpDir)
12121212
cmd += "bcftools index -t -f %s/post/%s.indel.vcf.gz " % (
12131213
tmpDir, args.name)
12141214
runCommand(cmd)

0 commit comments

Comments
 (0)