Better output from MelSimCoverage

Track per-base coverage, rather than number of reads.

Output now includes positions in the table.

Author: petercombs

MelSimCoverage.py

@@ -1,4 +1,5 @@
 from progressbar import ProgressBar
+import pandas as pd
 try:
     import matplotlib.pyplot as mpl
     mpl.figure()
@@ -15,18 +16,21 @@
     mel = Samfile('analysis/on_mel/mel_gdna_bowtie2_dedup.bam')
     sim = Samfile('analysis/on_mel/sim_gdna_bowtie2_dedup.bam')
 
+    posns = []
     mel_covs = []
     sim_covs = []
     prog = 0
     pbar = ProgressBar(max_value=sum(mel.lengths) + 1)
     for r, l in zip(mel.references, mel.lengths):
         for start in range(0, l, window_size):
-            mel_covs.append(mel.count(r, start, start+window_size))
-            sim_covs.append(sim.count(r, start, start+window_size))
+            posns.append('{}_{}'.format(r, start))
+            mel_covs.append(sum(sum(i) for i in
+                                mel.count_coverage(r, start, start+window_size)))
+            sim_covs.append(sum(sum(i) for i in
+                                sim.count_coverage(r, start, start+window_size)))
         prog += l
         pbar.update(prog)
     pbar.finish()
-
-
-
-
+    (pd.DataFrame(index=posns, data={'mel':mel_covs, 'sim':sim_covs})
+     .to_csv('analysis/on_mel/melsim_cov.tsv', sep='\t')
+    )