scripts/release/rethreshold_family.py
def generate_search_stats(family_dir, scores_file='species', tag_miRNA=True):
    """
    Function to generate useful search stats per family

    family_dir: A valid Rfam family checkout directory where pre-computed searches
    were ran
    scores_file: A string specifying the scores file to parse (outlist, species)

    return: report string
    """

    rfam_acc = os.path.basename(family_dir)

    # check point flags
    flag_curr = 0
    flag_rev = 0
    elements = None
    prev_line = None

    seen_ga = False
    seen_rev_before_ga = False
    ga_bit_score = 0.0
    rev_bit_score = 0.0
    ga_rev_seq_gap = 0  # gap in sequences between GA/REV thresholds
    is_miRNA = 0

    seed_above_ga = None
    last_seed_seen = None
    seq_position = 0
    last_seed_pos = 0
    seed_above_ga_pos = 0
    ga_position = 0
    rev_position = 0
    position = 1  # keeps hold of the index position with respect to the entire outlist
    review_family = False
    full_check = False

    unique_seeds = {}

    # initialization of counts
    counts = {"seed_above_ga": 0,
              "full_above_ga": 0,
              "full_below_ga": 0,
              "seed_below_ga": 0,
              "seed_below_rev": 0,
              "full_below_rev": 0,
              "other_below_ga": 0}

    # fetch miRNA accessions from the database
    miRNAs = {}
    if tag_miRNA is True:
        miRNAs.update(db.fetch_type_specific_rfam_accessions("miRNA", return_type="dict"))

    if rfam_acc in miRNAs:
        is_miRNA = 1

    # get some useful numbers from the database
    num_seed_seqs_db = db.get_number_of_seed_sequences(rfam_acc)
    num_full_hits_db = db.get_number_of_full_hits(rfam_acc)
    unique_ncbi_ids_db = db.get_family_unique_ncbi_ids(rfam_acc)

    seedoutlist = os.path.join(family_dir, "seedoutlist")
    missing_seed_seqs_so = get_missing_seeds_seedoutlist(num_seed_seqs_db, seedoutlist)

    scores_fp = open(os.path.join(family_dir, scores_file), 'r')

    # this will basically read the first line which is a header so no harm
    line = scores_fp.readline()
    prev_line = line

    ncbi_ids_from_hits = set()

    # generate stats
    for line in scores_fp:
        position += 1  # starts from 1 because we read the 1st line out of the loop

        # make flag_curr = 1 when we reach that line
        if line.find("CURRENT GA THRESHOLD") != -1:
            flag_curr = 1
            # if we reached this point, it means we saw GA
            seen_ga = True
            ga_position = position

            # get all the elements of the last score line above the GA threshold
            seed_above_ga = last_seed_seen
            seed_above_ga_pos = last_seed_pos

            # get GA threshold
            elements = line.split(' ')
            ga_bit_score = float(elements[-3])
            continue

        # when we reach the reversed sequence line set the flag to 1
        if line.find("BEST REVERSED") != -1:
            flag_rev = 1
            rev_position = position
            # check if GA is false at this point. If yes, this means we saw REV first.
            # setting flag to True
            if seen_ga is False:
                seen_rev_before_ga = True
                continue

        if line[0] != '#':
            # increase sequence position
            # seq_position += 1
            elements = [x for x in line.strip().split(' ') if x != '']

            # first line after hitting REV line
            if flag_rev == 1 and rev_bit_score == 0.0:
                rev_bit_score = float(elements[0])

            # add id to ncbi_ids
            ncbi_ids_from_hits.add(elements[5])

            # we are above the GA
            if flag_curr == 0 and flag_rev == 0:
                if elements[2] == "SEED":
                    # make sure the sequences is not in the dictionary and that it starts from 1
                    if elements[3] not in unique_seeds:
                        counts["seed_above_ga"] += 1
                        unique_seeds[elements[3]] = (elements[8], elements[9])

                elif elements[2] == "FULL" or elements[2] == "FULL-SEED":
                    counts["full_above_ga"] += 1

            # we are somewhere in between current threshold and reversed cutoff
            elif flag_curr == 1 and flag_rev == 0:
                if elements[2] == "SEED":
                    if elements[3] not in unique_seeds:
                        counts["seed_below_ga"] += 1
                        unique_seeds[elements[3]] = (elements[8], elements[9])

            elif flag_curr == 1 and flag_rev == 1:
                if elements[2] == "SEED":
                    if elements[3] not in unique_seeds:
                        counts["seed_below_rev"] += 1
                        unique_seeds[elements[3]] = (elements[8], elements[9])

            # if between GA and REV count sequences
            if ((flag_curr == 1 and flag_rev == 0) or (flag_curr == 0 and flag_rev == 1)):
                ga_rev_seq_gap += 1

            # always stores the last seed up to the current iteration point
            # at the end of the file, last_seed_seen will be holding the last SEED sequence
            # seen in the outlist file
            if elements[2] == "SEED":
                last_seed_seen = elements
                # sets position
                last_seed_pos = position

        # current line becomes previous at the end of each iteration
        prev_line = line

    scores_fp.close()

    # computes the number of any missing SEED sequences. That is SEEDs that do not appear in the outlist
    missing_seed_seqs_o = abs(
        num_seed_seqs_db - (counts["seed_above_ga"] + counts["seed_below_ga"] + counts["seed_below_rev"]))

    # compute the total number of ncbi_ids including
    # total_ncbi_ids_found = len(list(set(unique_ncbi_ids_db).union(ncbi_ids_from_hits)))

    # calulates the number of new ncbi ids added to the full region after a new search
    # new_ncbi_ids_found = abs(total_ncbi_ids_found - len(unique_ncbi_ids_db))

    # ABS(NFULL_OLD-NFULL_NEW) > 0.1 * NFULL_OLD
    # full_diff = abs(num_full_hits_db - (counts["full_above_ga"] + counts["full_below_ga"]))

    # compute GA/REV bit score difference
    ga_rev_bitscore_diff = abs(ga_bit_score - rev_bit_score)

    # if full_diff > (0.1 * num_full_hits_db):
    # full_check = True

    # constraints to be met for reviewing families
    if (seen_rev_before_ga or (counts["seed_below_ga"] > 0) or (counts["seed_below_rev"] > 0)):
        review_family = True

    fields = [rfam_acc, str(num_seed_seqs_db), str(counts["seed_above_ga"]), str(counts["seed_below_ga"]),
              str(counts["seed_below_rev"]), str(missing_seed_seqs_o), str(missing_seed_seqs_so),
              str(ga_bit_score), str(rev_bit_score), str(ga_rev_bitscore_diff), str(ga_rev_seq_gap),
              str(int(seen_rev_before_ga)), seed_above_ga[0], str(seed_above_ga_pos),
              str(ga_position), str(rev_position), last_seed_seen[0], str(last_seed_pos),
              str(int(review_family))]

    """
    fields = [rfam_acc, str(num_seed_seqs_db), str(counts["seed_above_ga"]), str(counts["seed_below_ga"]),
    str(counts["seed_below_rev"]), str(missing_seed_seqs_o), str(missing_seed_seqs_so),
    str(num_full_hits_db), str(counts["full_above_ga"]), str(len(unique_ncbi_ids_db)),
    str(new_ncbi_ids_found), str(ga_bit_score), str(rev_bit_score), str(ga_rev_bitscore_diff),
    str(ga_rev_seq_gap), str(int(seen_rev_before_ga)), seed_above_ga[0], str(seed_above_ga_pos),
    str(ga_position), str(rev_position), last_seed_seen[0], str(last_seed_pos), str(int(review_family))]
    """

    if tag_miRNA is True:
        fields.append(str(is_miRNA))

    print ('\t'.join(fields))

    Function to generate useful search stats per family

    family_dir: A valid Rfam family checkout directory where pre-computed searches
    were ran
    scores_file: A string specifying the scores file to parse (outlist, species)

    return: report string
    
Function to generate useful search stats per family
family_dir: A valid Rfam family checkout directory where pre-computed searches
were ran
scores_file: A string specifying the scores file to parse (outlist, species)

report string
Function to generate useful search stats per family family_dir : A valid Rfam family checkout directory where pre - computed searches were ran scores_file : A string specifying the scores file to parse ( outlist species ) report string
def generate_search_stats(family_dir, scores_file='species', tag_miRNA=True):
    rfam_acc = os.path.basename(family_dir)
    flag_curr = 0
    flag_rev = 0
    elements = None
    prev_line = None
    seen_ga = False
    seen_rev_before_ga = False
    ga_bit_score = 0.0
    rev_bit_score = 0.0
    ga_rev_seq_gap = 0  
    is_miRNA = 0
    seed_above_ga = None
    last_seed_seen = None
    seq_position = 0
    last_seed_pos = 0
    seed_above_ga_pos = 0
    ga_position = 0
    rev_position = 0
    position = 1  
    review_family = False
    full_check = False
    unique_seeds = {}
    counts = {"seed_above_ga": 0,
              "full_above_ga": 0,
              "full_below_ga": 0,
              "seed_below_ga": 0,
              "seed_below_rev": 0,
              "full_below_rev": 0,
              "other_below_ga": 0}
    miRNAs = {}
    if tag_miRNA is True:
        miRNAs.update(db.fetch_type_specific_rfam_accessions("miRNA", return_type="dict"))
    if rfam_acc in miRNAs:
        is_miRNA = 1
    num_seed_seqs_db = db.get_number_of_seed_sequences(rfam_acc)
    num_full_hits_db = db.get_number_of_full_hits(rfam_acc)
    unique_ncbi_ids_db = db.get_family_unique_ncbi_ids(rfam_acc)
    seedoutlist = os.path.join(family_dir, "seedoutlist")
    missing_seed_seqs_so = get_missing_seeds_seedoutlist(num_seed_seqs_db, seedoutlist)
    scores_fp = open(os.path.join(family_dir, scores_file), 'r')
    line = scores_fp.readline()
    prev_line = line
    ncbi_ids_from_hits = set()
    for line in scores_fp:
        position += 1  
        if line.find("CURRENT GA THRESHOLD") != -1:
            flag_curr = 1
            seen_ga = True
            ga_position = position
            seed_above_ga = last_seed_seen
            seed_above_ga_pos = last_seed_pos
            elements = line.split(' ')
            ga_bit_score = float(elements[-3])
            continue
        if line.find("BEST REVERSED") != -1:
            flag_rev = 1
            rev_position = position
            if seen_ga is False:
                seen_rev_before_ga = True
                continue
        if line[0] != '#':
            elements = [x for x in line.strip().split(' ') if x != '']
            if flag_rev == 1 and rev_bit_score == 0.0:
                rev_bit_score = float(elements[0])
            ncbi_ids_from_hits.add(elements[5])
            if flag_curr == 0 and flag_rev == 0:
                if elements[2] == "SEED":
                    if elements[3] not in unique_seeds:
                        counts["seed_above_ga"] += 1
                        unique_seeds[elements[3]] = (elements[8], elements[9])
                elif elements[2] == "FULL" or elements[2] == "FULL-SEED":
                    counts["full_above_ga"] += 1
            elif flag_curr == 1 and flag_rev == 0:
                if elements[2] == "SEED":
                    if elements[3] not in unique_seeds:
                        counts["seed_below_ga"] += 1
                        unique_seeds[elements[3]] = (elements[8], elements[9])
            elif flag_curr == 1 and flag_rev == 1:
                if elements[2] == "SEED":
                    if elements[3] not in unique_seeds:
                        counts["seed_below_rev"] += 1
                        unique_seeds[elements[3]] = (elements[8], elements[9])
            if ((flag_curr == 1 and flag_rev == 0) or (flag_curr == 0 and flag_rev == 1)):
                ga_rev_seq_gap += 1
            if elements[2] == "SEED":
                last_seed_seen = elements
                last_seed_pos = position
        prev_line = line
    scores_fp.close()
    missing_seed_seqs_o = abs(
        num_seed_seqs_db - (counts["seed_above_ga"] + counts["seed_below_ga"] + counts["seed_below_rev"]))
    ga_rev_bitscore_diff = abs(ga_bit_score - rev_bit_score)
    if (seen_rev_before_ga or (counts["seed_below_ga"] > 0) or (counts["seed_below_rev"] > 0)):
        review_family = True
    fields = [rfam_acc, str(num_seed_seqs_db), str(counts["seed_above_ga"]), str(counts["seed_below_ga"]),
              str(counts["seed_below_rev"]), str(missing_seed_seqs_o), str(missing_seed_seqs_so),
              str(ga_bit_score), str(rev_bit_score), str(ga_rev_bitscore_diff), str(ga_rev_seq_gap),
              str(int(seen_rev_before_ga)), seed_above_ga[0], str(seed_above_ga_pos),
              str(ga_position), str(rev_position), last_seed_seen[0], str(last_seed_pos),
              str(int(review_family))]
    if tag_miRNA is True:
        fields.append(str(is_miRNA))
    print ('\t'.join(fields))
def generate_search_stats ( family_dir , scores_file = 'species' , tag_miRNA = True ) : rfam_acc = os . path . basename ( family_dir ) flag_curr = 0 flag_rev = 0 elements = None prev_line = None seen_ga = False seen_rev_before_ga = False ga_bit_score = 0.0 rev_bit_score = 0.0 ga_rev_seq_gap = 0 is_miRNA = 0 seed_above_ga = None last_seed_seen = None seq_position = 0 last_seed_pos = 0 seed_above_ga_pos = 0 ga_position = 0 rev_position = 0 position = 1 review_family = False full_check = False unique_seeds = { } counts = { "seed_above_ga" : 0 , "full_above_ga" : 0 , "full_below_ga" : 0 , "seed_below_ga" : 0 , "seed_below_rev" : 0 , "full_below_rev" : 0 , "other_below_ga" : 0 } miRNAs = { } if tag_miRNA is True : miRNAs . update ( db . fetch_type_specific_rfam_accessions ( "miRNA" , return_type = "dict" ) ) if rfam_acc in miRNAs : is_miRNA = 1 num_seed_seqs_db = db . get_number_of_seed_sequences ( rfam_acc ) num_full_hits_db = db . get_number_of_full_hits ( rfam_acc ) unique_ncbi_ids_db = db . get_family_unique_ncbi_ids ( rfam_acc ) seedoutlist = os . path . join ( family_dir , "seedoutlist" ) missing_seed_seqs_so = get_missing_seeds_seedoutlist ( num_seed_seqs_db , seedoutlist ) scores_fp = open ( os . path . join ( family_dir , scores_file ) , 'r' ) line = scores_fp . readline ( ) prev_line = line ncbi_ids_from_hits = set ( ) for line in scores_fp : position += 1 if line . find ( "CURRENT GA THRESHOLD" ) != - 1 : flag_curr = 1 seen_ga = True ga_position = position seed_above_ga = last_seed_seen seed_above_ga_pos = last_seed_pos elements = line . split ( ' ' ) ga_bit_score = float ( elements [ - 3 ] ) continue if line . find ( "BEST REVERSED" ) != - 1 : flag_rev = 1 rev_position = position if seen_ga is False : seen_rev_before_ga = True continue if line [ 0 ] != '#' : elements = [ x for x in line . strip ( ) . split ( ' ' ) if x != '' ] if flag_rev == 1 and rev_bit_score == 0.0 : rev_bit_score = float ( elements [ 0 ] ) ncbi_ids_from_hits . add ( elements [ 5 ] ) if flag_curr == 0 and flag_rev == 0 : if elements [ 2 ] == "SEED" : if elements [ 3 ] not in unique_seeds : counts [ "seed_above_ga" ] += 1 unique_seeds [ elements [ 3 ] ] = ( elements [ 8 ] , elements [ 9 ] ) elif elements [ 2 ] == "FULL" or elements [ 2 ] == "FULL-SEED" : counts [ "full_above_ga" ] += 1 elif flag_curr == 1 and flag_rev == 0 : if elements [ 2 ] == "SEED" : if elements [ 3 ] not in unique_seeds : counts [ "seed_below_ga" ] += 1 unique_seeds [ elements [ 3 ] ] = ( elements [ 8 ] , elements [ 9 ] ) elif flag_curr == 1 and flag_rev == 1 : if elements [ 2 ] == "SEED" : if elements [ 3 ] not in unique_seeds : counts [ "seed_below_rev" ] += 1 unique_seeds [ elements [ 3 ] ] = ( elements [ 8 ] , elements [ 9 ] ) if ( ( flag_curr == 1 and flag_rev == 0 ) or ( flag_curr == 0 and flag_rev == 1 ) ) : ga_rev_seq_gap += 1 if elements [ 2 ] == "SEED" : last_seed_seen = elements last_seed_pos = position prev_line = line scores_fp . close ( ) missing_seed_seqs_o = abs ( num_seed_seqs_db - ( counts [ "seed_above_ga" ] + counts [ "seed_below_ga" ] + counts [ "seed_below_rev" ] ) ) ga_rev_bitscore_diff = abs ( ga_bit_score - rev_bit_score ) if ( seen_rev_before_ga or ( counts [ "seed_below_ga" ] > 0 ) or ( counts [ "seed_below_rev" ] > 0 ) ) : review_family = True fields = [ rfam_acc , str ( num_seed_seqs_db ) , str ( counts [ "seed_above_ga" ] ) , str ( counts [ "seed_below_ga" ] ) , str ( counts [ "seed_below_rev" ] ) , str ( missing_seed_seqs_o ) , str ( missing_seed_seqs_so ) , str ( ga_bit_score ) , str ( rev_bit_score ) , str ( ga_rev_bitscore_diff ) , str ( ga_rev_seq_gap ) , str ( int ( seen_rev_before_ga ) ) , seed_above_ga [ 0 ] , str ( seed_above_ga_pos ) , str ( ga_position ) , str ( rev_position ) , last_seed_seen [ 0 ] , str ( last_seed_pos ) , str ( int ( review_family ) ) ] """
    fields = [rfam_acc, str(num_seed_seqs_db), str(counts["seed_above_ga"]), str(counts["seed_below_ga"]),
    str(counts["seed_below_rev"]), str(missing_seed_seqs_o), str(missing_seed_seqs_so),
    str(num_full_hits_db), str(counts["full_above_ga"]), str(len(unique_ncbi_ids_db)),
    str(new_ncbi_ids_found), str(ga_bit_score), str(rev_bit_score), str(ga_rev_bitscore_diff),
    str(ga_rev_seq_gap), str(int(seen_rev_before_ga)), seed_above_ga[0], str(seed_above_ga_pos),
    str(ga_position), str(rev_position), last_seed_seen[0], str(last_seed_pos), str(int(review_family))]
    """ if tag_miRNA is True : fields . append ( str ( is_miRNA ) ) print ( '\t' . join ( fields ) )
Function to generate useful search stats per family
family_dir: A valid Rfam family checkout directory where pre-computed searches
were ran
scores_file: A string specifying the scores file to parse (outlist, species)
Function to generate useful search stats per family family_dir : A valid Rfam family checkout directory where pre - computed searches were ran scores_file : A string specifying the scores file to parse ( outlist species )
"""
    Function to generate useful search stats per family

    family_dir: A valid Rfam family checkout directory where pre-computed searches
    were ran
    scores_file: A string specifying the scores file to parse (outlist, species)

    return: report string
    """ # check point flags # gap in sequences between GA/REV thresholds # keeps hold of the index position with respect to the entire outlist # initialization of counts # fetch miRNA accessions from the database # get some useful numbers from the database # this will basically read the first line which is a header so no harm # generate stats # starts from 1 because we read the 1st line out of the loop # make flag_curr = 1 when we reach that line # if we reached this point, it means we saw GA # get all the elements of the last score line above the GA threshold # get GA threshold # when we reach the reversed sequence line set the flag to 1 # check if GA is false at this point. If yes, this means we saw REV first. # setting flag to True # increase sequence position # seq_position += 1 # first line after hitting REV line # add id to ncbi_ids # we are above the GA # make sure the sequences is not in the dictionary and that it starts from 1 # we are somewhere in between current threshold and reversed cutoff # if between GA and REV count sequences # always stores the last seed up to the current iteration point # at the end of the file, last_seed_seen will be holding the last SEED sequence # seen in the outlist file # sets position # current line becomes previous at the end of each iteration # computes the number of any missing SEED sequences. That is SEEDs that do not appear in the outlist # compute the total number of ncbi_ids including # total_ncbi_ids_found = len(list(set(unique_ncbi_ids_db).union(ncbi_ids_from_hits))) # calulates the number of new ncbi ids added to the full region after a new search # new_ncbi_ids_found = abs(total_ncbi_ids_found - len(unique_ncbi_ids_db)) # ABS(NFULL_OLD-NFULL_NEW) > 0.1 * NFULL_OLD # full_diff = abs(num_full_hits_db - (counts["full_above_ga"] + counts["full_below_ga"])) # compute GA/REV bit score difference # if full_diff > (0.1 * num_full_hits_db): # full_check = True # constraints to be met for reviewing families """
    fields = [rfam_acc, str(num_seed_seqs_db), str(counts["seed_above_ga"]), str(counts["seed_below_ga"]),
    str(counts["seed_below_rev"]), str(missing_seed_seqs_o), str(missing_seed_seqs_so),
    str(num_full_hits_db), str(counts["full_above_ga"]), str(len(unique_ncbi_ids_db)),
    str(new_ncbi_ids_found), str(ga_bit_score), str(rev_bit_score), str(ga_rev_bitscore_diff),
    str(ga_rev_seq_gap), str(int(seen_rev_before_ga)), seed_above_ga[0], str(seed_above_ga_pos),
    str(ga_position), str(rev_position), last_seed_seen[0], str(last_seed_pos), str(int(review_family))]
    """
family_dir
scores_file
tag_miRNA