Skip to content
Snippets Groups Projects
Commit a604f782 authored by Maximilian Press's avatar Maximilian Press
Browse files

Small changes fixing PE coverage

I am trying to identify a minimal set of code changes that correct outputs in the statistics.csv file. I think it might be the following:

1. apply utilities.correctEdge() to paired_whole_coverage array inside of readsCoverage().

2. Make statistics.csv report _always_ write out paired_whole_coverage rather than whole_coverage. They should be equivalent in single-end case, and whole_coverage is simply wrong in PE case. So just use paired_whole_coverage always, don't bother testing for it.

Trivially, (2) might fix the problem by just making PT report the coverage that it is actually using, if the conditional there is causing a problem. But since the statistics.csv file is our source of truth, it looks very wrong. Maybe this will fix it.
parent ae2c82fc
No related branches found
No related tags found
No related merge requests found
......@@ -340,7 +340,7 @@ def readsCoverage(inRawDArgs,refseq,inDArgs,fParms,return_dict, core_id,line_sta
print("WARNING: No read Match, please check your fastq file")
termini_coverage = correctEdge(termini_coverage, fParms.edge)
#paired_whole_coverage = correctEdge(whole_coverage, fParms.edge) #TODO: discuss with Julian and Max about the PE issue that Max reported.
paired_whole_coverage = correctEdge(whole_coverage, fParms.edge) #TODO: discuss with Julian and Max about the PE issue that Max reported.
whole_coverage = correctEdge(whole_coverage, fParms.edge)
phage_hybrid_coverage = correctEdge(phage_hybrid_coverage, fParms.edge)
if inDArgs.hostseq != "":
......@@ -694,19 +694,15 @@ def ExportStatistics(phagename, whole_coverage, paired_whole_coverage, termini_c
export = pd.DataFrame()
# ORGANIZE Column
export["Position"] = list(phage_plus_norm.sort_values("Position")["Position"])
if paired != "":
export["Coverage +"] = paired_whole_coverage[0]
else:
export["Coverage +"] = whole_coverage[0]
# coverage values for whole_coverage == paired_whole_coverage in the single-end case
export["Coverage +"] = paired_whole_coverage[0]
export["SPC +"] = termini_coverage[0]
export["T +"] = [format(x/100.0,'0.2') for x in list(phage_plus_norm.sort_values("Position")["SPC_std"])]
export["T + (close)"] = [format(x/100.0,'0.2') for x in list(phage_plus_norm.sort_values("Position")["SPC"])]
export["pvalue +"] = [format(x,'0.2e') for x in list(phage_plus_norm.sort_values("Position")["pval_gamma"])]
export["padj +"] = [format(x,'0.2e') for x in list(phage_plus_norm.sort_values("Position")["pval_gamma_adj"])]
if paired != "":
export["Coverage -"] = whole_coverage[1]
else:
export["Coverage -"] = paired_whole_coverage[1]
# coverage values for whole_coverage == paired_whole_coverage in the single-end case
export["Coverage -"] = paired_whole_coverage[1]
export["SPC -"] = termini_coverage[1]
export["T -"] = [format(x/100.0,'0.2') for x in list(phage_minus_norm.sort_values("Position")["SPC_std"])]
export["T - (close)"] = [format(x/100.0,'0.2') for x in list(phage_minus_norm.sort_values("Position")["SPC"])]
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment