<div dir="ltr"><div>Hi Damian,</div><div>As discussed in our preliminary communications working out an acceptable solution, we have the pieces put together correctly now. Instead of following your implementation, I tried to make it more efficient by limiting it to one dictionary and them use string methods to compare the sequence to see if all match 100%.</div><div><br></div><div>I am posting the solution here in case it is of use to anyone else.</div><div>The code is posted in a gist at <a href="https://gist.github.com/fomightez/abd953068a675f416ed3">https://gist.github.com/fomightez/abd953068a675f416ed3</a> .</div><div><br></div><div>Next, we can make a version with a setting for the minimum requirement for the length of consensus block necesary to get printed to the ouput file next and post it there too.</div><div><br></div><div>The main part of my approach follows:</div><div><br></div><div>def make_dict_of_alignments(fasta_seq_file):</div><div> '''</div><div> takes a file of aligned sequence records in fasta format and reads each one</div><div> placing itin a dictionary for easy access later using the id as a key.</div><div> returns the dictionary it produced.</div><div> also returns a string of the first sequence identified to use as the measure</div><div> of length and for establishing a default consensus sequence for comparison.</div><div> '''</div><div> # initialize needed variables</div><div> aligned_seq_dict = {}</div><div> example_sequence = "NotDefined"</div><div> for seq_record in SeqIO.parse(fasta_seq_file, "fasta"):</div><div> # extract the needed info and make dictionary</div><div> name, sequence = str(<a href="http://seq_record.id">seq_record.id</a>), str(seq_record.seq)</div><div> aligned_seq_dict[name] = sequence</div><div> # grab first aligned sequence encountered as representative</div><div> if example_sequence == "NotDefined":</div><div> example_sequence = str(seq_record.seq)</div><div> return aligned_seq_dict, example_sequence</div><div><br></div><div># Need dictionary of the aligned sequences for accessing and will use</div><div># a representatice sequence from there to extract details such as length and</div><div># default consensus.</div><div>alignments_dict, representative_aligned_seq = make_dict_of_alignments(file_name)</div><div><br></div><div># initialize a string for storing the consensus. This is optional in addition to</div><div># the outout file being made</div><div>consensus_string = ""</div><div><br></div><div># Iterate over the representative sequence string using each character for</div><div># comparison. Idx holds the current index position and will make calling each</div><div># position in each sequence easy.</div><div># NOTE that the `else` is indented correctly and belongs to the `for` loop. It</div><div># will only be performed if the for loop completes without a `break` as</div><div># described at <a href="https://docs.python.org/2/tutorial/controlflow.html">https://docs.python.org/2/tutorial/controlflow.html</a></div><div>for idx, extracted_residue in enumerate(representative_aligned_seq):</div><div> for each_key in alignments_dict:</div><div> if alignments_dict[each_key][idx] != extracted_residue:</div><div> # A blank line is written to the line output file</div><div> output_file_handler.write("\n")</div><div> # consensus_string gets a blank character added to it</div><div> consensus_string += " "</div><div> break</div><div> else:</div><div> #if matches in all sequences, write the residue</div><div> output_file_handler.write(extracted_residue)</div><div> consensus_string += extracted_residue</div><div><br></div><div><br></div><div><br></div><div>Wayne</div><div class="gmail_extra"><br></div><div class="gmail_extra"><br></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Dec 15, 2015 at 5:34 AM, <span dir="ltr"><<a href="mailto:biopython-request@mailman.open-bio.org" target="_blank">biopython-request@mailman.open-bio.org</a>></span> wrote:<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
------------------------------<br>
<br>
Message: 2<br>
Date: Mon, 14 Dec 2015 16:45:22 -0900<br>
From: Damian Menning <<a href="mailto:dmenning@mail.usf.edu">dmenning@mail.usf.edu</a>><br>
To: Biopython Mailing List <<a href="mailto:biopython@mailman.open-bio.org">biopython@mailman.open-bio.org</a>><br>
Subject: [Biopython] Parsing through multiple DNA sequences one<br>
nucleotide at a time<br>
Message-ID:<br>
<CAA_KZtXdU07B1ro873i6dit4WE2q+at--Y8-WEQEsC=O6=<a href="mailto:RAeg@mail.gmail.com">RAeg@mail.gmail.com</a>><br>
Content-Type: text/plain; charset="utf-8"<br>
<br>
Hello everyone,<br>
<br>
I have multiple fasta sequences that have been aligned and cut to the<br>
same length. I want to go through each nucleotide position one at a time<br>
to see if it is the same for all of the sequences. I have written a script<br>
that can do one nucleotide position at a time but I can't figure out how to<br>
get it to loop through the entire length of the sequence.<br>
<br>
name = "test.fas"<br>
handle = open(name, 'r')<br>
conout = open("consensus.txt", 'a')<br>
<br>
val = []<br>
from Bio import SeqIO<br>
pos = 45 #test position<br>
seqlength = 1542<br>
<br>
for seq_record in SeqIO.parse(handle, "fasta"): #parses fasta file<br>
# seqlength = len(seq_record.seq)<br>
val.append(seq_record.seq[pos]) #creates dictionary with<br>
key = sequence#, value = nucleotide at position (pos)<br>
length = len(val) #determines the total<br>
number of key/value pairs<br>
<br>
y=0<br>
z=1<br>
<br>
for x in val: #parses through position<br>
values<br>
if val[y] == val[z]: #checks to see if adjacent<br>
values are equal<br>
# print val[z]<br>
z=z+1<br>
if z == length: #if all values are the<br>
same, writes value (nucleotide) to file<br>
print "position " + str(pos+1) + " equals " + val[y]<br>
conout.write(val[y])<br>
break<br>
else: #if all values are not the<br>
same, writes newline to file<br>
print "position " + str(pos+1) + " does not have a common<br>
nucleotide"<br>
conout.write('\n')<br>
break<br>
<br>
The way I have written the script, if all of the nucleotides in the same<br>
position are the same it will write the nucleotide to a file. If they are<br>
not it will write a newline. What I want is for the script to go through<br>
the length of the DNA sequence (1542 bp) and write this information to a<br>
text file so that I will end up with essentially all of the consensus<br>
sequences that I can then check for potential primer locations.<br>
<br>
When I try to put in a for or while loop I end up getting the first<br>
nucleotide repeated for each position. I think I just need to clear/reset<br>
the val dictionary after each run but it doesn't seem to work.<br>
<br>
Any help would be greatly appreciated.<br>
<br>
Damian<br>
<br>
--<br>
Damian Menning, Ph.D.<br>
<br>
"There are two types of academics. Those who use the Oxford comma, those<br>
who don't and those who should."<br><br>
</blockquote></div><br></div></div>