Harvard:Biophysics 101/Notebook:ZS/2007-2-6

From OpenWetWare

Jump to: navigation, search

Contents

Assignment 1, due 2/6/07

Link to Assignment 1

Original

Output of original code (for Homo sapiens topoisomerase (DNA) II binding protein 1, mRNA link here):

 GenBank id: BC126209.1
 total sequence length: 5148
 method 1
 A 1669
 AA 427
 AAA 150
 AAAA 48
 AAAAA 16
 AAAAAA 8
 AAAAAAA 4
 AAAAAAAA 1
 AAAAAAAAA 0
 
 method 2
 A 1669
  AA 589
 AAA 218
 AAAA 76
 AAAAA 29
 AAAAAA 13
 AAAAAAA 5
 AAAAAAAA 1
 AAAAAAAAA 0

Different GI

For modified version to display subsection of cyanobacteria PCC7942 genome (GI: 23957787) : Image:2607 kaiC 1.txt


Output:

 GenBank id: AY120853.2
 total sequence length: 74355
 method 1
 A 16631
 AA 3436
 AAA 845
 AAAA 210
 AAAAA 57
 AAAAAA 16
 AAAAAAA 5
 AAAAAAAA 0
 AAAAAAAAA 0
 
 method 2
 A 16631
 AA 4327
 AAA 1117
 AAAA 288
 AAAAA 78
 AAAAAA 21
 AAAAAAA 5
 AAAAAAAA 0
 AAAAAAAAA 0

Poly-T instead of Poly-A

For code: Image:Code1 forT.txt

Output:

 GenBank id: BC126209.1
 total sequence length: 5148
 method 1
 T 1387
 TT 309
 TTT 102
 TTTT 31
 TTTTT 10
 TTTTTT 4
 TTTTTTT 1
 TTTTTTTT 0
 TTTTTTTTT 0
 
 method 2
 T 1387
 TT 418
 TTT 144
 TTTT 46
 TTTTT 15
 TTTTTT 5
 TTTTTTT 1
 TTTTTTTT 0
 TTTTTTTTT 0

Translating protein sequence

The tricky part with translating the sequence is that the coding sequence is buried in the data from NCBI; there needs to be a second parsing layer to find the start codon, extract the CDS, and find the end codon at which to terminate the signal. For this, an extra parsing layer needs to be made; the code is generalizable for any sequence with 1 CDS in it.

For code:Image:Code fortranslating.txt

Output (sorry for formatting, but it matches NCBI):

MSRNDKEPFFVKFLKSSDNSKCFFKALESIKEFQSEEYLQIITEEEALKIKENDRSLYICDPFSGVVFDHLKKLGCRIVGPQVVIFCMHHQRCVPRAEHPVYNMVMSDVTISCTSLEKEKREEVHKYVQMMGGRVYRDLNVSVTHLIAGEVGSKKYLVAANLKKPILLPSWIKTLWEKSQEKKITRYTDINMEDFKCPIFLGCIICVTGLCGLDRKEVQQLTVKHGGQYMGQLKMNECTHLIVQEPKGQKYECAKRWNVHCVTTQWFFDSIEKGFCQDESIYKTEPRPEAKTMPNSSTPTSQINTIDSRTLSDVSNISNINASCVSESICNSLNSKLEPTLENLENLDVSAFQAPEDLLDGCRIYLCGFSGRKLDKLRRLINSGGGVRFNQLNEDVTHVIVGDYDDELKQFWNKSAHRPHVVGAKWLLECFSKGYMLSEEPYIHANYQPVEIPVSHKPESKAALLKKKNSSFSKKDFAPSEKHEQADEDLLSQYENGSSTVVEAKTSEARPFNDSTHAEPLNDSTHISLQEENQSSVSHCVPDVSTITEEGLFSQKSFLVLGFSNENESNIANIIKENAGKIMSLLSRTVADYAVVPLLGCEVEATVGEVVTNTWLVTCIDYQTLFDPKSNPLFTPVPVMTGMTPLEDCVISFSQCAGAEKESLTFLANLLGASVQEYFVRKSNAKKGMFASTHLILKERGGSKYEAAKKWNLPAVTIAWLLETARTGKRADESHFLIENSTKEERSLETEITNGINLNSDTAEHPGTRLQTHRKTVVTPLDMNRFQSKAFRAVVSQHARQVAASPAVGQPLQKEPSLHLDTPSKFLSKDKLFKPSFDVKDALAALETPGRPSQQKRKPSTPLSEVIVKNLQLALANSSRNAVALSASPQLKEAQSEKEEAPKPLHKVVVCVSKKLSKKQSELNGIAASLGADYRWSFDETVTHFIYQGRPNDTNREYKSVKERGVHIVSEHWLLDCAQECKHLPESLYPHTYNPKMSLDISAVQDGRLCNSRLLSAVSSTKDDEPDPLILEENDVDNMATNNKESAPSNGSGKNDSKGVLTQTLEMRENFQKQLQEIMSATSIVKPQGQRTSLSRSGCNSASSTPDSTRSARSGRSRVLEALRQSRQTVPDVNTEPSQNEQIIWDDPTAREERARLASNLQWPSCPTQYSELQVDIQNLEDSPFQKPLHDSEIAKQAVCDPGNIRVTEAPKHPISEELETPIKDSHLIPTPQAPSIAFPLANPPVAPHPREKIITIEETHEELKKQYIFQLSSLNPQERIDYCHLIEKLGGLVIEKQCFDPTCTHIVVGHPLRNEKYLASVAAGKWVLHRSYLEACRTAGHFVQEEDYEWGSSSILDVLTGINVQQRRLALAAMRWRKKIQQRQESGIVEGAFSGWKVILHVDQSREAGFKRLLQSGGAKVLPGHSVPLFKEATHLFSDLNKLKPDDSGVNIAEAAAQNVYCLRTEYIADYLMQESPPHVENYCLPEAISFIQNNKELGTGLSQKRKAPTEKNKIKRPRVH

Raw record and final code

Final code

#!/usr/bin/env python

from Bio import GenBank, Seq
from Bio.Seq import Seq,translate

# We can create a GenBank object that will parse a raw record
# This facilitates extracting specific information from the sequences
record_parser = GenBank.FeatureParser()

# NCBIDictionary is an interface to Genbank
ncbi_dict = GenBank.NCBIDictionary('nucleotide', 'genbank', parser = record_parser)

# If you pass NCBIDictionary a GenBank id, it will download that record
parsed_record = ncbi_dict['116496646']

print "GenBank id:", parsed_record.id

# Extract the sequence from the parsed_record
s = parsed_record.seq.tostring()
print "total sequence length:", len(s)

max_repeat = 9

print "method 1"
for i in range(max_repeat):
    substr = ''.join(['T' for n in range(i+1)])
    print substr, s.count(substr)

print "\nmethod 2"
for i in range(max_repeat):
    substr = ''.join(['T' for n in range(i+1)])
    count = 0
    pos = s.find(substr,0)
    while not pos == -1:
        count = count + 1
        pos = s.find(substr,pos+1)
    print substr, count

#translate sequence ==> protein, searching for start codon

start_site = s.find('ATG') #position of start codon
totrans ='' #compiling string of CDS

counter = start_site
countcodon = 0;
for i in range(len(s)-4-start_site):
    totrans = totrans + s[counter]
    counter = counter + 1
    stoptest = s[counter]+s[counter+1]+s[counter+2]
    if countcodon == 2:
        if stoptest == 'TAA' or stoptest == 'TAG' or stoptest == 'TGA':
            totrans = totrans + s[counter+1] + s[counter+2]
            break
        countcodon = -1
    countcodon = countcodon + 1
        
    
trans = translate(totrans)

print 'protein translation is %s' % trans

# new ncbi dictionary
print '\n'
new_dict = GenBank.NCBIDictionary('nucleotide','genbank')
record =  new_dict['116496646']
print record
<\pre>
Personal tools