[BioPython] Getting SNPs from a given gene

[Bonis Sanz, Julio]

Hi, 

I am writing a script with biopython. The aim is that given a gi, the scripts returns a list of the SNPs (variations) in that given gene. (I work with mRNA).

I did it some months ago and seems to work, but now, I have found that for the gi-s I try (mRNA for Homo sapiens and using only refseq as source database) there are no "variation" features in the record!.

Any idea?

I attach an example of code.

    def populateByGI(self,gi):
        from Bio import GenBank
        rParser = GenBank.FeatureParser()
        rDict = GenBank.NCBIDictionary(parser = rParser)
        sq = rDict[gi]
        self.gi = gi
        sequence = ''
        for a in sq.seq:
            sequence = sequence + a
        self.seq = sequence
        for feature in sq.features:
            if (feature.type == 'gene'):
                self.names.append(feature.qualifiers['gene'][0])
                dbrefs = feature.qualifiers['db_xref']
                for dbref in dbrefs:
                    if (dbref[:3] == 'MIM'):
                        self.mim = dbref[4:]
                    if (dbref[:7] == 'LocusID'):
                        self.locusID = dbref[8:]
            if (feature.type == 'variation'):
                newVar = GBvariation()
                newVar.gi = gi
                newVar.position = feature.location._end.position
                newVar.dbsnp = feature.qualifiers['db_xref'][0][6:]
                newVar.allele1 = feature.qualifiers['replace'][0]
                newVar.allele2 = feature.qualifiers['replace'][1]
                self.variations.append(newVar)