#!usr/bin/env python
# -*- coding: utf-8 -*-

'''
This script plots the normalized coverage distribution per amplicon in a
given directory. The average coverage per target region is normalized by
dividing this value by the overall average coverage of each target sequences.
'''

#Load modules
import pybedtools
import re
import sys
import glob
from collections import namedtuple, defaultdict, OrderedDict
import numpy as np
import matplotlib.pyplot as plt
from operator import itemgetter
import seaborn as sns
import pandas as pd

#INTERVALS_BED = '/home/benflies/NGS_Data/TST_15-B-manifest.bed'

INTERVALS_BED = '/home/benflies/NGS_Data/00100-1407755742_Regions.bed'


COVERAGE_THRESHOLD = 1000

directory = '/home/benflies/NGS_Data/BAM/bam_hpx_unique'

fig,ax1 = plt.subplots()
plt.hold = True
boxes = []

bam_file_list = [f for f in glob.iglob(directory+"/*.bam")]

collected = defaultdict(list)

for file in bam_file_list:

    print file

    almnt = pybedtools.BedTool(file)

    IntervalColumns = namedtuple('bed', ['chr', 'start', 'end', 'gene'])
    intervals_list = []

    if INTERVALS_BED:
        with open(INTERVALS_BED, "r") as fin:
            for line in fin.readlines():
                line = line.rstrip('\n')
                line = re.split(r'\t+', line.rstrip('\t'))

                # For Agilent Haloplex BED files
                if len(line) == 4:
                    line[1] = int(line[1])
                    line[2] = int(line[2])
                    bed_line = IntervalColumns(*(line[0], int(line[1]), int(line[2]), line[3]))
                    intervals_list.append(bed_line)

                # For Illumina Trusight Tumor 15 BED files
                if len(line) == 12:
                    line[1] = int(line[1])
                    line[2] = int(line[2])
                    bed_line = IntervalColumns(*(line[0], int(line[1]), int(line[2]), line[3]))
                    intervals_list.append(bed_line)
    else:
        print("ERROR: Provide an interval list (bed format)")

    coverage_result = almnt.coverage(intervals_list)
    print coverage_result
