<div dir="ltr"><span style="font-size:12.8px">Hello Biopython users, </span><div style="font-size:12.8px"><br></div><div style="font-size:12.8px">I'm trying to create a script that will allow me to add 'N' to a final fasta of a template-guided assembly. I have already posted this question on BioStars and received some good feedback and with it I feel I'm 99% done with what I am trying to complete. </div><div style="font-size:12.8px"><br></div><div style="font-size:12.8px">I am using python3 and Biopython 1.66</div><div style="font-size:12.8px"><br></div><div style="font-size:12.8px">The link to the BioStars post is here: <a href="https://www.biostars.org/p/179755/" target="_blank">https://www.biostars.org/p/179755/</a></div><div style="font-size:12.8px"><br></div><div style="font-size:12.8px">The code seems to do fine on making a dictionary of the coverage file and then making a set containing only the "no coverage" areas. I can edit the fasta header information i.e. new id, new description however I am unable to make an editted version of the Seq object. I have tried using the str() function and the .tomutable() function in several ways but none of them seem to be giving me my expected output. I suspect the error is in how I'm trying to call/edit the Seq function but from the tutorial guide I'm not sure why it's not working. </div><div style="font-size:12.8px"><br></div><div style="font-size:12.8px">Any guidance would be appreciated.</div><div style="font-size:12.8px"><br></div><div style="font-size:12.8px">Thank you, </div><div style="font-size:12.8px">Sean</div><div style="font-size:12.8px"><br></div><div style="font-size:12.8px">my code is below:</div><blockquote class="gmail_quote" style="font-size:12.8px;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>from Bio.SeqRecord import SeqRecord<br>from Bio.Alphabet import IUPAC<br>from Bio.Seq import Seq<br>from Bio import SeqIO<br>import argparse <br></blockquote><div style="font-size:12.8px"> </div><blockquote class="gmail_quote" style="font-size:12.8px;margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">## this makes a dict of the samtools depth coverage input file of all the 0's<br>coverageDict = {}<br></blockquote><div style="font-size:12.8px"> </div><blockquote class="gmail_quote" style="font-size:12.8px;margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">## this loops over the input depth information and appends the dictionary with <br>## a key,value as genome position, depth of coverage<br>coverage = open("/home/sbrimer/Desktop/Columbia ST/1819_1/coverage.txt","r")<br>for line in coverage:<br>    coverages = line.split("\t")<br>    coverageDict[coverages[1]] = coverages[2]<br>coverage.close()<br></blockquote><div style="font-size:12.8px"> </div><blockquote class="gmail_quote" style="font-size:12.8px;margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">## this should read in the sequence file as an oblect, change it to a mutable<br>## object, then anywhere the coverageDict value = 0. Replace that nt with "N"<br></blockquote><div style="font-size:12.8px"> </div><blockquote class="gmail_quote" style="font-size:12.8px;margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">missing = {index for index,value in coverageDict.items() if value == 0}<br>def filter_coverage(index,value):<br>    return 'N' if index in missing else value<br></blockquote><div style="font-size:12.8px"> </div><blockquote class="gmail_quote" style="font-size:12.8px;margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">newSeqRec = []<br>with open("newfile.fasta","w") as handle:<br>    for seq_record in SeqIO.parse("Salmonella enterica subsp.fasta","fasta"):<br>        refseq = seq_record.seq<br>        newSeq = "".join(filter_coverage(index,value) for index,value in enumerate(refseq.tomutable(), start=1))<br>        newSeqRec = SeqRecord(Seq(newSeq),id ="new_id", description = str(seq_record.description))<br>    SeqIO.write(newSeqRec, handle, "fasta")<br>handle.close()</blockquote></div>