scripts/preprocessing/relabel_seed.py
def relabel_seeds_from_rnacentral_urs_mapping(seed, expert_db=None, dest_dir=None, clean=False):
    """
    Relabels the accessions of a SEED alignment using RNAcentral
    identifiers. This is done by matching the seed sequences, with
    sequences existing in RNAcentral using md5 hashing.

    seed: A reformatted seed in Pfam format
    expert_db: An existing RNAcentral expert database
    dest_dir: The path to the destination directory. None by default
    clean: Specifies whether to clean seeds by ommitting sequences not found in RNAcentral

    return: The path to the relabelled SEED alignement
    """

    if dest_dir is None:
        # fetch path of seed alignment
        dest_dir = os.path.split(seed)[0]

    sequence_label = 0
    new_line = ''

    write_log = False

    filename = os.path.split(seed)[1].partition('.')[0]

    new_seed_filename = filename + '_relabelled'

    new_seed_loc = os.path.join(dest_dir, new_seed_filename)
    seed_fp = open(seed, 'r')
    new_seed_fp = open(new_seed_loc, 'w')
    log_fp = None
    sequence_count = 0

    seed_seq_id = ''

    unique_seed_accs = {}
    # dictionary to keep track of sequence mismatches
    sequence_mismatches = {}

    for line in seed_fp:
        # check if this is an actual sequence line
        if line[0] != '#' and len(line) > 1 and line[0:2] != '//':
            line_elements = [x for x in line.strip().split(' ') if x != '']
            sequence_count += 1

            # replace alignment characters
            seed_seq_id = line_elements[0].split('/')[0]
            seed_seq = line_elements[1].replace('.', '')
            seed_seq = seed_seq.replace('-', '').replace('T', 'U').replace('t', 'u').upper()

            rnacentral_id = map_rnacentral_urs_wirh_db_accessions(seed_seq_id, expert_db)

            # if any of the sequences isn't found prepare to log missing accessions
            # sets write_log to only open the file once
            if rnacentral_id is None and write_log is False:
                seed_filename = os.path.basename(seed).partition('.')[0]
                log_fp = open(os.path.join(dest_dir, seed_filename)+'.log', 'w')
                write_log = True

            # write log file if one exists and jump to next iteration
            if rnacentral_id is None:
                log_fp.write("RNACENTRAL ID MISMATCH: %s\n" % seed_seq_id)
                continue

            rnacentral_sequence = fetch_sequence_from_rnacentral(rnacentral_id)
            coordinates = fetch_seed_sequence_coordinates(seed_seq, rnacentral_sequence)

            new_label = ''

            # make sure subsequence was found
            if coordinates[1] != 0:
                new_label = rnacentral_id + '/' + str(coordinates[0]) + '-' + str(coordinates[1])

                if new_label not in unique_seed_accs:
                    unique_seed_accs[new_label] = ''

                # if we found a duplicated id, write log and skip this sequence
                else:
                    if write_log is False:
                        seed_filename = os.path.basename(seed).partition('.')[0]
                        log_fp = open(os.path.join(dest_dir, seed_filename) + '.log', 'w')
                        write_log = True

                    log_fp.write("RNACENTRAL DUPLICATED ID: %s\n" % new_label)
                    continue

            # if we reached this point, it means there was an id match, but no exact sequence match
            else:

                # try mapping SEED sequences in smaller segments
                if write_log is False:
                    seed_filename = os.path.basename(seed).partition('.')[0]
                    log_fp = open(os.path.join(dest_dir, seed_filename) + '.log', 'w')
                    write_log = True

                log_fp.write("RNACENTRAL SEQ MISMATCH: %s\n" % seed_seq_id)
                sequence_mismatches[rnacentral_id] = rnacentral_sequence

                continue
                #else:
                    # here we need to generate the new label with coords and
                #   print (sequence)

            new_line = "\t".join([new_label, line_elements[1], '\n'])

        else:
            new_line = line

        new_seed_fp.write(new_line)

    seed_fp.close()
    new_seed_fp.close()

    sequence_count = len(unique_seed_accs)

    # checks if dictionary isn't empty
    if not bool(sequence_mismatches) is False:

        # write RNAcentral sequences to fasta file to be used with cmalign
        fasta_filename = filename + '_rnac'

        fasta_file = write_fasta_seed_file(sequence_mismatches, fasta_filename, dest_dir)

        if fasta_file is None:
            sys.exit("FILE ERROR: Fasta file %s could not be generated\n" % fasta_filename)

        # generate CM file from original seed
        cmfile = build_temporary_cm_from_seed(seed, dest_dir)

        if cmfile is None:
            sys.exit("FILE ERROR: CM file for seed %s could not be generated\n" % seed_filename)

        # align fasta file to covariance model
        cmaligned_sequences = align_sequences_to_cm(cmfile, fasta_file, dest_dir)

        if cmaligned_sequences is None:
            sys.exit("FILE ERROR: New cmaligned SEED for family %s could not be generated\n" % seed_filename)

        # if the number of sequences in the original SEED is larger than the sequence mismatches
        # this means we have two smaller MSAs that need to be merged
        if sequence_count > len(sequence_mismatches.keys()):
            filename = os.path.split(seed)[1].partition(".")[0]

            cmaligned_relabelled_seed = align_sequences_to_cm(cmfile, new_seed_loc)
            merged_seed = merge_seeds(cmaligned_relabelled_seed, cmaligned_sequences, filename, dest_dir)
            new_seed_loc = merged_seed

        else:
            os.remove(new_seed_loc)
            new_seed_loc = cmaligned_sequences

        # renames files
        os.rename(new_seed_loc, os.path.join(dest_dir, new_seed_filename))
        new_seed_loc = os.path.join(dest_dir, new_seed_filename)
        #final_seed = remove_all_gap_columns(new_seed_loc, filename, dest_dir)

        #if final_seed is None:
        #    sys.exit("FILE ERROR: SEED reformatting failed\n")

    # close log file if one exists
    if write_log is True:
        log_fp.close()

    return new_seed_loc

    Relabels the accessions of a SEED alignment using RNAcentral
    identifiers. This is done by matching the seed sequences, with
    sequences existing in RNAcentral using md5 hashing.

    seed: A reformatted seed in Pfam format
    expert_db: An existing RNAcentral expert database
    dest_dir: The path to the destination directory. None by default
    clean: Specifies whether to clean seeds by ommitting sequences not found in RNAcentral

    return: The path to the relabelled SEED alignement
    
Relabels the accessions of a SEED alignment using RNAcentral
identifiers. This is done by matching the seed sequences, with
sequences existing in RNAcentral using md5 hashing.

A reformatted seed in Pfam format
expert_db: An existing RNAcentral expert database
dest_dir: The path to the destination directory. None by default
clean: Specifies whether to clean seeds by ommitting sequences not found in RNAcentral

The path to the relabelled SEED alignement
Relabels the accessions of a SEED alignment using RNAcentral identifiers . This is done by matching the seed sequences with sequences existing in RNAcentral using md5 hashing . A reformatted seed in Pfam format expert_db : An existing RNAcentral expert database dest_dir : The path to the destination directory . None by default clean : Specifies whether to clean seeds by ommitting sequences not found in RNAcentral The path to the relabelled SEED alignement
def relabel_seeds_from_rnacentral_urs_mapping(seed, expert_db=None, dest_dir=None, clean=False):
    if dest_dir is None:
        dest_dir = os.path.split(seed)[0]
    sequence_label = 0
    new_line = ''
    write_log = False
    filename = os.path.split(seed)[1].partition('.')[0]
    new_seed_filename = filename + '_relabelled'
    new_seed_loc = os.path.join(dest_dir, new_seed_filename)
    seed_fp = open(seed, 'r')
    new_seed_fp = open(new_seed_loc, 'w')
    log_fp = None
    sequence_count = 0
    seed_seq_id = ''
    unique_seed_accs = {}
    sequence_mismatches = {}
    for line in seed_fp:
        if line[0] != '#' and len(line) > 1 and line[0:2] != '//':
            line_elements = [x for x in line.strip().split(' ') if x != '']
            sequence_count += 1
            seed_seq_id = line_elements[0].split('/')[0]
            seed_seq = line_elements[1].replace('.', '')
            seed_seq = seed_seq.replace('-', '').replace('T', 'U').replace('t', 'u').upper()
            rnacentral_id = map_rnacentral_urs_wirh_db_accessions(seed_seq_id, expert_db)
            if rnacentral_id is None and write_log is False:
                seed_filename = os.path.basename(seed).partition('.')[0]
                log_fp = open(os.path.join(dest_dir, seed_filename)+'.log', 'w')
                write_log = True
            if rnacentral_id is None:
                log_fp.write("RNACENTRAL ID MISMATCH: %s\n" % seed_seq_id)
                continue
            rnacentral_sequence = fetch_sequence_from_rnacentral(rnacentral_id)
            coordinates = fetch_seed_sequence_coordinates(seed_seq, rnacentral_sequence)
            new_label = ''
            if coordinates[1] != 0:
                new_label = rnacentral_id + '/' + str(coordinates[0]) + '-' + str(coordinates[1])
                if new_label not in unique_seed_accs:
                    unique_seed_accs[new_label] = ''
                else:
                    if write_log is False:
                        seed_filename = os.path.basename(seed).partition('.')[0]
                        log_fp = open(os.path.join(dest_dir, seed_filename) + '.log', 'w')
                        write_log = True
                    log_fp.write("RNACENTRAL DUPLICATED ID: %s\n" % new_label)
                    continue
            else:
                if write_log is False:
                    seed_filename = os.path.basename(seed).partition('.')[0]
                    log_fp = open(os.path.join(dest_dir, seed_filename) + '.log', 'w')
                    write_log = True
                log_fp.write("RNACENTRAL SEQ MISMATCH: %s\n" % seed_seq_id)
                sequence_mismatches[rnacentral_id] = rnacentral_sequence
                continue
            new_line = "\t".join([new_label, line_elements[1], '\n'])
        else:
            new_line = line
        new_seed_fp.write(new_line)
    seed_fp.close()
    new_seed_fp.close()
    sequence_count = len(unique_seed_accs)
    if not bool(sequence_mismatches) is False:
        fasta_filename = filename + '_rnac'
        fasta_file = write_fasta_seed_file(sequence_mismatches, fasta_filename, dest_dir)
        if fasta_file is None:
            sys.exit("FILE ERROR: Fasta file %s could not be generated\n" % fasta_filename)
        cmfile = build_temporary_cm_from_seed(seed, dest_dir)
        if cmfile is None:
            sys.exit("FILE ERROR: CM file for seed %s could not be generated\n" % seed_filename)
        cmaligned_sequences = align_sequences_to_cm(cmfile, fasta_file, dest_dir)
        if cmaligned_sequences is None:
            sys.exit("FILE ERROR: New cmaligned SEED for family %s could not be generated\n" % seed_filename)
        if sequence_count > len(sequence_mismatches.keys()):
            filename = os.path.split(seed)[1].partition(".")[0]
            cmaligned_relabelled_seed = align_sequences_to_cm(cmfile, new_seed_loc)
            merged_seed = merge_seeds(cmaligned_relabelled_seed, cmaligned_sequences, filename, dest_dir)
            new_seed_loc = merged_seed
        else:
            os.remove(new_seed_loc)
            new_seed_loc = cmaligned_sequences
        os.rename(new_seed_loc, os.path.join(dest_dir, new_seed_filename))
        new_seed_loc = os.path.join(dest_dir, new_seed_filename)
    if write_log is True:
        log_fp.close()
    return new_seed_loc
def relabel_seeds_from_rnacentral_urs_mapping ( seed , expert_db = None , dest_dir = None , clean = False ) : if dest_dir is None : dest_dir = os . path . split ( seed ) [ 0 ] sequence_label = 0 new_line = '' write_log = False filename = os . path . split ( seed ) [ 1 ] . partition ( '.' ) [ 0 ] new_seed_filename = filename + '_relabelled' new_seed_loc = os . path . join ( dest_dir , new_seed_filename ) seed_fp = open ( seed , 'r' ) new_seed_fp = open ( new_seed_loc , 'w' ) log_fp = None sequence_count = 0 seed_seq_id = '' unique_seed_accs = { } sequence_mismatches = { } for line in seed_fp : if line [ 0 ] != '#' and len ( line ) > 1 and line [ 0 : 2 ] != '//' : line_elements = [ x for x in line . strip ( ) . split ( ' ' ) if x != '' ] sequence_count += 1 seed_seq_id = line_elements [ 0 ] . split ( '/' ) [ 0 ] seed_seq = line_elements [ 1 ] . replace ( '.' , '' ) seed_seq = seed_seq . replace ( '-' , '' ) . replace ( 'T' , 'U' ) . replace ( 't' , 'u' ) . upper ( ) rnacentral_id = map_rnacentral_urs_wirh_db_accessions ( seed_seq_id , expert_db ) if rnacentral_id is None and write_log is False : seed_filename = os . path . basename ( seed ) . partition ( '.' ) [ 0 ] log_fp = open ( os . path . join ( dest_dir , seed_filename ) + '.log' , 'w' ) write_log = True if rnacentral_id is None : log_fp . write ( "RNACENTRAL ID MISMATCH: %s\n" % seed_seq_id ) continue rnacentral_sequence = fetch_sequence_from_rnacentral ( rnacentral_id ) coordinates = fetch_seed_sequence_coordinates ( seed_seq , rnacentral_sequence ) new_label = '' if coordinates [ 1 ] != 0 : new_label = rnacentral_id + '/' + str ( coordinates [ 0 ] ) + '-' + str ( coordinates [ 1 ] ) if new_label not in unique_seed_accs : unique_seed_accs [ new_label ] = '' else : if write_log is False : seed_filename = os . path . basename ( seed ) . partition ( '.' ) [ 0 ] log_fp = open ( os . path . join ( dest_dir , seed_filename ) + '.log' , 'w' ) write_log = True log_fp . write ( "RNACENTRAL DUPLICATED ID: %s\n" % new_label ) continue else : if write_log is False : seed_filename = os . path . basename ( seed ) . partition ( '.' ) [ 0 ] log_fp = open ( os . path . join ( dest_dir , seed_filename ) + '.log' , 'w' ) write_log = True log_fp . write ( "RNACENTRAL SEQ MISMATCH: %s\n" % seed_seq_id ) sequence_mismatches [ rnacentral_id ] = rnacentral_sequence continue new_line = "\t" . join ( [ new_label , line_elements [ 1 ] , '\n' ] ) else : new_line = line new_seed_fp . write ( new_line ) seed_fp . close ( ) new_seed_fp . close ( ) sequence_count = len ( unique_seed_accs ) if not bool ( sequence_mismatches ) is False : fasta_filename = filename + '_rnac' fasta_file = write_fasta_seed_file ( sequence_mismatches , fasta_filename , dest_dir ) if fasta_file is None : sys . exit ( "FILE ERROR: Fasta file %s could not be generated\n" % fasta_filename ) cmfile = build_temporary_cm_from_seed ( seed , dest_dir ) if cmfile is None : sys . exit ( "FILE ERROR: CM file for seed %s could not be generated\n" % seed_filename ) cmaligned_sequences = align_sequences_to_cm ( cmfile , fasta_file , dest_dir ) if cmaligned_sequences is None : sys . exit ( "FILE ERROR: New cmaligned SEED for family %s could not be generated\n" % seed_filename ) if sequence_count > len ( sequence_mismatches . keys ( ) ) : filename = os . path . split ( seed ) [ 1 ] . partition ( "." ) [ 0 ] cmaligned_relabelled_seed = align_sequences_to_cm ( cmfile , new_seed_loc ) merged_seed = merge_seeds ( cmaligned_relabelled_seed , cmaligned_sequences , filename , dest_dir ) new_seed_loc = merged_seed else : os . remove ( new_seed_loc ) new_seed_loc = cmaligned_sequences os . rename ( new_seed_loc , os . path . join ( dest_dir , new_seed_filename ) ) new_seed_loc = os . path . join ( dest_dir , new_seed_filename ) if write_log is True : log_fp . close ( ) return new_seed_loc
Relabels the accessions of a SEED alignment using RNAcentral
identifiers.
Relabels the accessions of a SEED alignment using RNAcentral identifiers .
"""
    Relabels the accessions of a SEED alignment using RNAcentral
    identifiers. This is done by matching the seed sequences, with
    sequences existing in RNAcentral using md5 hashing.

    seed: A reformatted seed in Pfam format
    expert_db: An existing RNAcentral expert database
    dest_dir: The path to the destination directory. None by default
    clean: Specifies whether to clean seeds by ommitting sequences not found in RNAcentral

    return: The path to the relabelled SEED alignement
    """ # fetch path of seed alignment # dictionary to keep track of sequence mismatches # check if this is an actual sequence line # replace alignment characters # if any of the sequences isn't found prepare to log missing accessions # sets write_log to only open the file once # write log file if one exists and jump to next iteration # make sure subsequence was found # if we found a duplicated id, write log and skip this sequence # if we reached this point, it means there was an id match, but no exact sequence match # try mapping SEED sequences in smaller segments #else: # here we need to generate the new label with coords and #   print (sequence) # checks if dictionary isn't empty # write RNAcentral sequences to fasta file to be used with cmalign # generate CM file from original seed # align fasta file to covariance model # if the number of sequences in the original SEED is larger than the sequence mismatches # this means we have two smaller MSAs that need to be merged # renames files #final_seed = remove_all_gap_columns(new_seed_loc, filename, dest_dir) #if final_seed is None: #    sys.exit("FILE ERROR: SEED reformatting failed\n") # close log file if one exists
seed
expert_db
dest_dir
clean