2ndscore - find the best hairpin anchored at each position.
2ndscore in.fasta > out.hairpins
For every position in the sequence this will output a line:
-0.6 52 .. 62 TTCCTAAAGGTTCCA GCG CAAAA TGC CATAAGCACCACATT
(score) (start .. end) (left context) (hairpin) (right contenxt)
For positions near the ends of the sequences, the context may be padded with 'x'
characters. If no hairpin can be found, the score will be 'None'.
Multiple fasta files can be given and multiple sequences can be in each fasta
file. The output for each sequence will be separated by a line starting with
'>' and containing the FASTA description of the sequence.
Because the hairpin scores of the plus-strand and minus-strand may differ (due
to GU binding in RNA), by default 2ndscore outputs two sets of hairpins for
every sequence: the FORWARD hairpins and the REVERSE hairpins. All the forward
hairpins are output first, and are identified by having the word 'FORWARD' at
the end of the '>' line preceding them. Similarly, the REVERSE hairpins are
listed after a '>' line ending with 'REVERSE'. If you want to search only
one or the other strand, you can use:
--no-fwd Don't print the FORWARD hairpins
--no-rvs Don't print the REVERSE hairpins
You can set the energy function used, just as with transterm with the --gc,
--au, --gu, --mm, --gap options. The --min-loop, --max-loop, and --max-len
options are also supported.
The columns for the .bag files are, in order:
1. gene_name
2. terminator_start
3. terminator_end
4. hairpin_score
5. tail_score
6. terminator_sequence
7. terminator_confidence: a combination of the hairpin and tail score that
takes into account how likely such scores are in a random sequence. This
is the main "score" for the terminator and is computed as described in
the paper.
8. APPROXIMATE_distance_from_end_of_gene: The *approximate* number of base
pairs between the end of the gene and the start of the terminator. This
is approximate in several ways: First, (and most important) TransTermHP
doesn't always use the real gene ends. Depending on the options you give
it may trim some off the ends of genes to handle terminators that
partially overlap with genes. Second, where the terminator "begins"
isn't that well defined. This field is intended only for a sanity check
(terminators reported to be the best near the ends of genes shouldn't be
_too far_ from the end of the gene).
TransTermHP uses known gene information for only 3 things: (1) tagging the
putative terminators as either "inside genes" or
"intergenic," (2) choosing the background GC-content percentage to
compute the scores, because genes often have different GC content than the
intergenic regions, and (3) producing slightly more readable output. Items (1)
and (3) are not really necessary, and (2) has no effect if your genes have
about the same GC-content as your intergenic regions.
Unfortunately, TransTermHP doesn't yet have a simple option to run without an
annotation file (either .ptt or .coords), and requires at least 2 genes to be
present. The solution is to create fake, small genes that flank each
chromosome. To do this, make a fake.coords file that contains only these two
lines:
fakegene1 1 2 chome_id
fakegene2 L-1 L chrom_id
where L is the length of the input sequence and L-1 is 1 less than the length of
the input sequence. "chrom_id" should be the word directly following
the ">" in the .fasta file containing your sequence. (If, for
example, your .fasta file began with ">seq1", then chrom_id =
seq1).
This creates a "fake" annotation with two 1-base-long genes flanking
the sequence in a tail-to-tail arrangement: --> <--. TransTermHP can
then be run with:
transterm -p expterm.dat sequence.fasta fake.coords
If the G/C content of your intergenic regions is about the same as your genes,
then this won't have too much of an effect on the scores terminators receive.
On the other hand, this use of TransTermHP hasn't been tested much at all, so
it's hard to vouch for its accuracy.
transterm(1)