Malheureusement, vous pouvez charger un tel fichier GFF3 dans R (en utilisant le package GenomicFeatures) et accéder aux UTR, mais ils ne sont pas ensuite enregistrés si vous utilisez la fonction export de rtracklayer (cela aurait été la solution la plus simple).
Comme je ne sais rien pour faire ça de façon désinvolte, j'ai passé quelques minutes à en écrire une. L'utilisation serait addGFF3UTRs.py input.gff output.gff . Vous pouvez modifier un certain nombre d'étiquettes pour qu'elles apparaissent comme vous le souhaitez. Cela n'a pas été testé de manière approfondie, mais il semble produire les bons résultats dans un petit test que je viens d'exécuter.
Le code est hébergé sur GitHub au lien que j'ai donné ci-dessus et j'inclus également ici:
#! / usr / bin / env pythonimport argparsedef parseAttributes (kvps): d = dict () pour kvp dans kvps: k, v = kvp.split ("=") d [k] = v return ddef parseGFF3 (fname, topLevel = "gene", secondLevel = "mRNA", CDS = "CDS", exon = "exon"): f = open (fname) genes = dict () transcripts = dict () g2t = dict () # Un dictionnaire avec les ID de gène associés à leurs ID de transcription pour la ligne dans f: if line.startswith ("#"): continue cols = line.strip () .split ("\ t") attribs = parseAttributes (cols [8] .split (";")) ID = attribs ["ID"] if cols [2] == topLevel: genes [ID] = ligne g2t [ID] = [] elif cols [2] == secondLevel: parent = attribs ["Parent"] g2t [parent ] .append (ID) # Génère une liste de [ligne, [entrées exon], [entrées CDS]] si l'ID ne figure pas dans les transcriptions: transcriptions [ID] = [ligne, [], []] sinon: transcriptions [ID] [0] = line else: parent = attribs ["Parent"] si le parent n'est pas dans les transcriptions: transcriptions [parent] = [None, [], []] if cols [2] == exon: transcripts [parent] [1 ] .append ((int (cols [3]), int (cols [4]), ligne))
elif cols [2] == CDS: transcriptions [parent] [2] .append ((int (cols [3]), int (cols [4]), ligne)) f.close () retourne les gènes, transcriptions, g2tdef findUTRs (transcripts, fivePrime = True): pour k, t dans transcripts.items (): # S'il n'y a pas d'entrées CDS alors sautez si len (t [2]) == 0: continue # Récupère le brin ("." sera traité comme "+") strand = t [0] .split ("\ t") [6] # Pour 3 'UTR, permutez simplement le brin sinon cinqPrime: strand = "+" if strand == "- "else" - "exons = [(s, e) pour s, e, _ dans t [1]] CDSs = [(s, e) pour s, e, _ dans t [2]] exons.sort () CDSs.sort () UTRs = [] if strand! = "-": final = CDSs [0] [0] for s, e in exons: if e < final: UTRs.append ((s, e)) elif s < final: UTRs.append ((s, final - 1)) else: break else: final = CDSs [-1] [- 1] pour s, e dans les exons: if e < final: continue elif s > final: UTRs.append ((s, e)) else: UTRs.append ((final + 1, e) ) t.append (UTRs) def saveGFF (oname, genes, transcripts, g2t, fivePrime, threePrime): o = open (oname, "w") o.write ("## gff-version 3 \ n") pour geneID , geneLine dans genes.items (): o.write (geneLine) # Parcourez chaque transcrit pour transID dans g2t [geneID]: t = transcripts [transID] si t [0] est None: continue o.write (t [0 ]) cols = t [0] .strip (). split ("\ t") # Ecrire les exons, puis CDS, puis 5'UTR, puis 3'UTR pour l'exon dans t [1]: o.write (exon [2]) pour CDS dans t [2]: o.write (CDS [2]) pour idx, UTR dans enumerate (t [3]):
o.write ("{} \ t {} \ t {} \ t {} \ t {} \ t. \ t {} \ t. \ t" .format (cols [0], cols [1], fivePrime , UTR [0], UTR [1], cols [6])) o.write ("ID = {}. FivePrimeUTR {}; Parent = {} \ n" .format (transID, idx, transID)) pour idx , UTR dans enumerate (t [4]): o.write ("{} \ t {} \ t {} \ t {} \ t {} \ t. \ T {} \ t. \ T" .format ( cols [0], cols [1], threePrime, UTR [0], UTR [1], cols [6])) o.write ("ID = {}. threePrimeUTR {}; Parent = {} \ n". format (transID, idx, transID)) o.close () parser = argparse.ArgumentParser (description = "Analyser un fichier GFF3 sans entrées UTR et les ajouter. Notez que ce programme suppose un fichier GFF3 raisonnablement formaté") parser.add_argument ("--fiveUTRname", default = "five_prime_UTR", help = "Le libellé des 5 entrées UTR (par défaut:% (default) s)") parser.add_argument ("- threeUTRname", default = "three_prime_UTR", help = "Le libellé pour 3 'entrées UTR (par défaut:% (default) s)") parser.add_argument ("- topLevelID", default = "gene", help = "Le' type 'désignant l'entrée de niveau supérieur (par défaut:% (par défaut) s) ") parser .add_argument ("- secondLevelID", default = "mRNA", help = "Le 'type' désignant l'entrée de deuxième niveau, généralement quelque chose comme 'mRNA' ou 'transcript' (par défaut:% (par défaut) s)") parser.add_argument ("- exonID", default = "exon", help = "Le 'type' désignant les exons (par défaut:% (default) s)") parser.add_argument ("- CDSID", default = "CDS ", help =" Le 'type' désignant CDS (par défaut:% (valeur par défaut) s) ") parser.add_argument (" input ", metavar =" input.gff3 ", help =" Input file name ") parser.add_argument ( "output", metavar = "output.gff3", help = "Output file name") args = parser.parse_args () genes, transcripts, g2t = parseGFF3 (args.input, topLevel = args.topLevelID, secondLevel = args.secondLevelID , exon = args.exonID, CDS = args.CDSID) # Ajouter les UTRsfindUTRs (transcriptions) findUTRs (transcripts, fivePrime = False) # Ecrire les sortiesaveGFF (args.output, genes, transcripts, g2t, args.fiveUTRname, args.threeUTRname )