4
$\begingroup$

I'm given a Fasta file, containing a large DNA(over 115,000 long) sequence, and I am tasked with finding a single large open reading frame contained within the DNA sequence using Biopython.

I'm aware this has been asked before consistently however i can't find anything for actually finding the source of my ORF, but apologies for any overlap.

Now from other sources and the Biopython cookbook I've translated my sequence and found six open reading frames (three for each strand) and their positions within the sequence;

def find_orfs_with_trans(seq, trans_table, min_protein_length):
    answer = []
    seq_len = len(seq)
    for strand, nuc in [(+1, seq), (-1, seq.reverse_complement())]:
        for frame in range(3):
            trans = str(nuc[frame:].translate(trans_table))
            trans_len = len(trans)
            aa_start = 0
            aa_end = 0
            while aa_start < trans_len:
                aa_end = trans.find("*", aa_start)
                if aa_end == -1:
                    aa_end = trans_len
                if aa_end-aa_start >= min_protein_length:
                    if strand == 1:
                        start = frame+aa_start*3
                        end = min(seq_len,frame+aa_end*3+3)
                    else:
                        start = seq_len-frame-aa_end*3-3
                        end = seq_len-frame-aa_start*3
                    answer.append((start, end, strand,
                                   trans[aa_start:aa_end]))
                aa_start = aa_end+1
    answer.sort()
    return answer

orf_list = find_orfs_with_trans(record.seq, table, min_pro_len)

for start, end, strand, pro in orf_list:
print("%s...%s - length %i, strand %i, %i:%i" \
      % (pro[:30], pro[-3:], len(pro), strand, start, end))

OUTPUT

 TNRQVYGGTLQSLRTGTGIYSRLASSPTNR...CEA - length 160, strand 1, 7950:8433
 GIPPGRTEGLGRYVHGESIFLEATLVPEPQ...VQS - length 171, strand -1, 19275:19791
 ISVHLRRYFSVLLRAPVALNADSRVAGPLD...ITP - length 190, strand 1, 34079:34652
 GCVNGNFPDYRVAVDDPGALVVGGELGQAL...VNT - length 771, strand -1, 39335:41651
 LLKLLASRLTWLLVLSAALGFLTSVSYRLG...SPG - length 235, strand -1, 39358:40066
 PGDVIVSKVPADNLRPRMSEINDFTNSIII...VLV - length 826, strand 1, 39362:41843
 GPGRQSSSAYERDQRLHQQHHHQRERHQVR...QRS - length 764, strand 1, 39385:41680
 LLLSSTSSASCWLKKFFFRLLLLSRSVALL...LAL - length 158, strand -1, 40459:40936
 TRLIIRPPSHTTLYYDPYPQPLKVSSSVPC...RIT - length 159, strand -1, 54663:55143
 IGAFSRGYNQGWVLIEAQEVRTTSPLFPRS...RDQ - length 201, strand 1, 56104:56710
 KARGSSLGAIYCRRLKRRSLPSWQIESFNC...TNS - length 150, strand -1, 62497:62950
 RRPARAGAELNKDDIRDTCLLSCSEVDNKV...VTL - length 168, strand 1, 72339:72846
 SGPTTVRTSAPQRCWRTSLKPKLSLSDSTE...PVA - length 177, strand -1, 114494:115028
 ERELWLETGPPTPLWGTGPDCSRAALNRQR...PRR - length 157, strand 1, 114953:115427

However I'm unsure if these open reading frames contain the ribosome binding site(consensus sequence of AAGGAGGTG and occurring between 6 to 9 nucleotides upstream from the positive 5' strand). Is there anyway to add to my code to include this? or would it already be included with Biopython?

My specification says to find the single open reading frame however since none of these begin the the M start codon i'm a bit lost how they can be so, I have no information on where my sequence is from (something I have to find out) so using the standard translation table 1.

$\endgroup$
0
2
$\begingroup$

From what I can tell, the program output does not indicate whether the open reading frames (ORFs) are preceded by ribosome binding sites (RBSs). I'm not sure what parameters are used when establishing ORFs. To verify your RBS requirement, index the input DNA (using given ORF locations from biopython output, pay attention to strands, which are designated by 1/-1) and search upstream for a sequence similar to your consensus. You'd have to pre-define your criteria for similarity to consensus in order to find it programmatically (or you can do it visually to better avoid false negatives). Perhaps try a search that permits ambiguous/non-matching third nucleotides in codons (like a regular expression such as "AA.GA.GT." where "." can be any nucleotide). Biologically, third nucleotide changes are least probable to drastically alter the translated residue (more likely to encode amino acid with similar chemical properties).

$\endgroup$
2
$\begingroup$

You should start by finding all instances of AAGGAGGTG and then look for ORFs that begin 6 to 9 nucleotides downstream of them. By begin, I mean find an AUG codon, and then extend until the next in-frame stop codon.

There's also no point in translating, just look for long stretches of sequence between AUG and UGA, UAA or UAG, those will define your reading frames. You seem to be starting the reading frame at any non-stop codon instead of starting from AUGs. This may well make biological sense in some contexts (there are proteins that don't start with M) but is probably beyond the scope of what sounds like a simple assignment.

$\endgroup$
2
  • $\begingroup$ Thanks for the heads up, yeah I think I can safely assume that yeah the start codon will be an M. Apologies for the follow up question, is there anywhere you can point me in the right direction to in terms of achieving this? I'm struggling to find anything related. $\endgroup$
    – daenwaels
    Dec 13 '17 at 19:28
  • $\begingroup$ @maltedg sorry, I just saw your comment, but I don't know what you mean. This is pointing you in the right direction. As I said, start by finding AAGGAGGTG, and then look for ORFs that start after it. If you can assume that they'll always start with an AUG, then look for AUGs close to the AAGGAGGTG and extend the reading frame from there. $\endgroup$
    – terdon
    Dec 18 '17 at 13:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.