Discussion :

[mycaweb]

Bonjour, j'ai hérité d'un script super util pour récupérer des séquences d'ADN depuis GenBank (une base de données) à partir des numéros d'accessions. Le script fonctionnait encore il y a peu, mais j'ai l'impression qu'une mise à jour de BioPython a modifié l'utilisation de la fonction SeqIO.read, mais je ne vois pas trop la solution.

Voici le code, il est relativement simple et se lance comme ça dans le terminal :

./code.py JN403374-JN406266 mon_mail@gmail.com

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#!/usr/bin/python
#-*- coding: utf-8 -*-
 
import re
import sys
from Bio import Entrez
from Bio import SeqIO
 
 
###variables########
Entrez.email=sys.argv[2]##ton adresse mail
nums=sys.argv[1]##JN403374-JN406266 pour récupérer (par exemple) les séquences nucléotidiques JN403374, JN403375, .. ,  JN406266
################
 
var=""
split=nums.split("-")
head=re.match("^([^\d]+)(\d+)$",split[0]).group(1)
start=int(re.match("^([^\d]+)(\d+)$",split[0]).group(2))
stop=int(re.match("^([^\d]+)(\d+)$",split[1]).group(2))
 
for i in range (start,stop+1):
	gi="{0}{1}".format(head,i)
	handle=Entrez.efetch(db='nucleotide',id=gi,rettype="gb")# Accession id works, returns genbank format, looks in the 'nucleotide' database
	record = SeqIO.read(handle, "genbank")
	handle.close()	
	var=var+">{0}\n{1}\n".format(gi,record.seq)
 
 
######menu au choix######
##output dans terminal##
#print var
##output dans fichier###
file=open("{0}.fasta".format(nums),'w')
file.write(var)
file.close()
######################

Désolé pour la question un peu c** c**
Myca.

[utopman]

Bonjour mycaweb,
j'ai testé ton code, il semble bon.
l'erreur dit que le "handle" n'est pas bon, ce qui veut dire que c'est la ligne

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
handle=Entrez.efetch(db='nucleotide',id=gi,rettype="gb")

qui pose problème. les paramètres sont ils bon et d'actualité ? (nucleotide et gb)

ces elements : JN403374-JN406266, sont ils bien référencés dans la base de données?

Cordialement.

[mycaweb]

Bonjour,
depuis des changements récents de la part de GenBank, il faut modifier

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

handle=Entrez.efetch(db='nucleotide',id=gi,rettype="gb")

par

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

handle=Entrez.efetch(db='nucleotide',id=gi,rettype="gb",retmode="text")

Merci à l'auteur du code pour sa réactivité et à utopman de sa réponse,
Myca