How To Calculate The Identify Of Each Site For A Sequence Alignment?

Question

[biolab]

Dear all, I have an alignment of several sequences. I need ot make statistics ot the identify of each nt site. Would you please give a hand? Thank you very much.

[samsara]

Can you make clear what '..calculate the identify ..' means ?

[biolab]

Sorry for unclear descirption. Actually my analysis is to calculate the frequency of each site for an alignment. For example, 3 sequences have been aligned as shown below. I want to know the frequency of these four sites (No.1: 100% A, NO.2: 100% T, NO3: 100%G, NO. 4: 33.3%C 66.6%G). THANKS!

>seq1
ATGC
>seq2
ATGG
>seq3
ATGG

[biolab]

i further make some comments: my sequences are the same in length. So I just need to get the nt composition information for each position. THANKS.

[samsara]

do you allow gaps in your alignment ?

[biolab]

no gaps in the alignment. Actually my sequences are all 21 nt miRNAs.

[samsara]

Then it would be easy. Put all sequences in array, make loop with increment counter by array indices, then calculate similarity percentage in each loop count

[biolab]

Hi World, I have encounted problem with perl. I am a true beginner of perl. I tried my best to write a script, but not working. I show it below. This is an example: I have one alignment file that contains 4 sequences, each of which is 21 nt long. Actually I also need to run more alignment files in batch, but do not know how to use @ARGV. Could you please take a look at the following script? THANKS!

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

open SEQ, "seqaln.txt"; #one file only, need to batch input many files
my @seq=<SEQ>;

    my $j = 0; 
    while ($j < 21){
    my $o = substr($seq[2],$j,1);
    my $p = substr($seq[4],$j,1);
    my $q = substr($seq[6],$j,1);
    my $r = substr($seq[8],$j,1);
    $j++;
print "$o\t$p\t$q\t$r";
}

close SEQ;

[biolab]

Note I don't know how to calculate similarity percentage, so just print the four nts. Please modify my script. THANKS

[biolab]

I finally worked out a script for displaying the nt at each position for 4 sequences. for example

>s1
ATGCACTGACTGTATGACTG
>S2
ATGGTGACTGTGACTGACTG
>S3
ATGGACCATGACTGCATGTG
>S4
ATCCACTGTGACGTGCAACA

My script is as follows:

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

open SEQ, "seqaln.txt"; 
my @seq=<SEQ>;

    my $j = 0; 
    while ($j < length($seq[1])){
    my $o = substr($seq[1],$j,1);
    my $p = substr($seq[3],$j,1);
    my $q = substr($seq[5],$j,1);
    my $r = substr($seq[7],$j,1);
    $j ++;
print "$j\t$o\t$p\t$q\t$r\n";
}

close SEQ;

my output result is:

1       A       A       A       A
2       T       T       T       T
3       G       G       G       C
4       C       G       G       C
5       A       T       A       A
6       C       G       C       C
7       T       A       C       T
8       G       C       A       G
9       A       T       T       T
10      C       G       G       G
11      T       T       A       A
12      G       G       C       C
13      T       A       T       G
14      A       C       G       T
15      T       T       C       G
16      G       G       A       C
17      A       A       T       A
18      C       C       G       A
19      T       T       T       C
20      G       G       G       A
21
22

As a beginner i tried best to produce this output. However, I need to calculate in each position (here from position 1 to 20) the frequency of ATGC. How to achiev that? Could anyone give a hand? THANKS!

Answer 1

[samsara]

Here is python script. I think this is what you want.

from numpy import matrix
from collections import defaultdict

Seq1="ATCCTTAGGT"
Seq2="ACCCTACAAT"
Seq3="TTCTGGCAGC"

SequenceLength=len(Seq1)

A = matrix( [list(Seq1),list(Seq2),list(Seq3)])   

for i in range(0,SequenceLength):
    dic = defaultdict(int)
    for j in range(0,len(A)):
        dic[A[j,i]] += 1
    Str="Position:"+str(i+1)
    if dic['A'] !=0:
        Str = Str + " A: " + str(float(dic['A'])/len(A) *100)

    if dic['T'] !=0:
        Str = Str + " T: " + str(float(dic['T'])/len(A) *100 )

    if dic['G'] !=0:
        Str = Str + " G: " + str(float(dic['G'])/len(A) *100)

    if dic['C'] !=0:
        Str = Str + " C: "+ str( float(dic['C'])/len(A) *100)

    print Str

Output:

Position:1 A: 66.6666666667 T: 33.3333333333
Position:2 T: 66.6666666667 C: 33.3333333333
Position:3 C: 100.0
Position:4 T: 33.3333333333 C: 66.6666666667
Position:5 T: 66.6666666667 G: 33.3333333333
Position:6 A: 33.3333333333 T: 33.3333333333 G: 33.3333333333
Position:7 A: 33.3333333333 C: 66.6666666667
Position:8 A: 66.6666666667 G: 33.3333333333
Position:9 A: 33.3333333333 G: 66.6666666667
Position:10 T: 66.6666666667 C: 33.3333333333

[biolab]

THANKS a lot for the python script!

[Andreas]

At least upvote if you like the answer

[biolab]

I have upvoted these answers, which are very helpful! (Sorry, i am new in Biostar and have not noticed upvote). THANKS again for the python and perl scripts!

Answer 2

[PoGibas]

biolab, know what you're looking for...
This "statistics ot the identify of each nt site" is called Position weight matrix (or PSWM, PSSM).
I strongly recommend going trough this tutorial (you can analyze your sequences their).
Also, similar question were asked on biostars before: Parsing a Fasta file by columns.

Answer 3

[swbarnes2]

Your script is close. What you want is something that keeps count of how often every letter comes up in each position. I'd use a hash of hashes. It's sort of like a three dimensional hash, where instead of one key having one value, each key is its own hash, with more keys and values. So

$hash{14}{'A'}++

would increment up the count of A's at position 14 by 1. $hash{14} is actually a reference to another hash that has the counts for all 4 letters. So:

$hash{$j}{$o}++

and so on for each sequence.

Then you use two nested foreach loops to get the data out

foreach my $position (keys %hash) {
     foreach my $letter ('A', 'C', 'G', 'T') {
          print "$position: $letter: ", $hash{$position}{$letter} / 4, "\n";
     }
}

Once you get your script working, you next job is to generalize it, so it works on more than 4 sequences.

[biolab]

THANKS for your suggestions, really helpful!

[biolab]

Dear swbarnes2,

I tried my best to write a code, as shown below. However, I still get problem. could you please check it and indicate the erros for me? THANK YOU VERY MUCH.

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

open SEQ, $ARGV[0];
my %hash = 0;
while (<SEQ>){
chomp;
my @seq = split (/ /, $_);
    for (my $i=0; $i<=$#seq; $i++){
    $hash {$i}{$seq[$i]} +=1;
    }
    }
close SEQ;
    foreach my $posi (keys %hash){
        foreach my $nt ('A', 'T', 'G', 'C'){
    print "$posi", "$nt", %hash{$posi}{$nt}/4, "\n";
        }
    }

Error message

syntax error at aln comp_2.pl line 18, near "%hash{"
Global symbol "$nt" requires explicit package name at aln comp_2.pl line 18.
syntax error at aln comp_2.pl line 18, near "},"
syntax error at aln comp_2.pl line 20, near "}"
Execution of aln comp_2.pl aborted due to compilation errors.