Question: Text Histogram From Sequence Lengths Of Fasta File?
2
8.2 years ago by
United Kingdom
141341254653464453.5k wrote:

I swear I have seen this somewhere but I can't find it again. Is there a script somewhere to produce text Stem-and-Leaf Histograms from sequence lengths of fasta file?

fasta • 4.3k views
modified 3.9 years ago by Biostar ♦♦ 20 • written 8.2 years ago by 141341254653464453.5k
5
8.2 years ago by
Leonor Palmeira3.7k
Liège, Belgium
Leonor Palmeira3.7k wrote:

You can do this in R with a combination of the seqinR and the aplpack packages:

library(seqinr)
l=getLength(x)

library(aplpack)
stem.leaf(l)


There are a number of options to suit your needs:

stem.leaf(l, unit=100)


Just have a look here for more output options and details:

?stem.leaf


This functions extends far more than the simple stem() function in the graphics package.

2
8.2 years ago by
Jeroen Van Goey2.3k
Ghent, Belgium
Jeroen Van Goey2.3k wrote:

My initial idea was to use Biopython for the parsing of the fasta file and matplotlib for the drawing of the stem-leaf plot. I couldn't find a function that draws a text histogram in matplotlibto my amazement, so here is my take:

from Bio import SeqIO
from collections import defaultdict

input_file = r'path\to\file.fasta'
lengths = [len(record.seq) for record in SeqIO.parse(open(input_file), 'fasta')]
d = defaultdict(list)
for length in lengths:
d[length/10].append(length%10)

for nr in range(min(d.keys()), max(d.keys())+1):
print "{0:>3} | {1}".format(nr, ' '.join(map(str, sorted(d[nr]))) if nr in d else '')


For this fictitious input file it produces this histogram:

  5 | 3 4 8
6 | 0 0 6 7 9
7 | 0
8 | 6
9 | 0
10 | 9
11 | 2 5 6
12 | 0 0 0 0 5 8

ADD COMMENTlink modified 8.2 years ago • written 8.2 years ago by Jeroen Van Goey2.3k
1
8.2 years ago by
Kansas City

There is something in the kent source tree called textHistogram. It takes a file with one number per line, which you can get from fasta using awk...

textHistogram <(awk < test.fa '$0 !~ /^>/ {print length($0)}')