Mercurial > repos > artbio > mircounts
diff mature_mir_gff_translation.py @ 9:2a08a6eb471c draft
planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/mircounts commit 6013aaf29ff7aa2d1aab434f2355da327c7ef102
author | artbio |
---|---|
date | Wed, 25 Apr 2018 12:48:27 -0400 |
parents | 3f62272192f9 |
children | de227b7307cf |
line wrap: on
line diff
--- a/mature_mir_gff_translation.py Mon Apr 23 13:21:16 2018 -0400 +++ b/mature_mir_gff_translation.py Wed Apr 25 12:48:27 2018 -0400 @@ -13,105 +13,90 @@ return args -def get_gff_header(gff_input_file): - string_list = [] - for line in open(gff_input_file, "r"): - if line[0] == '#': - string_list.append(line) - string_list.append('# generated by mature_mir_gff_translation.py %s\n#\n' % - str(datetime.now())) - return ''.join(string_list) +def convert_and_print_gff(gff_input_file, output): + def get_gff_header(gff_input_file): + string_list = [] + for line in open(gff_input_file, "r"): + if line[0] == '#': + string_list.append(line) + string_list.append('# generated by mature_mir_gff_translation.py \ + %s\n#\n' % str(datetime.now())) + return ''.join(string_list) -def load_gff_in_dict(gff_input_file): - ''' - Reads the gff3 file and return a dictionary of dictionaries - with keys equal to standard gff3 fields (9) - Note that the key of the primary dictionary is the ID - ''' gff_dict = {} + # see https://github.com/ARTbio/tools-artbio/issues/246 + # currently fixed by perl pretreatment or the gff3 file for line in open(gff_input_file, "r"): if line[0] == "#": continue gff_fields = line[:-1].split("\t") ID = gff_fields[8].split("ID=")[1].split(";")[0] - gff_dict[ID] = {} - gff_dict[ID]["seqid"] = gff_fields[0] - gff_dict[ID]["source"] = gff_fields[1] - gff_dict[ID]["type"] = gff_fields[2] - gff_dict[ID]["start"] = gff_fields[3] - gff_dict[ID]["end"] = gff_fields[4] - gff_dict[ID]["score"] = gff_fields[5] - gff_dict[ID]["strand"] = gff_fields[6] - gff_dict[ID]["phase"] = gff_fields[7] - gff_dict[ID]["attributes"] = gff_fields[8] - if "erives_from" in gff_dict[ID]["attributes"]: - parent_primary_transcript = gff_dict[ID]["attributes"].split( - "erives_from=")[1] - parent_primary_transcript = gff_dict[parent_primary_transcript][ - "attributes"].split("Name=")[1] - gff_dict[ID]["attributes"] = "%s;Parent_mir_Name=%s" % ( - gff_dict[ID]["attributes"], parent_primary_transcript) - return gff_dict - - -def genome_to_mir_gff(gff_dict, output, header): - ''' - Converts seqid field from chromosome to item Name - Then converts coordinates relative to "miRNA_primary_transcript" - Note that GFF files are 1-based coordinates - ''' - for key in gff_dict: - name = gff_dict[key]["attributes"].split("Name=")[1].split(";")[0] - gff_dict[key]["seqid"] = name - if "erives_from=" in gff_dict[key]["attributes"]: - parent_ID = gff_dict[key]["attributes"].split( - "erives_from=")[1].split(";")[0] - gff_dict[key]["start"] = str(int(gff_dict[key]["start"])-int( - gff_dict[parent_ID]["start"])+1) - gff_dict[key]["end"] = str(int(gff_dict[key]["end"])-int( - gff_dict[parent_ID]["start"])+1) - hairpins = {} - matures = {} - # treats miRNA_primary_transcript coordinates - # in a second loop to avoid errors in conversion - for key in gff_dict: - if gff_dict[key]["type"] == "miRNA_primary_transcript": - gff_dict[key]["end"] = str(int(gff_dict[key]["end"])-int( - gff_dict[key]["start"]) + 1) - gff_dict[key]["start"] = '1' - # now, do a dict[ID]=Name but only for miRNA_primary_transcript - hairpins[key] = gff_dict[key]["attributes"].split( - "Name=")[1].split( - ";")[0] + if gff_fields[2] == "miRNA_primary_transcript": + gff_dict[ID] = {} + gff_dict[ID]["premir_name"] = gff_fields[8].split( + "Name=")[1].split(";")[0] + gff_dict[ID]["primary"] = line[:-1] + gff_dict[ID]["miRNAs"] = [] + elif gff_fields[2] == "miRNA": + parent_ID = gff_fields[8].split("erives_from=")[1] + gff_dict[parent_ID]["miRNAs"].append(line[:-1]) + # Now reorganise features and recalculate coordinates of premirs and mirs + gff_list = [] + for ID in sorted(gff_dict, key=lambda x: (gff_dict[x]['premir_name'])): + # delete premir and their mir with ID containing "_" + if "_" in ID: + del gff_dict[ID] else: - matures[key] = gff_dict[key]["attributes"].split( - "Name=")[1].split( - ";")[0] + primary_fields = gff_dict[ID]["primary"].split('\t') + seqid = primary_fields[8].split("Name=")[1].split(";")[0] + source = primary_fields[1] + type = primary_fields[2] + start = primary_fields[3] + newstart = "1" + end = primary_fields[4] + newend = str(int(end)-int(start)+1) + score = primary_fields[5] + strand = primary_fields[6] + phase = primary_fields[7] + attributes = primary_fields[8] + gff_list.append('%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s' % (seqid, + source, type, newstart, newend, score, strand, + phase, attributes)) + # ensure their is only 2 child miRNAs at best + if len(gff_dict[ID]["miRNAs"]) > 2: + gff_dict[ID]["miRNAs"] = gff_dict[ID]["miRNAs"][:2] + # sort child miRNAs 5p first 3p second + if gff_dict[ID]["miRNAs"][0].find('5p') == -1: + gff_dict[ID]["miRNAs"] = gff_dict[ID]["miRNAs"][::-1] + for mir in gff_dict[ID]["miRNAs"]: + mir_fields = mir.split('\t') + mir_seqid = mir_fields[8].split("Name=")[1].split(";")[0] + mir_source = mir_fields[1] + mir_type = mir_fields[2] + mir_start = mir_fields[3] + mir_end = mir_fields[4] + new_mir_start = str(int(mir_start)-int(start)+1) + new_mir_end = str(int(mir_end)-int(start)+1) + mir_score = mir_fields[5] + mir_strand = mir_fields[6] + mir_phase = mir_fields[7] + mir_attributes = mir_fields[8] + mir_sfx = ";Parent_mir_Name=%s" % gff_dict[ID]["premir_name"] + gff_list.append( + '%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s%s' % ( + mir_seqid, mir_source, mir_type, + new_mir_start, new_mir_end, mir_score, + mir_strand, mir_phase, mir_attributes, + mir_sfx)) with open(output, "w") as output: - output.write(header) - for ID in sorted(hairpins, key=hairpins.get): - output.write("\t".join([gff_dict[ID]["seqid"], - gff_dict[ID]["source"], gff_dict[ID]["type"], - gff_dict[ID]["start"], gff_dict[ID]["end"], - gff_dict[ID]["score"], gff_dict[ID]["strand"], - gff_dict[ID]["phase"], gff_dict[ID]["attributes"]])) - output.write("\n") - for id in sorted(matures, key=matures.get, reverse=True): - if ID in gff_dict[id]["attributes"]: - output.write("\t".join([gff_dict[id]["seqid"], - gff_dict[id]["source"], gff_dict[id]["type"], - gff_dict[id]["start"], gff_dict[id]["end"], - gff_dict[id]["score"], - gff_dict[id]["strand"], - gff_dict[id]["phase"], - gff_dict[id]["attributes"]])) - output.write("\n") + output.write('%s' % get_gff_header(gff_input_file)) + output.write('\n'.join(gff_list)) + output.write('\n') def main(infile, outfile): - gff_dict = load_gff_in_dict(infile) - genome_to_mir_gff(gff_dict, outfile, get_gff_header(infile)) + convert_and_print_gff(infile, outfile) if __name__ == "__main__":