#*****************************************************************************
# Copyright (C) 2007 Marshall Hampton hamptonio at gmail.com
#
# Distributed under the terms of the GNU General Public License (GPL)
# http://www.gnu.org/licenses/
#*****************************************************************************
from Bio import SeqIO
from Bio import GenBank
def cg_imagemap(sterm, imagename, htmlfilename, verbose = False):
'''
Takes the top hit in GenBank (Core Nucleotide) to the search term sterm, and constructs an HMTL image map of the G+C content linking to the original record. At some point this should be refactored into smaller, cleaner code chunks. Note that the examples below use some hardcoded paths that would have to be modified for your system.
EXAMPLES:
sage: cg_imagemap("NC_001807","/Volumes/D/hum1.png","/Volumes/D/hum1.html", verbose = True)
REVIEWED REFSEQ: This record has been curated by NCBI staff. The
reference sequence was derived from AF347015.
On Dec 27, 2001 this sequence version replaced gi:13959823.
Please be aware that the NCBI Reference sequence represents the
mtDNA from an African (Yoruba) individual. The modified version of
the original Cambridge Reference Sequence (GenBank Acc: J01415,
Anderson et al (1981)) is NCBI alternative Refseq AC_000021.
COMPLETENESS: full length.
C+G averages computed
sage: cg_imagemap("Synechococcus cyanophage syn9, complete genome","/Volumes/D/tim1.png","/Volumes/D/tim1.html")
'''
gi_list = GenBank.search_for(sterm)
record_parser = GenBank.RecordParser()
ncbi_nuc_dict = GenBank.NCBIDictionary(database='nucleotide', format='genbank', parser=record_parser)
record = ncbi_nuc_dict[gi_list[0]]
if verbose: print record.comment
seq = record.sequence
cgavg = []
for i in range(len(seq)):
lower = i - 50
upper = i + 50
if lower < 0:
tempseq = seq[lower:-1] + seq[0:upper]
elif upper < len(seq):
tempseq = seq[lower:upper]
else:
tempseq = seq[lower:len(seq)] + seq[0:upper-len(seq)]
tempcg = 0
tempcg += tempseq.count('C') + tempseq.count('G')
cgavg.append([i,float(tempcg)/100.0])
if verbose: print "C+G averages computed"
circdata = []
for pair in cgavg:
loc = float(2*pi*pair[0])/len(seq)
offset = pair[1]-.5
cpoint = [(1+offset)*cos(loc),(1+offset)*sin(loc)]
circdata.append(cpoint)
imgA = point(circdata[0],pointsize=1,hue=((circdata[0][0])^2+circdata[0][1]^2)^(1/2)-.5)
for i in range(0,len(circdata),len(circdata)/3600):
imgA += point(circdata[i],pointsize=1,hue=((circdata[i][0])^2+circdata[i][1]^2)^(1/2)-.5)
imgB = line([circdata[i] for i in range(0,len(circdata),len(circdata)/3600)], thickness=.1)
img = imgA + imgB
img.save(filename=imagename, axes = False, figsize = [4,4])
htmlstr1 = '''G+C content'''
htmlstr2=''''''
tempfile = file(htmlfilename,'a')
slen = len(seq)
local_imagename = imagename.split('/')[-1]
htmlmap2 = '''
'''
tempfile.write(htmlstr1)
tempfile.write(htmlmap2)
tempfile.write(htmlstr2)
tempfile.close()