Skip to content
Snippets Groups Projects
Commit 3dcfdb96 authored by Laros's avatar Laros
Browse files

Created branch lrg. 


git-svn-id: https://humgenprojects.lumc.nl/svn/mutalyzer/branches/lrg@63 eb6bd6ab-9ccd-42b9-aceb-e2899b4a52f1
parents 3c895565 8eef6214
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Db.txt
+ 15
− 0
View file @ 3dcfdb96
......@@ -26,6 +26,7 @@ In database mutalyzer:
---
CREATE TABLE GBInfo (
AccNo CHAR(20) PRIMARY KEY,
GI CHAR(13) UNIQUE NOT NULL,
hash CHAR(32) UNIQUE NOT NULL,
ChrAccVer CHAR(20),
ChrStart INT(12),
......@@ -49,6 +50,20 @@ CREATE TABLE BatchJob (
FromHost Char(255) NOT NULL
);
CREATE TABLE Link (
mrnaAcc CHAR(20) PRIMARY KEY,
protAcc CHAR(20) UNIQUE NOT NULL
);
INSERT INTO Link
SELECT DISTINCT acc, protAcc
FROM hg18.map
WHERE protAcc != ''
UNION
SELECT DISTINCT acc, protAcc
FROM hg19.map
WHERE protAcc != '';
CREATE TABLE Var (
Name CHAR(20) PRIMARY KEY,
Value Char(255) NOT NULL
......
src/Modules/Config.py
+ 9
− 7
View file @ 3dcfdb96
......@@ -20,7 +20,7 @@ class Config() :
Mutator ; Container for the Mutator configuration variables.
Scheduler ; Container for the Scheduler configuration variables.
File ; Container for the File configuration variables.
GenRecord ; Container for the File configuration variables.
GBparser ; Container for the File configuration variables.
Special Methods:
__init__ ; Read the configuration file and initialise the
......@@ -118,9 +118,9 @@ class Config() :
pass
#File
class GenRecord() :
class GBparser() :
"""
Container class for the GenRecord configuration variables.
Container class for the GBparser configuration variables.
Public variables:
upstream ; Number of upstream nucleotides when searching for a
......@@ -160,7 +160,6 @@ class Config() :
# Set the variables needed by the Db module.
self.Db.internalDb = config["internalDb"]
self.Db.dbNames = config["dbNames"]
self.Db.LocalMySQLuser = config["LocalMySQLuser"]
self.Db.LocalMySQLhost = config["LocalMySQLhost"]
self.Db.RemoteMySQLuser = config["RemoteMySQLuser"]
......@@ -191,9 +190,12 @@ class Config() :
self.File.header = config["header"]
self.File.tempDir = config["tempDir"]
# Set the variables needed by the GenRecord module.
self.File.upstream = int(config["upstream"])
self.File.downstream = int(config["downstream"])
# Set the variables needed by the GBparser module.
self.GBparser.email = config["email"]
## Set the variables needed by the File module.
#self.File.upstream = int(config["upstream"])
#self.File.downstream = int(config["downstream"])
#__init__
#Config
......
src/Modules/Db.py
+ 45
− 0
View file @ 3dcfdb96
......@@ -1040,6 +1040,51 @@ class Cache(Db) :
return ret[0][0]
return None
#getGI
def getProtAcc(self, mrnaAcc) :
"""
"""
statement = """
SELECT protAcc
FROM Link
WHERE mrnaAcc = %s;
""", mrnaAcc
ret = self.query(statement)
if ret :
return ret[0][0]
return None
#getProtAcc
def getmrnaAcc(self, protAcc) :
"""
"""
statement = """
SELECT mrnaAcc
FROM Link
WHERE protAcc = %s;
""", protAcc
ret = self.query(statement)
if ret :
return ret[0][0]
return None
#getProtAcc
def insertLink(self, mrnaAcc, protAcc) :
"""
"""
statement = """
INSERT INTO Link
VALUES (%s, %s);
""", (mrnaAcc, protAcc)
self.query(statement)
#insertLink
#Cache
class Batch(Db) :
......
src/Modules/GBparser.py
+ 249
− 176
View file @ 3dcfdb96
from Bio import SeqIO # read()
import bz2 # BZ2Compressor(), BZ2File()
from Bio import SeqIO, Entrez # read()
import bz2 # BZ2Compressor(), BZ2File()
from GenRecord import PList, Locus, Gene, Record, GenRecord
import Db
"""
Module contains one public function createGBRecord which returns a
mutalyzer GenRecord.Record populated with data from a GenBank file.
"""
def __location2pos(location) :
class GBparser() :
"""
Convert a location object to a tuple of integers.
Arguments:
location ; A location object (see the BioPython documentation).
Returns:
List ; A tuple of integers.
"""
ret = []
if not str(location.start).isdigit() or \
not str(location.end).isdigit() :
return None
#if
ret.append(location.start.position + 1)
ret.append(location.end.position)
return ret
#__location2pos
def __locationList2posList(locationList) :
"""
Convert a list of locations to a list of integers.
Arguments:
locationList ; A list of locations
(see the BioPython documentation).
Returns:
List ; A list (of even length) of integers.
"""
ret = []
if not str(locationList.location.start).isdigit() or \
not str(locationList.location.end).isdigit() :
return None
#if
for i in locationList.sub_features :
if i.ref : # This is a workaround for a bug in BioPython.
ret = None
break
def __init__(self) :
"""
"""
import Config
config = Config.Config()
Entrez.email = config.Retriever.email
self.__database = Db.Cache(config.Db)
#__init__
def __location2pos(self, location) :
"""
Convert a location object to a tuple of integers.
Arguments:
location ; A location object (see the BioPython documentation).
Returns:
List ; A tuple of integers.
"""
ret = []
if not str(location.start).isdigit() or \
not str(location.end).isdigit() :
return None
#if
temp = __location2pos(i.location)
ret.append(temp[0])
ret.append(temp[1])
#for
return ret
#__locationList2posList
def createGBRecord(filename):
"""
Create a GenRecord.Record from a GenBank file
Arguments:
filename ; The full path to the compresed GenBank file
Returns:
record ; A GenRecord.Record instance
"""
# first create an intermediate genbank record with BioPython
file_handle = bz2.BZ2File(filename, "r")
biorecord = SeqIO.read(file_handle, "genbank")
file_handle.close()
record = Record()
record.seq = biorecord.seq
mRNAProducts = []
CDSProducts = []
for i in biorecord.features :
if i.qualifiers :
if i.qualifiers.has_key("gene") :
if i.type == "mRNA" :
if i.qualifiers.has_key("product") :
mRNAProducts.append(i.qualifiers["product"][0])
if i.type == "CDS" :
if i.qualifiers.has_key("product") :
CDSProducts.append(i.qualifiers["product"][0])
#if
for i in biorecord.features :
if i.qualifiers :
if i.type == "source" :
if i.qualifiers.has_key("mol_type") :
if i.qualifiers["mol_type"][0] == "mRNA" :
record.molType = 'n'
else :
record.molType = 'g'
#if
if i.qualifiers.has_key("organelle") :
record.organelle = i.qualifiers["organelle"][0]
if record.organelle == "mitochondrion" :
record.molType = 'm'
#if
fakeGene = Locus("001")
record.source.transcriptList.append(fakeGene)
fakeGene.CDS = PList()
fakeGene.CDS.location = __location2pos(i.location)
ret.append(location.start.position + 1)
ret.append(location.end.position)
return ret
#__location2pos
def __locationList2posList(self, locationList) :
"""
Convert a list of locations to a list of integers.
Arguments:
locationList ; A list of locations
(see the BioPython documentation).
Returns:
List ; A list (of even length) of integers.
"""
ret = []
if not str(locationList.location.start).isdigit() or \
not str(locationList.location.end).isdigit() :
return None
#if
for i in locationList.sub_features :
if i.ref : # This is a workaround for a bug in BioPython.
ret = None
break
#if
if i.qualifiers.has_key("gene") :
gene = i.qualifiers["gene"][0]
GeneInstance = record.findGene(gene)
if not GeneInstance :
GeneInstance = Gene(gene)
record.geneList.append(GeneInstance)
#if
if i.type == "gene" :
if i.strand :
GeneInstance.orientation = i.strand
GeneInstance.location = __location2pos(i.location)
#if
# Look if there is a locus tag present, if not, give it the
# default tag `001'.
locusName = "001"
locusTag = None
if i.qualifiers.has_key("locus_tag") :
locusTag = i.qualifiers["locus_tag"][0]
locusName = locusTag[-3:]
#if
LocusInstance = GeneInstance.findLocus(locusName)
if not LocusInstance :
LocusInstance = Locus(locusName)
GeneInstance.transcriptList.append(LocusInstance)
#if
if i.type == "mRNA" :
PListInstance = PList()
LocusInstance.mRNA = PListInstance
posList = __locationList2posList(i)
if posList != None :
PListInstance.location = \
__location2pos(i.location)
PListInstance.positionList = posList
temp = self.__location2pos(i.location)
ret.append(temp[0])
ret.append(temp[1])
#for
return ret
#__locationList2posList
def __transcriptToProtein(self, transcriptAcc) :
"""
"""
proteinAcc = self.__database.getProtAcc(transcriptAcc)
if not proteinAcc :
handle = Entrez.esearch(db = "nucleotide", term = transcriptAcc)
result = Entrez.read(handle)
handle.close()
transcriptGI = result["IdList"][0]
handle = Entrez.elink(dbfrom = "nucleotide", db = "protein",
id = transcriptGI)
result = Entrez.read(handle)
handle.close()
proteinGI = result[0]["LinkSetDb"][0]["Link"][0]["Id"]
handle = Entrez.efetch(db = "protein", id = proteinGI,
rettype = "acc")
proteinAcc = handle.read().split('.')[0]
handle.close()
self.__database.insertLink(transcriptAcc, proteinAcc)
#if
return proteinAcc
#__transcriptToProtein
def createGBRecord(self, filename):
"""
Create a GenRecord.Record from a GenBank file
Arguments:
filename ; The full path to the compresed GenBank file
Returns:
record ; A GenRecord.Record instance
"""
# first create an intermediate genbank record with BioPython
file_handle = bz2.BZ2File(filename, "r")
biorecord = SeqIO.read(file_handle, "genbank")
file_handle.close()
record = Record()
record.seq = biorecord.seq
#mRNAProducts = []
#CDSProducts = []
#for i in biorecord.features :
# if i.qualifiers :
# if i.qualifiers.has_key("gene") :
# if i.type == "mRNA" :
# if i.qualifiers.has_key("product") :
# mRNAProducts.append(i.qualifiers["product"][0])
# if i.type == "CDS" :
# if i.qualifiers.has_key("product") :
# CDSProducts.append(i.qualifiers["product"][0])
# #if
#print mRNAProducts
#print CDSProducts
for i in biorecord.features :
if i.qualifiers :
if i.type == "source" :
if i.qualifiers.has_key("mol_type") :
if i.qualifiers["mol_type"][0] == "mRNA" :
record.molType = 'n'
else :
record.molType = 'g'
#if
if i.qualifiers.has_key("transcript_id") :
LocusInstance.transcriptID = \
i.qualifiers["transcript_id"][0]
LocusInstance.locusTag = locusTag
#if
if i.type == "CDS" :
PListInstance = PList()
LocusInstance.CDS = PListInstance
PListInstance.location = __location2pos(i.location)
PListInstance.positionList = \
__locationList2posList(i)
if i.qualifiers.has_key("transl_table") :
LocusInstance.txTable = \
int(i.qualifiers["transl_table"][0])
if i.qualifiers.has_key("protein_id") :
LocusInstance.proteinID = \
i.qualifiers["protein_id"][0]
LocusInstance.locusTag = locusTag
#if
if i.type == "exon" :
if not LocusInstance.exon :
LocusInstance.exon = PList()
LocusInstance.exon.positionList.extend(
__location2pos(i.location))
if i.qualifiers.has_key("organelle") :
record.organelle = i.qualifiers["organelle"][0]
if record.organelle == "mitochondrion" :
record.molType = 'm'
#if
fakeGene = Locus("001")
record.source.transcriptList.append(fakeGene)
fakeGene.CDS = PList()
fakeGene.CDS.location = self.__location2pos(i.location)
#if
#if
#if
#for
return record
#parseRecord
if i.qualifiers.has_key("gene") :
gene = i.qualifiers["gene"][0]
GeneInstance = record.findGene(gene)
if not GeneInstance :
GeneInstance = Gene(gene)
record.geneList.append(GeneInstance)
#if
if i.type == "gene" :
if i.strand :
GeneInstance.orientation = i.strand
GeneInstance.location = self.__location2pos(i.location)
if not GeneInstance.location :
GeneInstance.transcribe = False
#if
# RESOLV
LocusInstance = None
locusTag = None
if i.qualifiers.has_key("locus_tag") :
locusTag = i.qualifiers["locus_tag"][0]
locusName = locusTag[-3:]
LocusInstance = GeneInstance.findLocus(locusTag[-3:])
#if
else :
if i.qualifiers.has_key("transcript_id") :
LocusInstance = GeneInstance.findLink(
self.__transcriptToProtein(
i.qualifiers["transcript_id"][0].split('.')[0]))
if i.qualifiers.has_key("protein_id") :
LocusInstance = GeneInstance.findLink(
i.qualifiers["protein_id"][0].split('.')[0])
#else
if not LocusInstance and (i.type == "mRNA" or i.type == "CDS") :
if record.molType != 'n' :
LocusInstance = Locus(GeneInstance.newLocusTag())
GeneInstance.transcriptList.append(LocusInstance)
else :
LocusInstance = GeneInstance.transcriptList[0]
if not LocusInstance and i.type == "exon" :
if GeneInstance.transcriptList :
LocusInstance = GeneInstance.transcriptList[0]
else :
LocusInstance = Locus(GeneInstance.newLocusTag())
GeneInstance.transcriptList.append(LocusInstance)
# /RESOLV
if i.type == "mRNA" :
PListInstance = PList()
LocusInstance.mRNA = PListInstance
posList = self.__locationList2posList(i)
if posList != None :
PListInstance.location = \
self.__location2pos(i.location)
PListInstance.positionList = posList
#if
if i.qualifiers.has_key("transcript_id") :
LocusInstance.transcriptID = \
i.qualifiers["transcript_id"][0]
LocusInstance.link = self.__transcriptToProtein(
LocusInstance.transcriptID.split('.')[0])
LocusInstance.locusTag = locusTag
#if
if i.type == "CDS" :
PListInstance = PList()
LocusInstance.CDS = PListInstance
PListInstance.location = self.__location2pos(i.location)
PListInstance.positionList = \
self.__locationList2posList(i)
if i.qualifiers.has_key("transl_table") :
LocusInstance.txTable = \
int(i.qualifiers["transl_table"][0])
if i.qualifiers.has_key("protein_id") :
LocusInstance.proteinID = \
i.qualifiers["protein_id"][0]
LocusInstance.link = \
LocusInstance.proteinID.split('.')[0]
LocusInstance.locusTag = locusTag
#if
if i.type == "exon" :
if not LocusInstance.exon :
LocusInstance.exon = PList()
LocusInstance.exon.positionList.extend(
self.__location2pos(i.location))
#if
#if
#if
#for
return record
#parseRecord
#GBparser
src/Modules/GenRecord.py
+ 25
− 2
View file @ 3dcfdb96
......@@ -2,6 +2,7 @@
import Crossmap
import Bio
import Db
"""
Module to convert a GenBank record to a nested dictionary consisting of
......@@ -91,8 +92,11 @@ class Locus(object) :
self.proteinDescription = "?"
self.proteinRange = []
self.locusTag = None
self.link = None
self.transLongName = ""
self.protLongName = ""
self.transcribe = True
self.translate = True
#__init__
def addToDescription(self, rawVariant) :
......@@ -134,8 +138,18 @@ class Gene(object) :
self.transcriptList = []
self.location = []
self.longName = ""
self.__locusTag = "000"
#__init__
def newLocusTag(self) :
"""
"""
self.__locusTag = "%03i" % (int(self.__locusTag) + 1)
return self.__locusTag
#newLocusTag
def findLocus(self, name) :
"""
"""
......@@ -155,6 +169,16 @@ class Gene(object) :
ret.append(i.name)
return ret
#listLoci
def findLink(self, protAcc) :
"""
"""
for i in self.transcriptList :
if i.link == protAcc :
return i
return None
#findCDS
#Gene
class Record(object) :
......@@ -243,11 +267,10 @@ class GenRecord() :
checkRecord() ; Check and repair self.record
"""
def __init__(self, config, output) :
def __init__(self, output) :
"""
"""
self.__config = config
self.__output = output
self.record = None
#__init__
......
src/Modules/Retriever.py
+ 2
− 1
View file @ 3dcfdb96
......@@ -561,7 +561,8 @@ class GenBankRetriever(Retriever):
if filename is None: return None
# Now we have the file, so we can parse it.
record = GBparser.createGBRecord(filename)
GenBankParser = GBparser.GBparser()
record = GenBankParser.createGBRecord(filename)
return record
#loadrecord
......
src/Mutalyzer.py
+ 84
− 46
View file @ 3dcfdb96
......@@ -356,7 +356,7 @@ def findFrameShift(str1, str2) :
lcp = __lcp(str1, str2)
return "p.(%s%i%sfs*%i)" % (seq3(str1[lcp]), lcp + 1, seq3(str2[lcp]),
len(str2) - lcp)
len(str2) - lcp + 1)
#findFrameShift
def __toProtDescr(CDSStop, orig, trans) :
......@@ -491,9 +491,27 @@ def checkDeletionDuplication(start_g, end_g, mutationType, MUU,
"""
roll = __roll(MUU.orig, start_g, end_g)
if roll[1] :
newStart = start_g + roll[1]
newStop = end_g + roll[1]
shift = roll[1]
if GenRecordInstance.record.molType == 'n' :
mRNA = GenRecordInstance.record.geneList[0].transcriptList[0
].mRNA.positionList
print mRNA
print end_g - roll[0], end_g + roll[1]
for i in mRNA :
if end_g <= i and end_g + roll[1] > i :
print "ALARM"
shift = i - end_g
print shift
break
#if
#for
#if
if shift :
newStart = start_g + shift
newStop = end_g + shift
O.addMessage(__file__, 2, "WROLL",
"Sequence %s at position %i_%i was given, however, " \
"the HGVS notation prescribes that it should be %s at " \
......@@ -506,7 +524,7 @@ def checkDeletionDuplication(start_g, end_g, mutationType, MUU,
else :
MUU.dupM(start_g, end_g)
GenRecordInstance.name(start_g, end_g, mutationType, "", "",
roll)
(roll[0], shift))
#checkDeletionDuplication
def checkInversion(start_g, end_g, MUU, GenRecordInstance, O) :
......@@ -553,23 +571,42 @@ def checkInsertion(start_g, end_g, Arg1, MUU, GenRecordInstance, O) :
newStart = MUU.shiftpos(start_g)
newStop = MUU.shiftpos(start_g) + insertionLength
roll = __roll(MUU.mutated, newStart + 1, newStop)
if roll[0] + roll[1] >= insertionLength :
#roll = __roll(MUU.orig, start_g, end_g)
shift = roll[1]
if GenRecordInstance.record.molType == 'n' :
mRNA = GenRecordInstance.record.geneList[0].transcriptList[0
].mRNA.positionList
print mRNA
print newStop - roll[0], newStop + roll[1]
for i in mRNA :
if newStop <= i and newStop + roll[1] > i :
print "ALARM"
shift = i - newStop
print shift
break
#if
#for
#if
if roll[0] + shift >= insertionLength :
O.addMessage(__file__, 2, "WINSDUP",
"Insertion of %s at position %i_%i was given, " \
"however, the HGVS notation prescribes that it should be a " \
"duplication of %s at position %i_%i." % (
MUU.mutated[newStart:newStop], start_g, start_g + 1,
MUU.mutated[newStart:newStop], start_g + roll[1],
start_g + roll[1] + insertionLength - 1))
end_g += roll[1] - 1
MUU.mutated[newStart:newStop], start_g + shift,
start_g + shift + insertionLength - 1))
end_g += shift - 1
start_g = end_g - insertionLength + 1
GenRecordInstance.name(start_g, end_g, "dup", "", "",
(roll[0] + roll[1] - insertionLength, 0))
(roll[0] + shift - insertionLength, 0))
#if
else :
GenRecordInstance.name(start_g, start_g + 1, "ins",
MUU.mutated[newStart + roll[1]:newStop + roll[1]] , "",
roll)
GenRecordInstance.name(start_g, start_g + 1, "ins",
MUU.mutated[newStart + shift:newStop + shift] , "",
(roll[0], shift))
#checkInsertion
def __rv(MUU, RawVar, GenRecordInstance, parts, O, transcript) :
......@@ -828,7 +865,7 @@ def process(cmd, C, O) :
del retriever
del D
GenRecordInstance = GenRecord.GenRecord(C.GenRecord, O)
GenRecordInstance = GenRecord.GenRecord(O)
GenRecordInstance.record = record
GenRecordInstance.checkRecord()
#NOTE: GenRecordInstance is carrying the sequence in .record.seq
......@@ -839,38 +876,39 @@ def process(cmd, C, O) :
# PROTEIN
for i in GenRecordInstance.record.geneList :
if i.location :
for j in i.transcriptList :
if not ';' in j.description and j.CDS :
#print i.name, j.name, j.CDS.positionList, j.CDS.location
cds = Seq(str(__splice(MUU.orig, j.CDS.positionList)),
IUPAC.unambiguous_dna)
cdsm = Seq(str(__nsplice(MUU.mutated,
MUU.newSplice(j.mRNA.positionList),
MUU.newSplice(j.CDS.location),
j.CM.orientation)),
IUPAC.unambiguous_dna)
if j.CM.orientation == -1 :
cds = Bio.Seq.reverse_complement(cds)
cdsm = Bio.Seq.reverse_complement(cdsm)
#if
#if '*' in cds.translate()[:-1] :
# O.addMessage(__file__, 3, "ESTOP",
# "In frame stop codon found.")
# return
##if
orig = cds.translate(table = j.txTable, cds = True,
to_stop = True)
trans = cdsm.translate(table = j.txTable,
to_stop = True)
cdsLen = __cdsLen(MUU.newSplice(j.CDS.positionList))
#print cdsLen, len(__splice(MUU.orig, MUU.newSplice(j.CDS.positionList)))
j.proteinDescription = __toProtDescr(cdsLen, orig,
trans)
#print i.location, i.transcriptList
#if i.location :
for j in i.transcriptList :
if not ';' in j.description and j.CDS :
#print i.name, j.name, j.CDS.positionList, j.CDS.location
cds = Seq(str(__splice(MUU.orig, j.CDS.positionList)),
IUPAC.unambiguous_dna)
cdsm = Seq(str(__nsplice(MUU.mutated,
MUU.newSplice(j.mRNA.positionList),
MUU.newSplice(j.CDS.location),
j.CM.orientation)),
IUPAC.unambiguous_dna)
if j.CM.orientation == -1 :
cds = Bio.Seq.reverse_complement(cds)
cdsm = Bio.Seq.reverse_complement(cdsm)
#if
#if '*' in cds.translate()[:-1] :
# O.addMessage(__file__, 3, "ESTOP",
# "In frame stop codon found.")
# return
##if
orig = cds.translate(table = j.txTable, cds = True,
to_stop = True)
trans = cdsm.translate(table = j.txTable,
to_stop = True)
cdsLen = __cdsLen(MUU.newSplice(j.CDS.positionList))
#print cdsLen, len(__splice(MUU.orig, MUU.newSplice(j.CDS.positionList)))
j.proteinDescription = __toProtDescr(cdsLen, orig,
trans)
# /PROTEIN
......
templates/gbupload.html
+ 5
− 2
View file @ 3dcfdb96
......@@ -89,8 +89,11 @@
<input type='text' name="start"><br>
Stop Position<br>
<input type='text' name="stop"><br>
Orientation (1 or -1)<br>
<input type='text' name="orientation"><br>
Orientation<br>
<select>
<option value="1">Forward</option>
<option value="-1">Reverse</option>
</select>
</span>
<input type='hidden' name="method" value="gbupload">
<input type='hidden' name="client" value="human">
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment