# determine longest poly-purine tract

how to determine/ find the longest poly-purine tract in any genome and this needs to be on the E. coli genome . is it to figure out the polypurine tract and then figure out the longest chain ? or is it to splice the introns and exons away from the DNA ? since E. coli's genome is 4.6 million BP long, i need some help in breaking this down ?

• One comment: Bacteria have no intron/exons, only eukaryotes have them. Otherwise I don't understand your question. For what reason are you looking for the "longest poly-purine"? – Chris Aug 8 '14 at 20:57
• i need to find out a longest poly-purine tract in a E.coli genome – harry Aug 8 '14 at 21:27
• Is a sequenced genome on the internet good enough, or do you need to take real cells, extract the DNA, and analyze it directly for polypurines? – user137 Aug 8 '14 at 21:40
• yes sequenced genome on the internet is good enough – harry Aug 8 '14 at 21:42
• are you on windows, mac, or linux? Please say linux. – user137 Aug 8 '14 at 23:27

Code Golf!

This is my Python script, avoiding regex and not at all prettified.

def gettract(strand):
#if strand is 1 then replace C with T and split on T
#else replace A with G and split on G

if strand==1:
print 'top strand purine tracts:'
aseq=seq.replace('C','T')
tracts=aseq.split('T')
elif strand==2:
print 'top strand pyrimidine tracts:'
aseq=seq.replace('A','G')
tracts=aseq.split('G')
sortedtracts=sorted(tracts,key=len,reverse=True)
for i in range(5):
print len(sortedtracts[i]), sortedtracts[i]

f=open('ec_dna.txt','r')
print 'Sequence analysed is %s' % l[1:] #the header line
f.close()
seq=''.join([i.rstrip() for i in l[1:]])
print 'sequence length is %i' %len(seq)
print
gettract(1)
print
gettract(2)


And here is the output:

Sequence analysed is E. coli K-12 MG1655 U00096.2 (1 to 4639675 = 4639675 bp)

sequence length is 4639675

top strand purine tracts:

28 AGAGAGAAAAGAAGAAGAGAAAAGAAAA

22 GGGGGAGAAAGGGGGGAAAAGG

20 AAAGAGGAAAGGGGGGGGGG

20 GAAGAGGAAAAGAAAGAAGG

20 AAAAAAGAAGAGAAAAAAGA

top strand pyrimidine tracts:

29 TTTTTCTCTTTCCCTTTTCCCTTCCCTCT

22 CCCTCCTCTTCTCTCCCTCCCC

21 TTCTCCTCCTTCTCTTTTCTT

21 TTTTTCTTCTCTCTTTTCCTC

21 CTCCCCTTTTCCCTTTTCTCC

The longest purine tract is 29 bases long, and is on the bottom strand.

This script doesn't handle the problem that the E. coli genome is a circular molecule, so if the longest tract was split at the origin it would be missed - I leave that and obtaining the coordinates of the longest tract as an exercise.

Why do you need to do this?

This is pretty simple with the right computer tools and some knowledge of how to use regular expressions, but the fact that you are daunted by a 4.6 Mbase file suggests to me you don't have that.

Using an old copy of E.coli DH10B, and Lasergene's SeqBuilder, I found a stretch of 29 Cs and T's, so in the reverse stand, that would be a poly-purine. I don't see a longer stretch in my fasta.

• thank you for the reply . can you share your code or steps followed ? iam trying to use regular expressions . what are the recommend tools needed ? – harry Aug 8 '14 at 22:16
• Can I get the DNA sequence you used? I'd like to test it against my code. – user137 Aug 9 '14 at 0:50
• I used commercial Windows software, just to get the answer fast. If I was going to do it in Perl, I'd make sure newlines were stripped out, and just search for /[AG/{30}]/ and if I didn't find it, I'd keep dropping the number until it hit. And the same for CT strings. You might not be able to search the whole 4.6 million letter string at once, but cutting it into overlapping first 10k long pieces should work. Grep in Unix might be able to handle the whole string at once. – swbarnes2 Aug 9 '14 at 19:29
• I meant /[AG]{30}/ or /[CT]{30}/ – swbarnes2 Aug 9 '14 at 20:09
• search for /[AG/{30}]/ and if I didn't find it, I'd keep dropping the number until it hit... I like this method.. It is likely to be fast. But the issue is that it wont be able to find all instances- but a small modification in the code would fix it. Keep removing the part rightwards from the first match position and continue regex search. An awk example here – WYSIWYG Aug 11 '14 at 6:26

Since you said an internet sequence is good enough, I will propose a computer programming approach. The basic idea is to enter the sequence as a string of text and then have the program read the text, keeping track of the length of the current polypurine tract, and the longest encountered so far. When it's done, it reports the longest one. I'm not sure if you know any programming, so I'll attempt to write up an example in C++ and post it. There are probably other programs that can do this, but I don't know about them. Here is the code:

## main.cpp

#include <fstream>
#include <iostream>

// this program reads a DNA sequence from a text file and then looks for
// the longest polypurine sequence. Please make sure there isn't non DNA
// text in the text file. Please call your text file "DNA.txt" and place
// in the same folder as this program.

std::string base_complement(std::string in){
std::string return_me = "";
for(int i = 0; i < in.size(); i++){
if(in.c_str()[i] == 'a' || in.c_str()[i] == 'a'){
return_me += 'T';
}
if(in.c_str()[i] == 'g' || in.c_str()[i] == 'G'){
return_me += 'C';
}
if(in.c_str()[i] == 't' || in.c_str()[i] == 'T'){
return_me += 'A';
}
if(in.c_str()[i] == 'c' || in.c_str()[i] == 'C'){
return_me += 'G';
}
}
return return_me;
}

int main(){

std::cout << "WELCOME TO POLYPURINE COUNTER" << std::endl;
std::cout << "WRITTEN BY SAM CROWLEY, SAMUEL-CROWLEY@UIOWA.EDU" << std::endl;
std::cout << "WILL READ THE DNA IN DNA.txt AND REPORT LONGEST PURINE STRETCH" << std::endl;
std::cout << std::endl;
// open the file and read the contents into the string
std::string DNA;    // this string will hold the entire DNA sequence
std::string line;   // this will hold one line of the file at a time
std::ifstream text_file("DNA.txt");

while(!text_file.eof()){        // read the file until we get to the end
getline(text_file, line);   // get one line at a time in the line string
DNA += line;                // add that line to the DNA string
}

//std::cout << DNA << std::endl;  // output the DNA string just to check

int current_tract = 0;          // length of current tract
int longest_tract = 0;          // length of longest tract so far
// true = sense strand, false = opposite strand
bool tract = (DNA.c_str() == 'a' || DNA.c_str() == 'g');
char c = DNA.c_str();
std::string purine = "";
std::string pyrimidine = "";
std::string longest = "";
if(c == 'a' || c == 'g' || c == 'A' ||c == 'G'){
purine += c;
}else{
pyrimidine += c;
}
for(int i = 1; i < DNA.size(); i++){    // iterate through the DNA string
if(c == 'a' || c == 'g' || c == 'A' || c == 'G'){
if(DNA.c_str()[i] == 'a' || DNA.c_str()[i] == 'g' || DNA.c_str()[i] == 'A' || DNA.c_str()[i] == 'G'){
purine += DNA.c_str()[i];
}else{
std::cout << purine << std::endl;
if(purine.size() > longest.size()){
longest = purine;
}
purine = "";
c = DNA.c_str()[i];
pyrimidine += c;
}
}else{
if(DNA.c_str()[i] == 't' || DNA.c_str()[i] == 'c' || DNA.c_str()[i] == 'T' || DNA.c_str()[i] == 'C'){
pyrimidine += DNA.c_str()[i];
}else{
std::cout << pyrimidine << std::endl;
if(pyrimidine.size() > longest.size()){
longest = pyrimidine;
}
pyrimidine = "";
c = DNA.c_str()[i];
purine += c;
}
}
}
if(purine.size() > 0){
std::cout << purine << std::endl;
if(purine.size() > longest.size()){
longest = purine;
}
}
if(pyrimidine.size() > 0){
std::cout << pyrimidine << std::endl;
if(pyrimidine.size() > longest.size()){
longest = pyrimidine;
longest = base_complement(longest);
}
}
longest_tract = longest.size();
std::cout << std::endl;
std::cout << "LONGEST POLYPURINE STRETCH IS " << longest_tract << " BASES LONG." << std::endl;
std::cout << "LONGEST STRETCH IS " << longest << std::endl;
}


If you're running linux put the file in a folder and use a terminal to go that folder, then run g++ main.cpp -o polypurine.exe Then run it with ./polypurine.exe

The DNA sequence goes in a file called DNA.txt in the same folder. Make sure it only contains a, t, g, c, so remove any header data or other annotations. If it contains an X, N, or other symbol to represent ambiguous bases I don't know what will happen, but the results may not be accurate. I forget which is used for purines and which for pyrimidines, but if you tell me I can modify the code to make it work with some ambiguous bases.

If you use windows or mac this same code should work but compiling it into a working program will be different. I don't program on those machines so I don't know how to make this work on them.

Also, I ran a test on the cDNA for bee venom PLA2 and got a longest purine tract of 17 bases. I think I can post that DNA here:

## Bee Venom PLA2 cDNA

atgattatctatccgggtacactgtggtgtggtcatggtaataaaagcagcggtccgaacgaactgggtcgttttaaacataccgatgcatgttgtcgtacccatgatatgtgtccggatgttatgagtgccggtgaaagcaaacatggtctgaccaataccgcaagccatacccgtctgagctgtgattgtgatgataaattctatgattgcctgaaaaacagcgcagataccatcagcagctattttgttggcaaaatgtacttcaacctgatcgacaccaaatgctataaactggaacatccggttaccggttgtggtgaacgtaccgaaggtcgttgtctgcattataccgttgataaaagcaaaccgaaagtgtatcagtggtttgatctgcgcaaatatggatccatggataaagttttccggaattccgcaaaaaagaagagaaaggtagaagaccccaaggactttccttcagaattgctaagttttttgagtccaagcgcggccgca


If you use this program for publications or anything, let me know. I put it here so you can use it, but I'd at least like to know if it's used.

Just search for /[AG/{30}]/ and if I didn't find it, I'd keep dropping the number until it hit - swbarnes2

Code for swbarnes2's method (will work on Linux terminal, Mac Console and cygwin terminal)

Usually in a genome file a single fasta sequence will represent one chromosome. In E.coli there will be only one. You would need to pop the fasta header (first line) out when analyzing your sequence. (Assuming that you have done so).

This code below will find all polypurine tracts of a given length in the genome. Basically you have to keep dropping the part from the sequence start to the match location and redo the regex search.

cat ecoli_genome.fa | tr -d "\n" | awk '{for(i=30;i>=10;i--){
pp="[A|G]{"i"}"
s=\$0
while (match(s,pp)){
print substr(s, RSTART, RLENGTH)
s=substr(s, RSTART+RLENGTH)
}
exit
}
}'


Regex searches of type /pattern{N}/ may not work on all versions of awk. It does work in gawk >=4.0 but doesn't work in mawk.

This will load the entire genome sequence in the RAM (this is not huge at all). If you have RAM limitations and the genome file is really huge then you can stream read letter by letter (but you can't do regex search in this case). Stream read and regex are two strong features of awk.