<commit_before>#!/usr/bin/env python

import pysam
import sys
import optparse

def option_parsing():
    parser = optparse.OptionParser(description='Count miRNAs')
    parser.add_option('-p', '--pre_mirs', action='store_true', dest='pre_mirs', help='Count pre-miRNAs', default=False)
    parser.add_option('-m', '--mirs', action='store_true', dest='mirs', help='Count mature miRNAs', default=False)
    inputs = optparse.OptionGroup(parser, 'Inputs')
    inputs.add_option('--alignment', type='string', dest='alignment_file', help='Alignment tabular or sam file')
    inputs.add_option('--gff', type='string', dest='gff_file', help='GFF3 describing both pre-miRNAs and mature miRNAs', metavar='FILE')
    inputs.add_option('--quality_threshold', type='int', dest='quality_threshold', help='Quality of the base alignment to take in to consideration when counting coverage (default=10)', default=10)
    parser.add_option_group(inputs)
    outputs = optparse.OptionGroup(parser, 'Outputs')
    outputs.add_option('--pre_mirs_output', type='string', dest='output_pre_mirs', help='Output file containing table containing the number of hits per pre-miRNA', metavar='FILE', default='output_pre_mirs_count.tab')
    outputs.add_option('--mature_mirs_output', type='string', dest='output_mature_mirs', help='Output file containing the number of hits per mature miRNA', metavar='FILE', default='output_mirs_count.tab')
    outputs.add_option('--lattice', type='string', dest='lattice', help="Output file for the lattice dataframe (if none given it won't be outputed)")
    parser.add_option_group(outputs)
    (options, args) = parser.parse_args()
    """ Check if the options were correctly passed """
    if len(args) > 0:
        parser.error('Please use the options flags to pass your arguments')
    if (not options.alignment_file or not options.gff_file):
        parser.error("Missing file. Both '--alignment' and '--gff' files are needed")
    return options

def get_pre_mir_counts(bamfile, quality_th):
    """
    Takes a AlignmentFile object and and returns a count of reads
    hiting 'miRNA_primary_transcript' positions and the coverage of said region
    """
    count = dict()
    reference_lengths = bamfile.lengths
    it = 0
    """
    Need reference legths beacuse when no end parameter is given to count_coverage an error raises
    This works beacuse : 
    lengths
    tuple of the lengths of the reference sequences.
    This is a read-only attribute. The lengths are in the same order as pysam.AlignmentFile.references
    (from : http://pysam.readthedocs.io/en/latest/api.html)
    """
    for ref in bamfile.references:
        count[ref] = [bamfile.count(ref),
	              bamfile.count_coverage(reference=ref,start=0,end=reference_lengths[it], quality_threshold=quality_th)]
        it += 1
    return count

def get_mir_counts(bamfile, gff_file, quality_th):
    """
    Takes a AlignmentFile and a gff file and computes for
    each 'miRNA' region of the gff the number of reads that hit it
    and the coverage of each of its positions
    and returns a dictionnary : dict[reference name] = [polarity, coverage as [[],[],[],[]]]
    """
    counts = dict()
    try:
        gff = open(gff_file, 'r')
        for line in gff.readlines():
            if line[0] != '#':
                gff_fields = line[:-1].split("\t")
                if gff_fields[2] == 'miRNA':
                    mir_name = gff_fields[0]
                    mir_start = int(gff_fields[3])
                    mir_end = int(gff_fields[4])
                    counts[mir_name] = [bamfile.count(reference=mir_name, start=mir_start, end=mir_end),
                                        bamfile.count_coverage(reference=mir_name, start=mir_start,
                                                               end=mir_end, quality_threshold=quality_th)]
        gff.close()
    except IOError as e:
        sys.stderr.write("Error while reading file %s\n" % gff_file)
        sys.stderr.write("I/O error(%s): %s\n" % (e.errno, e.strerror))
    return counts

def write_dataframe(mirs, outfile, sample):
    """
    Takes a dictionnary dict[reference name] = [count,
                                                [[Hits of each position for 'A'],
                                                 [Hits of each position for 'C'],
                                                 [Hits of each position for 'G'],
                                                 [Hits of each position for 'T']]]
    And prints a dataframe with columns: sample, mir, offset, offsetNorm, counts, countsNorm
    in the outfile
    """
    dataframe = []
    dataframe.append("sample\tmir\toffset\toffsetNorm\tcounts\tcountsNorm")
    for ref in mirs.keys():
        """ For each reference name in mirs write the coverage of each of its positions """
        maximum = 0
        coverage_arrays = mirs[ref][1]
        reference_length = len(coverage_arrays[1])
        for i in range(4):
            if max(coverage_arrays[i]):
                maximum = max(coverage_arrays[i])<commit_msg>Added verification on mature mirs counting to resolve problem of searching for a hit on a mir read on gff but not hit in alignment<commit_after>#!/usr/bin/env python

import pysam
import sys
import optparse

def option_parsing():
    parser = optparse.OptionParser(description='Count miRNAs')
    parser.add_option('-p', '--pre_mirs', action='store_true', dest='pre_mirs', help='Count pre-miRNAs', default=False)
    parser.add_option('-m', '--mirs', action='store_true', dest='mirs', help='Count mature miRNAs', default=False)
    inputs = optparse.OptionGroup(parser, 'Inputs')
    inputs.add_option('--alignment', type='string', dest='alignment_file', help='Alignment tabular or sam file')
    inputs.add_option('--gff', type='string', dest='gff_file', help='GFF3 describing both pre-miRNAs and mature miRNAs', metavar='FILE')
    inputs.add_option('--quality_threshold', type='int', dest='quality_threshold', help='Quality of the base alignment to take in to consideration when counting coverage (default=10)', default=10)
    parser.add_option_group(inputs)
    outputs = optparse.OptionGroup(parser, 'Outputs')
    outputs.add_option('--pre_mirs_output', type='string', dest='output_pre_mirs', help='Output file containing table containing the number of hits per pre-miRNA', metavar='FILE', default='output_pre_mirs_count.tab')
    outputs.add_option('--mature_mirs_output', type='string', dest='output_mature_mirs', help='Output file containing the number of hits per mature miRNA', metavar='FILE', default='output_mirs_count.tab')
    outputs.add_option('--lattice', type='string', dest='lattice', help="Output file for the lattice dataframe (if none given it won't be outputed)")
    parser.add_option_group(outputs)
    (options, args) = parser.parse_args()
    """ Check if the options were correctly passed """
    if len(args) > 0:
        parser.error('Please use the options flags to pass your arguments')
    if (not options.alignment_file or not options.gff_file):
        parser.error("Missing file. Both '--alignment' and '--gff' files are needed")
    return options

def get_pre_mir_counts(bamfile, quality_th):
    """
    Takes a AlignmentFile object and and returns a count of reads
    hiting 'miRNA_primary_transcript' positions and the coverage of said region
    """
    count = dict()
    reference_lengths = bamfile.lengths
    it = 0
    """
    Need reference legths beacuse when no end parameter is given to count_coverage an error raises
    This works beacuse : 
    lengths
    tuple of the lengths of the reference sequences.
    This is a read-only attribute. The lengths are in the same order as pysam.AlignmentFile.references
    (from : http://pysam.readthedocs.io/en/latest/api.html)
    """
    for ref in bamfile.references:
        count[ref] = [bamfile.count(ref),
	              bamfile.count_coverage(reference=ref,start=0,end=reference_lengths[it], quality_threshold=quality_th)]
        it += 1
    return count

def get_mir_counts(bamfile, gff_file, quality_th):
    """
    Takes a AlignmentFile and a gff file and computes for
    each 'miRNA' region of the gff the number of reads that hit it
    and the coverage of each of its positions
    and returns a dictionnary : dict[reference name] = [polarity, coverage as [[],[],[],[]]]
    """
    counts = dict()
    refs = list()
    for h in header['SQ']:
        refs.append(h['SN'])
    try:
        gff = open(gff_file, 'r')
        for line in gff.readlines():
            if line[0] != '#':
                gff_fields = line[:-1].split("\t")
                if gff_fields[2] == 'miRNA':
                    mir_name = gff_fields[0]
                    mir_start = int(gff_fields[3])
                    mir_end = int(gff_fields[4])
                    if mir_name in refs:
                        counts[mir_name] = [bamfile.count(reference=mir_name, start=mir_start, end=mir_end),
                                            bamfile.count_coverage(reference=mir_name, start=mir_start,
                                                                   end=mir_end, quality_threshold=quality_th)]
        gff.close()
    except IOError as e:
        sys.stderr.write("Error while reading file %s\n" % gff_file)
        sys.stderr.write("I/O error(%s): %s\n" % (e.errno, e.strerror))
    return counts

def write_dataframe(mirs, outfile, sample):
    """
    Takes a dictionnary dict[reference name] = [count,
                                                [[Hits of each position for 'A'],
                                                 [Hits of each position for 'C'],
                                                 [Hits of each position for 'G'],
                                                 [Hits of each position for 'T']]]
    And prints a dataframe with columns: sample, mir, offset, offsetNorm, counts, countsNorm
    in the outfile
    """
    dataframe = []
    dataframe.append("sample\tmir\toffset\toffsetNorm\tcounts\tcountsNorm")
    for ref in mirs.keys():
        """ For each reference name in mirs write the coverage of each of its positions """
        maximum = 0
        coverage_arrays = mirs[ref][1]
        reference_length = len(coverage_arrays[1])
        for i in range(4):
            if max(coverage_arrays[i]):
                maximum = max(coverage_arrays[i])
{'fraction_non_alphanumeric': 0.11263259402121505, 'fraction_numerical': 0.0036644165863066536, 'mean_word_length': 2.7278936017253774, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7632502', 'n_tokens_mistral': 3063, 'n_tokens_neox': 2968, 'n_words': 872}