analyzebed.py

# Copyright (C) 2012 INRA
# 
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

from ng6.analysis import Analysis
from weaver.function import ShellFunction
from weaver.abstraction import Map
from jflow import utils
from collections import Counter
import os.path
import numpy
import itertools
import subprocess
import re

class Analyzebed (Analysis):
    
    def define_parameters(self, files, p, li, i_list, options):
        self.add_input_file_list( "input_files", "The files to be archived.", default=files, required=True )
        self.add_output_file_list( "output_links", "Link to the files.", pattern='link2_{basename}' , items=self.input_files )
        self.add_output_file( "software_merge", "TODO", filename="software_merge.bed")
        self.add_parameter("overlap_percent", "Minimum threshold of overlap percent between 2 SV", type=float, default=p)
        self.add_parameter("li_range", "Range of localisation interval around breakpoints", type=int, default=li)
        self.add_parameter("param_sum", "Summary of used parameters", type=str, default=options)
        self.add_parameter_list("indiv_list", "String list of studied individuals name", default=i_list, type=str, required=True)
    
    def get_version(self):
        pindel_help = subprocess.Popen(["pindel","-h"], universal_newlines=True, stdout=subprocess.PIPE)
        pindel_version = re.search("Pindel version ([0-9\.a-b]+),", pindel_help.stdout.read())
        cnvnator_help = subprocess.Popen(["cnvnator"], universal_newlines=True, stderr=subprocess.PIPE)
        cnvnator_version = re.search("CNVnator v([0-9\.]+)", cnvnator_help.stderr.read())
        delly_help = subprocess.Popen(["delly"], universal_newlines=True, stdout=subprocess.PIPE)
        delly_version = re.search("Delly \(Version: ([0-9\.]+)\)", delly_help.stdout.read())
        breakdancer_help = subprocess.Popen(["breakdancer-max"], universal_newlines=True, stderr=subprocess.PIPE)
        breakdancer_version = re.search("breakdancer-max version ([0-9\.]+)", breakdancer_help.stderr.read())
        versions = "Pindel: {0} / CNVnator: {1} / Delly: {2} / Breakdancer: {3}".format(pindel_version.groups()[0],
                     cnvnator_version.groups()[0], delly_version.groups()[0], breakdancer_version.groups()[0])
        return versions
        
    def define_analysis(self):
        self.options = self.param_sum
        self.name = "SV Detection"
        self.description = "A few stats on SV repartition"
        self.software = "Pindel, CNVnator, Delly, Breakdancer"
        
    def post_process(self):
        for i, bedfile in enumerate(self.output_links):
            sample = os.path.splitext(os.path.basename(bedfile))[0]
            soft = self.indiv_list[i]
            #for boxplot
            sizes_list = [] 
            sizes_by_type = {}
            #for accumulation curve
            n_list = []
            n_by_type = {}
            
            #PARSE FILE
            with open(bedfile) as result_file:
                for line in result_file:
                    sv = line.split()
                    size = int(sv[2])-int(sv[1])
                    sv_type = sv[4]
                    members = sv[7].split(";")
                    n = len(members)
                    
                    sizes_list.append(size)
                    n_list.append(n)
                    if sv_type in sizes_by_type: #fill
                        sizes_by_type[sv_type].append(size)
                        n_by_type[sv_type].append(n)
                    else: #if not, create
                        sizes_by_type[sv_type] = [size]
                        n_by_type[sv_type] = [n]

            #GENERATE DB INFOS
            #add "ALL" type
            sizes_by_type["ALL"] = sizes_list
            n_by_type["ALL"] = n_list
            for type in sizes_by_type:
                sizes = numpy.array(sizes_by_type[type])
                boxplot = self.boxplot(sizes)
                accu = self.accumulation(n_by_type[type])
                db_entry = {'boxplot':boxplot[0],
                            'outliers':boxplot[1],
                            'n_tag':accu[0],
                            'n_count':accu[1],
                            'mean':round(numpy.mean(sizes),1),
                            'min':min(sizes),
                            'max':max(sizes),
                            'sd': round(numpy.std(sizes),1),
                            'number':len(sizes),
                            'median':int(numpy.median(sizes))
                            }
                #fill db on soft-specific data
                self._add_result_element(sample, "soft", soft, type)
                for key in db_entry:
                    self._add_result_element(sample, key, db_entry[key], type)
        
        
        #CREATE VENN TABLE
        venn_table = {}
        for key in self.powerset(self.indiv_list):
            if key == ():
                venn_table["all"] = 0
            else:
                venn_table[";".join(key)] = 0
        
        venn_table_by_type = {}
        
        with open(self.software_merge) as merge_list:
            for line in merge_list:
                used_software = []
                svr = line.split()
                svr_type = svr[4]
                
                for column in svr[6].split(";"):
                    #get members
                    used_software.append("".join(column.split("-")[:-1]))
                if svr_type not in venn_table_by_type:  # create the dic
                    for key in self.powerset(self.indiv_list):
                        if key == ():
                            venn_table_by_type[svr_type] = {"all":0}
                        else:
                            venn_table_by_type[svr_type][";".join(key)] = 0
                #update the dic
                venn_table_by_type[svr_type][";".join(self.unique(used_software))] += 1
                venn_table[";".join(self.unique(used_software))] += 1
        #add "ALL" type
        venn_table_by_type['ALL'] = venn_table
        
        #GENERATE DB INFOS
        for type in venn_table_by_type:
            v_table = venn_table_by_type[type]
            #good format
            db_venn = "".join([y for y in str(v_table)[1:-1] if y not in ["\'", " "]])
            #fill the db for inter-soft data
            self._add_result_element("software_merge", "venn_count", db_venn, type)
        
        #output files
        result_files = [self.output_links, self.software_merge]
        self._create_and_archive(self.output_links, None) #TODO pass result_files

    def process(self):
        #create links
        link = ShellFunction(self.get_exec_path("ln") + " -s $1 $2", cmd_format='{EXE} {IN} {OUT}')
        Map(link, inputs=self.input_files, outputs=self.output_links)
        # MERGE INTERSOFT
        mrg = ShellFunction(self.get_exec_path("merge_sv.py") + \
                                    " --output $1 \
                                      --tag $2\
                                      --p $3\
                                      --li $4" \
                                    + " --names" + utils.get_argument_pattern(self.indiv_list, 5)[0] \
                                    + " --files" + utils.get_argument_pattern(self.input_files, len(self.indiv_list)+5)[0],
                                    cmd_format='{EXE} {OUT} {ARG} {IN}')
        mrg(outputs=self.software_merge, inputs=self.input_files, arguments=["soft_merge", self.overlap_percent, self.li_range, self.indiv_list])
        
    def js_prepross(self, my_list):
        """
        facilitate data transfer and comprehension for javascript via jflow & php
        
        :Example: self.js_prepross(['rock', 'scissor', 'paper'])
        >>> "rock;scissor;paper" 
        """
        my_string = ";".join(map(str,my_list))
        return my_string
        
    def accumulation(self, all_n):
        """
        Count how many SV are present in at least 1..N individuals
        
        :return: 1..N & associated count as strings
        
        :Example: self.accumulation([1,1,1,1,2,3,3])
        >>> ["1;2;3", "4;1;2"]
        """
        n_count = []
        all_n = Counter(all_n)
        for i,n in enumerate(all_n): # Rare bug if one n is missing in 1..N --> bug in GUI
            tmp = 0
            for j in list(reversed(range(i+1, len(all_n)+1))):
                tmp = tmp + all_n[j]
            n_count.append(tmp)
        n_tag = self.js_prepross(all_n.keys())
        n_count = self.js_prepross(n_count)
        return([n_tag, n_count])
        
    def boxplot(self, sizes):
        """
        give the values needed for drawing a boxplot
        
        :return: ["lower whisker;lower quartile;median;upper quartile;upper whisker", number of outliers]
        
        :Example: self.boxplot([13799,12699,101,...])
        >>> ["101;509;1099;5349;12299", 4]
        """
        #boxplot
        median = numpy.median(sizes)
        upper_quartile = numpy.percentile(sizes, 75)
        lower_quartile = numpy.percentile(sizes, 25)
        iqr = upper_quartile - lower_quartile
        upper_whisker = sizes[sizes<=upper_quartile+1.5*iqr].max()
        lower_whisker = sizes[sizes>=lower_quartile-1.5*iqr].min()
        boxplot = self.js_prepross([lower_whisker,lower_quartile,median,upper_quartile,upper_whisker])
        #outliers
        outliers = sizes[sizes > upper_whisker]
        numpy.append(outliers, sizes[sizes < lower_whisker])
        outliers = len(outliers)
        return([boxplot, outliers])

    def unique(self, seq, idfun=None):
        """
            unique([1,2,2,3,27,1,2,4]) --> [1,2,3,27,4]
            keep order of appearance intact
        """
        # found on peterbe.com
        if idfun is None:
            def idfun(x): return x
        seen = {}
        result = []
        for item in seq:
            marker = idfun(item)
            if marker in seen: continue
            seen[marker] = 1
            result.append(item)
        return result

    def powerset(self, iterable):
        """ 
        return all possible sub-sets from an iterable
        
        :Example: powerset("abc")
        >>> (''),('a'),('b'),('c'),('a','b'),('a','c'),('b','c'),('a','b','c')
        """
        s = list(iterable)
        return itertools.chain.from_iterable(itertools.combinations(s, r) for r in range(len(s)+1))