Python Programming for Biology: Bioinformatics and Beyond

GC content of DNA

The  next  example  investigates  a  DNA  sequence  by  measuring  its  GC  content:  i.e.  the

percentage of the total base pairs that are G:C (rather than A:T). All we need to do for this

---

is  to  take  the  sequence  of  one  strand  of  DNA  and  simply  count  how  many  of  the

nucleotides are G or C. Measuring the GC content of DNA is biologically relevant because

regions of a chromosome that are rich in G and C give a hint that they might be coding for

genes.

To  make  things  more  interesting,  rather  than  just  report  the  final  GC  content  for  the

whole of a sequence, we will measure the GC content for every possible 10 residue sub-

sequence  and  then  plot  the  values  along  the  length  of  the  sequence  as  a  graph.  In  other

words we will perform the calculation on a sliding window of residues.

Firstly, we define a function that takes a DNA sequence and a window size (optionally)

as input and gives a list of numerical GC content values as output. We will take the output

data and use it to draw a graph using an external Python module called Matplotlib, which

is  very  useful  for  plotting  numerical  data  (see  Chapter  9  for  more  details).  As  with  the

profile  search  above,  we  will  use  a  for  loop  to  scan  through  the  sequence,  while  taking

care to avoid falling off the end. However, this time, because we don’t need the position of

each  nucleotide  within  the  search  window  (we  only  needed  this  before  to  get  the  right

position  in  a  profile),  we  can  find  the  number  of  G  and  C  letters  by  using  the  .count()

method that is built into Python strings and lists.

def calcGcContent(seq, winSize=10):

gcValues = []

for i in range(len(seq)-winSize):

subSeq = seq[i:i+winSize]

numGc = subSeq.count('G') + subSeq.count('C')

value = numGc/float(winSize)

gcValues.append(value)

return gcValues

Each of the measurements for each sliding position are added to the output list, which

we  can  then  plot  as  a  graph  as  follows.  Note  that  this  example  assumes  that  we  have

installed  the  Matplotlib  module  (see

http://www.cambridge.org/pythonforbiology

 for

download links), otherwise you will get an error from the import command.

from matplotlib import pyplot

gcResults = calcGcContent(dnaSeq)

pyplot.plot(gcResults)

pyplot.show()