Question: How can I calculate the number of sequences starting with different nucleotides using perl
gravatar for MAPK
2.7 years ago by
Washington University in St. Louis
MAPK1.6k wrote:

I am not familiar with perl (I mostly work with R) and I was hoping someone would help me with this problem. I have this test.fasta file below. I need to get this table from this fasta file where I have sequence length starting from minimum sequence length to maximum sequence length. Then I want to get the counts for sequences with the start bases ( first base in 5') along with the number of sequences for that given length. For example, In the given fasta below, seq1, seq5 and seq6 are 12 bases long, so in the result table below I have length 12, number of sequences of 12 bases length = 3 and the sequences starting with A= 2 and sequences starting with T =1 and so forth. How do I get this done in perl? Thank you for your help.



Result I want (shown only for the sequence of length 12bases).

length number A   T   G   C
12      3     2   1    0  0
perl fasta • 1.5k views
ADD COMMENTlink modified 2.7 years ago by Neilfws49k • written 2.7 years ago by MAPK1.6k
gravatar for lieven.sterck
2.7 years ago by
VIB, Ghent, Belgium
lieven.sterck8.7k wrote:

this script will do it:

#!/usr/bin/env perl
use strict;
use warnings;

my %counts;

while (<STDIN>){
        next if ($_ =~ /^>/ || $_ =~ /^\n/);
        my $len = length($_);
        my $fb = substr($_,0,1);
        if (!$counts{$len}{'A'}) {$counts{$len}{'A'} =0;}
        if (!$counts{$len}{'C'}) {$counts{$len}{'C'} =0;}
        if (!$counts{$len}{'G'}) {$counts{$len}{'G'} =0;}
        if (!$counts{$len}{'T'}) {$counts{$len}{'T'} =0;}


print "length\tnumber\tA\tT\tG\tC\n";
foreach my $c (sort { $a <=> $b } keys %counts){
        print "$c\t$counts{$c}{'count'}\t$counts{$c}{'A'}\t$counts{$c}{'T'}\t$counts{$c}{'G'}\t$counts{$c}{'C'}\n";


reads from STDIN outputs to STDOUT, execute as

perl < yourFile
ADD COMMENTlink modified 2.7 years ago • written 2.7 years ago by lieven.sterck8.7k

Thank you so much. How do I get the sorted values in the table (shortest to longest sequence length)?

ADD REPLYlink written 2.7 years ago by MAPK1.6k

Look into unix sort. You can sort on columns.

ADD REPLYlink written 2.7 years ago by genomax91k

I updated the script such that it now prints the list from shortest to longest.

or as genomax says : use unix sort (will allow you to sort on whatever column you want)

ADD REPLYlink modified 2.7 years ago • written 2.7 years ago by lieven.sterck8.7k
gravatar for Neilfws
2.7 years ago by
Sydney, Australia
Neilfws49k wrote:

As you mentioned that you normally work in R, I thought an R-based solution might be of interest.

Good packages for working with sequences are Bioconductor's Biostrings and the older seqinr.

Using seqinr, you can read in a fasta-formatted file like so:

myseqs <- read.fasta("test.fasta")

Then you can use functions from seqinr to get the sequence length and first character. I'd use the tidyverse packages dplyr and tidyr to manipulate the results. One issue is that both seqinr and dplyr have a count() function, so you need to specify dplyr::count() explicitly.

An example. In this case counts are not returned for C or G since no sequences in the test set start with those bases. You could fix that up later if desired, along with the NA values (zero may be preferred). I would also be inclined to leave the output in "long" format (_i.e._ don't use spread), but you may want it to be "wide" for a reason.

tibble(names = names(myseqs), 
       len = getLength(myseqs), 
       first = sapply(myseqs, function(x) getFrag.character(x, 1, 1))) %>% 
dplyr::count(len, first) %>% 
group_by(len) %>% 
mutate(total = sum(n)) %>% 
spread(first, n)

# A tibble: 4 x 4
# Groups:   len [4]
    len total     a     t
  <int> <int> <int> <int>
1    12     3     2     1
2    15     1    NA     1
3    16     1     1    NA
4    19     1    NA     1
ADD COMMENTlink written 2.7 years ago by Neilfws49k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1229 users visited in the last hour