Python lib.BtLog类代码示例

    def parseFasta(self, fasta_f, fasta_type):
        print BtLog.status_d['1'] % ('FASTA', fasta_f)
        self.assembly_f = abspath(fasta_f)
        if (fasta_type):
            # Set up CovLibObj for coverage in assembly header
            self.covLibs[fasta_type] = CovLibObj(fasta_type, fasta_type, fasta_f)

        for name, seq in BtIO.readFasta(fasta_f):
            blObj = BlObj(name, seq)
            if not blObj.name in self.dict_of_blobs:
                self.seqs += 1
                self.length += blObj.length
                self.n_count += blObj.n_count
                
                if (fasta_type):
                    cov = BtIO.parseCovFromHeader(fasta_type, blObj.name)
                    self.covLibs[fasta_type].cov_sum += cov
                    blObj.addCov(fasta_type, cov)

                self.order_of_blobs.append(blObj.name)
                self.dict_of_blobs[blObj.name] = blObj
            else:
                BtLog.error('5', blObj.name)
        
        if self.seqs == 0 or self.length == 0:
            BtLog.error('1')

def readCas(infile, order_of_blobs):
    seqs_total, reads_total, reads_mapped = checkCas(infile)
    progress_unit = int(len(order_of_blobs)/100)
    cas_line_re = re.compile(r"\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+.\d{2})\s+(\d+)\s+(\d+.\d{2})")
    command = "clc_mapping_info -n " + infile
    cov_dict = {}
    read_cov_dict = {}
    seqs_parsed = 0 
    if (runCmd(command)):
        for line in runCmd(command):
            cas_line_match = cas_line_re.search(line)
            if cas_line_match:
                idx = int(cas_line_match.group(1)) - 1 # -1 because index of contig list starts with zero 
                try:
                    name = order_of_blobs[idx]
                    reads = int(cas_line_match.group(3))
                    cov = float(cas_line_match.group(6))
                    cov_dict[name] = cov
                    read_cov_dict[name] = reads
                    seqs_parsed += 1
                except:
                    pass
            BtLog.progress(seqs_parsed, progress_unit, seqs_total)
        BtLog.progress(seqs_total, progress_unit, seqs_total)
    return cov_dict, reads_total, reads_mapped, read_cov_dict

def readBam(infile, fasta_headers):
    reads_total, reads_mapped = checkBam(infile)
    progress_unit = int(int(reads_total)/1000)
    base_cov_dict = {}
    cigar_match_re = re.compile(r"(\d+)M") # only gets digits before M's

    read_cov_dict = {}
    # execute samtools to get only mapped reads from primary alignment
    command = "samtools view -q " + str(mq) + " -F 256 -F 4 " + infile
    # only one counter since only yields mapped reads
    parsed_reads = 0
    for line in runCmd(command):
        match = line.split("\t")
        seq_name = match[2]
        if seq_name not in fasta_headers:
            print BtLog.warn_d['2'] % (seq_name, infile)
        else:
            read_cov_dict[seq_name] = read_cov_dict.get(seq_name, 0) + 1
            if not (no_base_cov_flag):
                base_cov = sum([int(matching) for matching in cigar_match_re.findall(match[5])])
                if (base_cov):
                    base_cov_dict[seq_name] = base_cov_dict.get(seq_name, 0) + base_cov
            parsed_reads += 1
        BtLog.progress(parsed_reads, progress_unit, reads_total)
    BtLog.progress(reads_total, progress_unit, reads_total)
    return base_cov_dict, read_cov_dict, reads_total, parsed_reads

def readBam(infile, set_of_blobs):
    reads_total, reads_mapped = checkBam(infile)
    progress_unit = int(int(reads_mapped)/1000) + 1 # lazy fix
    base_cov_dict = {}
    read_cov_dict = {}
    cigar_match_re = re.compile(r"(\d+)M") # only gets digits before M's
    # execute samtools to get only mapped reads
    command = "samtools view -F 4 " + infile
    # ADD flag picard -F 1028 to not consider optical duplicates
    #command = "samtools view -F 1028 " + infile
    # only one counter since only yields mapped reads
    parsed_reads = 0 
    for line in runCmd(command):
        match = line.split("\t")
        if match >= 11:
            seq_name = match[2]
            base_cov = sum([int(matching) for matching in cigar_match_re.findall(match[5])])
            if (base_cov):
                parsed_reads += 1
                if seq_name not in set_of_blobs:
                    print BtLog.warn_d['2'] % (seq_name, infile)
                else:
                    base_cov_dict[seq_name] = base_cov_dict.get(seq_name, 0) + base_cov 
                    read_cov_dict[seq_name] = read_cov_dict.get(seq_name, 0) + 1 
        BtLog.progress(parsed_reads, progress_unit, reads_total)
    BtLog.progress(reads_total, progress_unit, reads_total)
    if not int(reads_mapped) == int(parsed_reads):
        print warn_d['3'] % (reads_mapped, parsed_reads)
    return base_cov_dict, reads_total, parsed_reads, read_cov_dict

def parseCatColour(catcolour_f):
    catcolour_dict = {}
    with open(catcolour_f) as fh:
        for l in fh:
            try:
                seq_name, category = l.rstrip("\n").split(",")
                catcolour_dict[seq_name] = category
            except:
                BtLog.error('23', catcolour_f)
    return catcolour_dict

def writeNodesDB(nodesDB, nodesDB_f):
    nodes_count = nodesDB['nodes_count']
    i = 0
    with open(nodesDB_f, 'w') as fh:
        fh.write("# nodes_count = %s\n" % nodes_count) 
        for node in nodesDB:
            if not node == "nodes_count": 
                i += 1
                BtLog.progress(i, 1000, nodes_count)
                fh.write("%s\t%s\t%s\t%s\n" % (node, nodesDB[node]['rank'], nodesDB[node]['name'], nodesDB[node]['parent']))

def parseRefCov(refcov_f):
    refcov_dict = {}
    with open(refcov_f) as fh:
        for l in fh:
            try:
                cov_lib, reads_total_ref, reads_mapped_ref = l.split(",")
                refcov_dict[cov_lib] = {
                                        'reads_total' : int(reads_total_ref), 
                                        'reads_mapped' : int(reads_mapped_ref)
                                       }
            except:
                BtLog.error('21', refcov_f)
    return refcov_dict

def parseCovFile(cov_f):
    cov_dict = {}
    with open(cov_f) as fh:
        for l in fh:
            try:
                seq_name, cov = l.rstrip("\n").split("\t")
                if float(cov) < 0.02:
                    cov_dict[seq_name] = 0.02
                else:
                    cov_dict[seq_name] = float(cov)
            except:
                BtLog.error('25', cov_f)
    return cov_dict

 def computeTaxonomy(self, taxrules, nodesDB):
     tree_lists = BtTax.getTreeList(self.set_of_taxIds, nodesDB)
     self.lineages = BtTax.getLineages(tree_lists, nodesDB)
     self.taxrules = taxrules
     i = 0
     for blObj in self.dict_of_blobs.values():
         i += 1
         BtLog.progress(i, 100, self.seqs)
         for taxrule in taxrules:
             if (blObj.hits):
                 blObj.taxonomy[taxrule] = BtTax.taxRule(taxrule, blObj.hits, self.lineages)
             else:
                 blObj.taxonomy[taxrule] = BtTax.noHit()

def readNodesDB(nodesDB_f):
    nodesDB = {}
    nodes_count = 0
    i = 0
    with open(nodesDB_f) as fh:
        for line in fh:
            if line.startswith("#"):
                nodes_count = int(line.lstrip("# nodes_count = ").rstrip("\n"))
            else:
                i += 1
                node, rank, name, parent = line.rstrip("\n").split("\t")
                nodesDB[node] = {'rank' : rank, 'name' : name, 'parent' : parent}
                BtLog.progress(i, 1000, nodes_count)
    return nodesDB

def parse_labels(labels):
    label_d = {}
    name, groups = '', ''
    if (labels):
        try:
            for label in labels:
                name, groups = str(label).split("=")
                if "," in groups:
                    for group in groups.split(","):
                        label_d[group] = name
                else:
                    label_d[groups] = name
        except:
            BtLog.error('17', labels)
    return label_d

def checkBam(infile):
    print BtLog.status_d['10']
    if not (which('samtools')):
        BtLog.error('7')
    reads_mapped_re = re.compile(r"(\d+)\s\+\s\d+\smapped")
    reads_total_re = re.compile(r"(\d+)\s\+\s\d+\sin total")
    reads_total, reads_mapped = 0, 0
    output = ''
    command = "samtools flagstat " + infile
    for line in runCmd(command):
        output += line
    reads_mapped = int(reads_mapped_re.search(output).group(1))
    reads_total = int(reads_total_re.search(output).group(1))
    print BtLog.status_d['11'] % ('{:,}'.format(reads_mapped), '{:,}'.format(reads_total), '{0:.1%}'.format(reads_mapped/reads_total))
    return reads_total, reads_mapped

def readCov(infile, set_of_blobs):
    old_cov_line_re = re.compile(r"^(\S+)\t(\d+\.*\d*)")
    base_cov_dict = {}

    cov_line_re = re.compile(r"^(\S+)\t(\d+\.*\d*)\t(\d+\.*\d*)")
    reads_total = 0
    reads_mapped = 0
    read_cov_dict = {}

    seqs_parsed = 0
    progress_unit = 1
    old_format = 1
    with open(infile) as fh:
        for line in fh:
            if line.startswith("#"):
                old_format = 0
            if old_format == 0:
                if line.startswith("# Total Reads"):
                    reads_total = int(line.split(" = ")[1])
                elif line.startswith("# Mapped Reads"):
                    reads_mapped = int(line.split(" = ")[1])
                elif line.startswith("# Unmapped Reads"):
                    pass
                elif line.startswith("# Parameters"):
                    pass
                elif line.startswith("# contig_id"):
                    pass
                else:
                    match = cov_line_re.search(line)
                    if match:
                        seqs_parsed += 1
                        name, read_cov, base_cov = match.group(1), int(match.group(2)), float(match.group(3))
                        if name not in set_of_blobs:
                            print BtLog.warn_d['2'] % (name, infile)
                        read_cov_dict[name] = read_cov
                        base_cov_dict[name] = base_cov
            else:
                match = old_cov_line_re.search(line)
                if match:
                    seqs_parsed += 1
                    name, base_cov = match.group(1), float(match.group(2))
                    if name not in set_of_blobs:
                        print BtLog.warn_d['2'] % (name, infile)
                    base_cov_dict[name] = base_cov
            BtLog.progress(seqs_parsed, progress_unit, len(set_of_blobs))
        #BtLog.progress(len(set_of_blobs), progress_unit, len(set_of_blobs))
    return base_cov_dict, reads_total, reads_mapped, read_cov_dict

def checkCas(infile):
    print BtLog.status_d['12']
    if not (which('clc_mapping_info')):
        BtLog.error('20')
    seqs_total_re = re.compile(r"\s+Contigs\s+(\d+)")
    reads_total_re = re.compile(r"\s+Reads\s+(\d+)")
    reads_mapping_re = re.compile(r"\s+Mapped reads\s+(\d+)\s+(\d+.\d+)\s+\%")
    seqs_total, reads_total, reads_mapping, mapping_rate = 0, 0, 0, 0.0
    output = ''
    command = "clc_mapping_info -s " + infile
    for line in runCmd(command):
        output += line
    seqs_total = int(seqs_total_re.search(output).group(1))
    reads_mapped = int(reads_mapping_re.search(output).group(1))
    reads_total = int(reads_total_re.search(output).group(1))
    print BtLog.status_d['11'] % ('{:,}'.format(reads_mapped), '{:,}'.format(reads_total), '{0:.1%}'.format(reads_mapped/reads_total))
    return seqs_total, reads_total, reads_mapped

def readTax(infile, set_of_blobs):
    '''
    If more fields need to be parsed:
        - change hit_line_re
        - catch matches in variables
        - add as key-value pairs to hitDict
    '''
    hit_line_re = re.compile(r"^(\S+)\s+(\d+)[\;?\d+]*\s+(\d+\.*\d*)") # TEST TEST , if not split it afterwards
    with open(infile) as fh:
        for line in fh:
            match = hit_line_re.search(line)
            if match:
                hitDict = {
                    'name' : match.group(1),
                    'taxId' : match.group(2), # string because if int, conversion is a nightmare ...
                    'score' : float(match.group(3))
                    }
                if hitDict['name'] not in set_of_blobs:
                    BtLog.error('19', hitDict['name'], infile)
                if hitDict['taxId'] == 'N/A':
                    BtLog.error('22', infile)
                yield hitDict