USSs framed

I've finished the reading frame analysis and here are the results. Each logo shows the pattern for the within-gene USSs in a specific position and orientation with respect to the reading frame of the genes they are in. The light-blue arrow below each logo shows the direction in which those USSs' genes are translated into proteins The white boxes superimposed on the arrows show the codons, with the one-letter amino acid code for those amino acids that are more- or less- strongly specified by the USS bases.

So the top left logo summarizes the 49 USSs whose positions and reverse orientation in their genes has their cores specifying the amino acids lysine (K; AAA), valine (V; GTG) and arginine (R; CGG). The consensus of the flanking regions are strong enough that they also are likely to specify particular amino acids, in this case another lysine and a phenylalanine (f; TTT). I put these last two letters in lower case to indicate that the consensus is relatively weak.

So what does this mean? We're still thinking about that. The numbers of sites (n=#, in blue) are sufficiently large that the differences between the logos are significant. I hope we will be able to correlate them to the specific coding constraints, but that's a complex analysis I'm not quite ready for yet.

Framed?

In between thinking about the sxy manuscript (yes, it's getting closer to being done), I'm making progress on the motif analysis.

One big issue arises because about 70% of the USS are in sequences that also code for proteins. So these sequences have two different evolutionary forces acting on them. First, they must continue to code for proteins that function well - any new USS that compromise the function of the part of the proteins where they arise will be eliminated by natural selection. Second, the USS that function well in DNA uptake are preferentially replacing sequences that function less well.

We already know that USS in coding regions are preferentially found in particular reading frames. The figure to the left shows the number of perfect USS cores in each of the 6 possible reading frames of their proteins, and also shows the amino acids the USS would specify in each reading frame (source Karlin et al. 1996 NAR 24:4263-4272). You can see that most of the 395 'forward' USS are in reading frame 3, where the USS core specifies the tripeptide SAV (amino acids serine alanine valine). A similar fraction of the 571 'reverse' USSs are in reading frame -2, where their cores specify TAL (threonine alanine leucine). In both orientations the favoured frame specifies a central versatile amino acid flanked by one hydrophilic and one hydrophobic amino acid. Thus USS whose cores encode SAL and TAV have been suggested to be common because they are better tolerated than the others.

This analysis considered only perfect USS cores. Now I'm using the Gibbs sampler to get better datasets, and I want to find out whether being in coding regions affects the USS motif. I've finally managed to get the sampler to work on sets of genes, by splitting the gene sequence file into four short files it could handle and by specifying the expected motif structure with a strict 'prior' file.

The TIGR database provides links where I could download a file with all the coding sequences and a fine with all the non-coding sequences. Originally I thought I could only use these to get single motifs for (i) all the USS that code in the forward direction, (ii) all the USS that code in the reverse direction, and (iii) all the USS in non-coding parts of the genome. I know that it would be much more informative if I could sort out the sequences by reading frame (as in the drawing above), but I thought that would take some sophisticated bioinformatics skills.

But today I realized that it would actually be simple. The Gibbs output specifies, for each site that fits the motif, the position of its first base relative to the first base of the sequence. Because the first base is the start of the protein, I can find the reading frame simply by dividing the position number by 3. If the position of the first base of the motif is an even multiple of three (e.g. position = 93) then the motif is in frame 1. If it's not an even multiple but the remainder is 1, the motif is in frame 2. If the remainder is 2, the motif is in frame 3.

So all I need to do is paste the list of output sites into Excel (after some massaging in Word to get rid of excess spaces and separate the columns by tabs). Then I take the column of positions, use Excel's MOD function to get the remainder after division by 3, and SORT the sequences by the value in this column.

I had originally searched for USS coding in the reverse direction by searching the forward-coding sequences for the reverse motif. This isn't very satisfactory as it gives me a reverse-orientation motif, which is hard to compare with the forward one. But now I've found a web page that will take batches of sequences and convert them into their reverse complements (the sequence of the other strand), so I can analyze these for the forward motif. (Is your head spinning yet, dear reader?) Unfortunately it can only handle about 100 sequences at once, so I have to do the whole set (about 1730) in parts.

So that's what I'm about to do. I think I have all the Gibbs output files I need for this. The outputs aren't always centered on the same USS position so they're not all equally useful, but I think I have complete sets with the same center. And I have a few more queue'd up just in case.

How random is it?

My first Gibbs sampler run was 'out of time' at 12hr, but the next one finished in 17hr. I used its output to evaluate whether reducing the stringency of the plateau period would reduce the quality of the analysis. The result: future runs will be about 40% faster because they'll use a plateau of 100 cycles rather than one of 200 cycles. I already have some whole-genome results using 100 seeds and 100 cycles, so I've queue'd up enough more to give me four with the forward-direction genome sequence and four with the reverse-complement sequence. These should provide enough data for all the whole-genome analyses.

I now realize there's another analysis I should try to do - testing the randomness of the positions of USSs around the genome. This is an interesting feature because USS spacing should reflect the forces that maintain all these USSs in the genome.

USS spacing was first addressed in the first genome USS analysis (Smith et al 1995), but they only said it was 'essentially random'. The human eye is notoriously bad at detecting randomness, so Karlin et al. took a much more rigorous approach (Karlin is a famous Stanford mathematician), calculating something called the r-scan statistic, which he developed and which looks too hard for me to follow. Karlin et al concluded that USS spacings were more even than expected for a randomly-located sequence element. This non-randomness led the authors to suggest that USSs might

contribute to global genomic activities such as replication and repair (the DNA repair hypothesis), sites of membrane attachments in association with domain loops, sites of nucleating Okazaki fragments or helix unwinding and/or sites contribution to genome packaging. (Yes, the syntax seems a bit off to me too.)

I think these suggestions are wrong, for reasons I'll go into another time, but the lack of randomness may still be telling us something important about the forces that maintain USS.

Karlin et al.'s analysis used only the positions of perfect USS cores (AAGTGCGGT and reverse complement). I think I should now repeat it on my new unbiased USS data. Well, really what I mean is that I think I should either find a tame mathematician/statistician who can show me how to do it, or find a similar analysis that's easier for me to understand. (Hmmm, I think my neighbour at a lunch on Thursday was a bioinformatics statistician - maybe she can give me some advice.)