Question: Morning Gang,

I feel like I am just missing something in the documentation for this but here goes, I'm currently running some salmonellas through the works using a very well annotated reference and we have come across a large drop out of genes in one location. We went back and confirmed they are not there with a quick pcr so I know the data is correct, however when I go to make the consensus chromosome I can not figure out how to include the drop out in the file.

For the mapping I'm using TMAP (ion torrent data, single end reads), freebayes for the haploid vcf creation, and SnpEff for the SNP annotation.

I know if I want to create a new chromosome of our sample I can use the vcf file to make the SNP changes in the reference file, however that will retain the genes that have dropped out because I'm making the changes to the reference not the alignment.

Does anyone have a suggestion or a work around for this type of problem?

Thank you in advance,
Sean

***EDIT 15MAR2016***

So I have been working on this for some time now and I'm close to my solution but I'm having an issue the "new" sequence. Below is my code.

## this makes a dict of the samtools depth coverage input file of all the 0's
coverageDict = {}

## this loops over the input depth information and appends the dictionary with
## a key,value as genome position, depth of coverage
coverage = open("/home/sbrimer/Desktop/Columbia ST/1819_1/coverage.txt","r")
for line in coverage:
coverages = line.split("\t")
coverageDict[coverages[1]] = coverages[2]
coverage.close()

## this makes a set of the large dict of just the missing.
missing = {index for index,value in coverageDict.items() if value == 0}

def filter_coverage(index,value):
return 'N' if index in missing else value

ref = open("Salmonella enterica subsp.fasta","rU")
## this should read in the sequence file as an object, change it to a mutable
## object, then anywhere the coverageDict value = 0. Replace that nt with "N"
for seq_record in SeqIO.parse(ref,"fasta"):
refSeq = seq_record.seq
new_seq = "".join(filter_coverage(index,value) for index,value in enumerate(refSeq.tomutable(), start=1))
new_rec = SeqRecord(Seq(new_seq),id="new_id",description="new_description")
SeqIO.write(new_rec,"NewSeq.fasta","fasta")

ref.close()

Everything seems to work fine until I want to add the 'N' to the fasta file, I can get the generic "new_id" and "new_description" to show up in the test file, but not the newseq. I have also tried the following:

first this

refSeq = seq_reocrd.seq
refSeq_mut = refSeq.tomutable()

then

refSeq = str(seq_record.seq)

but for whatever reason I can not seem to change the seq_record.seq part of the record.
Has anyone else had this issue?

Answer: Hi All,

It turns out the error in this scrip is right at the top. When I make the dictionary the index and value are encode as strings and not numbers so when all the other stuff runs, it runs without error and returns the original sequence because number do not match words.

Seriously. That took hours for me to figure out.

Lesson learned.

Please find the corrected script below if for some reason anyone else needs to do something like this.

All that is left to do with it is ti set up the argparser for the options I want so it will work with any file I pass.

from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC
from Bio.Seq import Seq
from Bio import SeqIO
import argparse

## this makes a dict of the samtools depth coverage input file of all the 0's
coverageDict = {} ## first the empty dictionary

## this loops over the input depth information and appends the dictionary with
## a key,value as genome position, depth of coverage
coverage = open("/home/sbrimer/Desktop/Columbia ST/1819_1/coverage.txt","r")
for line in coverage:
coverages = line.strip().split('\t') #added strip, values also had newline character
coverageDict[int(coverages[1])] = int(coverages[2])
coverage.close()

## making a set of only the index of the missing
missing = {index for index,value in coverageDict.items() if value == 0}

def filter_coverage(index,value):
if index in missing:
return "N"
else:
return value

## this should read in the sequence file as an oblect, change it to a mutable
## object, then anywhere the coverageDict value = 0. Replace that nt with "N"
append_handle = open("NewFile.fasta","a")
for seq_record in SeqIO.parse("Salmonella enterica subsp.fasta","fasta"):
newSeq= "".join(filter_coverage(index,value)for index,value in enumerate(seq_record.seq, start= 1))
SeqIO.write(SeqRecord(Seq(newSeq),id="something",description="something_else"),append_handle,"fasta")
append_handle.close()