Changes for distributions

main.py

@@ -3,19 +3,22 @@
 import argparse
 from tmb.config import Config
 from tmb.bedreader import BedReader
-from tmb.distribution import calculate_tmb
+from tmb.distribution import get_tmb_from_files
 from multiprocessing import Pool
+import logging
 
 
-#TODO: add logger
 def main():
     config = parse_config()
+    set_logging(config)
     exons = BedReader(config.bed).exonList
 
     pool = Pool(processes=Config.th)
-    results = [pool.apply_async(calculate_tmb, args=(config.tumor, config.normal, exon)) for exon in exons]
+    results = [pool.apply_async(get_tmb_from_files, args=(config.tumor, config.normal, exon)) for exon in exons]
+
+    if Config.header:
+        print_header()
 
-    print_header()
     for output in results:
         tmb = output.get()
         print(tmb)
@@ -34,6 +37,9 @@ def parse_config():
                         help="Minimum coverage of base for nucleotide to be considered as somatic. "
                              "Default: 2")
     optional.add_argument('--th', type=int, help="Number of threads for multiprocessing mode. Default: 1")
+    optional.add_argument('--pvalue', type=int, help="Pvalue threshold for somatic sites. Default: 5*10E20")
+    optional.add_argument('--log', type=str, help="The level of logging from CRITICAL, ERROR, WARNING, INFO, DEBUG. ")
+    optional.add_argument('--header', type=str, help="Print header to the output")
     required.add_argument('--normal', required=True, type=str, help="Path to normal SAM/BAM/CRAM file.")
     required.add_argument('--tumor', required=True, type=str, help="Path to tumor SAM/BAM/CRAM file.")
     required.add_argument('--bed', required=True, type=str, help="Path to BED file.")
@@ -46,5 +52,15 @@ def print_header():
     print("\t".join(['Chr', 'Start', 'End', 'Gene', 'SomaticSites', 'TMB']))
 
 
+def set_logging(config):
+    log_level = getattr(logging, Config.log_level.upper(), None)
+    if not isinstance(log_level, int):
+        raise ValueError('Invalid log level: %s' % log_level)
+    logging.basicConfig(filename='tmb.log', level=log_level, filemode='w',
+                        format='%(asctime)s %(levelname)s:%(message)s ')
+    logging.info("Starting TMB.")
+    logging.info("Retrieving exons from %s.", config.bed)
+
+
 if __name__ == "__main__":
     main()

tests/test_distribution.py

@@ -1,18 +1,24 @@
 from numpy.testing import assert_almost_equal
 
 from tmb.bedreader import BedReader
-from tmb.distribution import calculate_tmb
+from tmb.bamreader import BamReader
+from tmb.exome import Exome
+from tmb.distribution import get_tmb_from_files
+from tmb.distribution import calculate_tmb_from_sites
+from tmb.distribution import fill_maps_reads
+from tmb.distribution import fill_maps_pileup
+from tmb.distribution import hypotesis_test
 
 import pytest
 
 
-def test_distribution():
+def test_distribution_by_pileup():
     normal = '../data/chr5_665281/normal_chr5_665281.bam'
     tumor = '../data/chr5_665281/tumour_chr5_665281.bam'
     exons = BedReader('../data/chr5_665281/dist.bed').exonList
 
     for exon in exons:
-        tmb = calculate_tmb(tumor, normal, exon)
+        tmb = get_tmb_from_files(tumor, normal, exon)
         assert tmb.somatic_sites == [665281, 665308, 665336]
         assert_almost_equal(tmb.tmb, 50847.46)
 
@@ -23,4 +29,43 @@ def test_panel_indexes():
     exons = BedReader('../data/panel_az_600/panel_az600_chr7_MET.bed').exonList
 
     for exon in exons:
-        tmb = calculate_tmb(tumor, normal, exon)
+        tmb = get_tmb_from_files(tumor, normal, exon)
+        print(tmb)
+
+
+def test_cosine_and_other_dist():
+    exon = Exome('chr1', 0, 2, "")
+    normal = {0: {'T': 77, 'DEL': 19}}
+    tumor = {0: {'T': 23, 'DEL': 10}}
+    somatic_sites = hypotesis_test(normal, tumor)
+    tmb = calculate_tmb_from_sites(exon, somatic_sites)
+
+    assert len(tmb.somatic_sites) == 1
+    assert tmb.somatic_sites == [0]
+    assert tmb.tmb == 500000
+
+
+def test_fill_maps_reads():
+    normal_file = '../data/chr5_665281/normal_chr5_665281.bam'
+    tumor_file = '../data/chr5_665281/tumour_chr5_665281.bam'
+    exons = BedReader('../data/chr5_665281/dist.bed').exonList
+
+    tumor = BamReader(tumor_file).file
+    normal = BamReader(normal_file).file
+
+    for exon in exons:
+        mapNormalR = fill_maps_reads(normal, exon)
+        mapNormalP = fill_maps_pileup(normal, exon)
+        mapTumorR = fill_maps_reads(tumor, exon)
+        mapTumorP = fill_maps_pileup(tumor, exon)
+
+        print("normal reads:", mapNormalR)
+        print("normal pileup:", mapNormalP)
+
+        print("  tumor reads:", mapTumorR)
+        print(" tumor pileup:", mapTumorP)
+
+
+def test_string():
+    str = 'AAGTGCGAGGTGACAGAGTCCCTGCTCTGGGGACACCGGCAGTGCCGCGAGGCATGGAGCGGGGGCTACTTGCCCCCTTGCACCTTCTGGTCGTAGGTGCCTGCGTGCTTGTAGCCGGACACATAGCCAGACTCGTCCACCAGATCCACG'
+    print(str[128])

tmb/bedreader.py

@@ -11,12 +11,13 @@ def __init__(self, path):
     def read_bed(self, path):
         exonlist = []
 
-        file = open(path, "r")
-        for x in file:
-            line = x.split()
-            if len(line) == 4:
-                exome = Exome(str(line[0]), int(line[1]), int(line[2]), str(line[3]))
-            else:
-                exome = Exome(str(line[0]), int(line[1]), int(line[2]), '.')
-            exonlist.append(exome)
+        with open(path, "r") as file:
+            for line in file:
+                if not line == '\n':
+                    line_cols = line.split()
+                    if len(line_cols) == 4:
+                        exome = Exome(str(line_cols[0]), int(line_cols[1]), int(line_cols[2]), str(line_cols[3]))
+                    else:
+                        exome = Exome(str(line_cols[0]), int(line_cols[1]), int(line_cols[2]), '.')
+                    exonlist.append(exome)
         return exonlist

tmb/config.py

@@ -7,6 +7,9 @@ class Config:
     baseq = 25
     mincov = 2
     th = 1
+    pvalue = 10E-60
+    log_level = "WARNING"
+    header = False
 
     @staticmethod
     def set_config(args):
@@ -20,6 +23,12 @@ def set_config(args):
             Config.mincov = args.mincov
         if args.th:
             Config.th = args.th
+        if args.pvalue:
+            Config.pvalue = args.pvalue
+        if args.log:
+            Config.log_level = args.log
+        if args.header:
+            Config.header = True
 
     def __init__(self, args):
         self.normal = args.normal