Codon usage analysis for whole genomes

Question

I am new to bioinformatics. So if these questions seem you to a bit childish please forgive me.

I have two queries.

1. I am intending to perform a codon usage analysis followed by correspondence analyses for multiple microbial whole genomes of one bacterial species to find the association with the isolation source and ST type. Is this the right approach?

2. I have concatenated all the CDS in a single genome and joined all the genomes (average genes per genome is 4100). The problem is all the CDS when concatenated also have the stop codons in it. CodonW, famous for this analyses, cannot begin the analyses in the presence of stop codons. How can I remove stop codons?

Kindly suggest me a solution to the problem and concept.

Answer 1

This python script will remove all stop codons from your fasta file. I called it remove_stops.py

#!/usr/bin/env python3
import sys
from Bio import SeqIO
stop_codons = ["TAG", "TGA", "TAA", "UAG", "UGA", "UAA"]

fasta_it = SeqIO.parse(open(sys.argv[1]), 'fasta')
for fasta in fasta_it:
    name, sequence = fasta.description, str(fasta.seq)
    try:
        assert(len(sequence) % 3 == 0)
    except:
        raise IOError("length of sequence {} was not divisible by three".format(name))
    print(">{}".format(name))
    print("".join([sequence[i:i+3] for i in range(0, len(sequence) - 3 + 1, 3)
                   if sequence[i:i+3] not in stop_codons]))

You can run it like so:

python remove_stops.py stop_codons.fasta > filtered.fasta

This program won't work if there are any fasta records that do not have a length divisible by three. This indicates a sliding translation stop and these sequences should be trimmed individually. Also, all the genes must be in the 5'->3' direction.

It works by looking at each codon [sequence[i:i+3] for i in range(0, len(sequence) - 3 + 1, 3), then just makes sure they are not a stop codon: if sequence[i:i+3] not in stop_codons]