<p>Don't use the unique sequences. The whole point of sequence logos is that the height of the letter represents its frequency at a given position. You remove this information when you make your sequences unique. So, while it makes sense to use unique sequences when aligning, you shouldn't when building your logo.</p>

<p>Generate a simple mapping file (<code>counts</code>) with the number of times each sequence was found:</p>

<pre><code>seq1 3
seq2 10
seq3 2
</code></pre>

<p>Then, given an aligned protein fasta file (<code>foo.aln</code>) like this:</p>

<pre><code>&gt;seq1
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
&gt;seq2
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
&gt;seq3
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL
</code></pre>

<p>You can simply repeat the sequences as needed (see my answer <a href="https://bioinformatics.stackexchange.com/a/2651/298">here</a> for the <code>FastaToTbl</code> and <code>TblToFasta</code> commands):</p>

<pre><code>$ awk '
        (NR==FNR){
            count[$1]=$2; 
            next
        } 
        {
            for(k=1;k&lt;=count[$1];k++){
                print $1"."k"\t"$2
            }
        }' counts &lt;(FastaToTbl foo.aln) | TblToFasta
&gt;seq1.1 
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
&gt;seq1.2 
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
&gt;seq1.3 
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
&gt;seq2.1 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
&gt;seq2.2 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
&gt;seq2.3 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
&gt;seq2.4 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
&gt;seq2.5 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
&gt;seq2.6 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
&gt;seq2.7 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
&gt;seq2.8 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
&gt;seq2.9 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
&gt;seq2.10 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
&gt;seq3.1 
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL
&gt;seq3.2 
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL
</code></pre>

<p>Note that the above assumes your shell supports the <code>&lt;(command)</code> syntax for <a href="http://www.tldp.org/LDP/abs/html/process-sub.html" rel="nofollow noreferrer">process substitution</a>. Bash, the default shell on most Linux and macOSX installations does, but if you are using a shell that doesn't, you might need to run <code>FastaToTbl foo.aln &gt; foo.tbl</code> first, and then use <code>foo.tbl</code> as input.  </p>

CC BY-SA 3.0
2018-01-29T12:27:11.010Z
That's great thanks. I have my tabular file containing the read counts and my aligned fasta file although I am only familar with using R not c++. Can this script be adapted/used in R?
clfrankling
CC BY-SA 3.0
2018-01-29T12:28:59.970Z
That isn't C++! It's just a shell command using [`awk`](https://en.wikipedia.org/wiki/AWK). Just copy/paste the command above directly into a *nix terminal. You will have to also take the two scripts from the other answer I link to.
terdon
CC BY-SA 3.0
2018-01-31T11:46:37.767Z
Ok sorry I did not recognise this. I have my two files: sequence.tbl (ID with aligned sequences already converted to tabular form) and counts.tbl (ID with read counts). I set up the directory to these files and entered your code into my ubuntu terminal and nothing happens. I doesn't print the output or come up with any errors? Do I need to add a command to print into an output tabular file?
clfrankling
CC BY-SA 3.0
2018-01-31T11:54:12.773Z
I copy/paste the exact command above just changing the last line to, 
        }' counts.tbl < sequence.tbl
clfrankling
CC BY-SA 3.0
2018-01-31T16:49:27.707Z
@clfrankling yeah, that's a big change :) The `<` is wrong. If you already have them in tbl format, you need the command to end with `}' counts.tbl  sequence.tbl`.  Ping me (`@terdon`) in [our chat room](https://chat.stackexchange.com/rooms/58858/bioinformatics) if you need more help. It will be easier to debug that way
terdon
CC BY-SA 3.0
2018-01-24T11:43:41.427Z
2018-01-24T11:43:41.427Z
terdon
<p>I have a large aligned amino acid sequences file (fasta file) from next generation sequencing data. 
The unique sequences have been aligned but these all have certain read counts (frequencies associated with them).</p>

<p>I would like to generate a weblogo image that has the true probability of amino acids at each position. 
Is there a way to input frequencies/read counts with the aligned files of sequences using weblogo?
For example, in the title of the fasta file or as a csv/tsv file?</p>

<p>So that when generating the weblogo image it takes into account the abundances of each sequence?</p>

CC BY-SA 3.0
2018-01-24T10:50:06.167Z
can't you just duplicate the sequences to match the original frequencies?
heathobrien
<phylogenetics>
How to generate a weblogo image that takes into account sequence frequencies?
# How to generate a weblogo image that takes into account sequence frequencies?

(Asked by: username_0 on 2018-01-24)

I have a large aligned amino acid sequences file (fasta file) from next generation sequencing data.
The unique sequences have been aligned but these all have certain read counts (frequencies associated with them).

I would like to generate a weblogo image that has the true probability of amino acids at each position.
Is there a way to input frequencies/read counts with the aligned files of sequences using weblogo?
For example, in the title of the fasta file or as a csv/tsv file?

So that when generating the weblogo image it takes into account the abundances of each sequence?

**Tags:** phylogenetics
> **Comments:**
> **username_1**: can't you just duplicate the sequences to match the original frequencies?  

---

## 1 answers

### Answer by username_2 on 2018-01-24. Score: 3 | Accepted answer

Don't use the unique sequences. The whole point of sequence logos is that the height of the letter represents its frequency at a given position. You remove this information when you make your sequences unique. So, while it makes sense to use unique sequences when aligning, you shouldn't when building your logo.
Generate a simple mapping file (counts) with the number of times each sequence was found:
seq1 3
seq2 10
seq3 2

Then, given an aligned protein fasta file (foo.aln) like this:
>seq1
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq2
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq3
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL

You can simply repeat the sequences as needed (see my answer here for the FastaToTbl and TblToFasta commands):
$ awk '
        (NR==FNR){
            count[$1]=$2; 
            next
        } 
        {
            for(k=1;k<=count[$1];k++){
                print $1"."k"\t"$2
            }
        }' counts <(FastaToTbl foo.aln) | TblToFasta
>seq1.1 
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq1.2 
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq1.3 
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq2.1 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.2 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.3 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.4 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.5 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.6 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.7 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.8 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.9 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.10 
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq3.1 
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL
>seq3.2 
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL

Note that the above assumes your shell supports the <(command) syntax for process substitution. Bash, the default shell on most Linux and macOSX installations does, but if you are using a shell that doesn't, you might need to run FastaToTbl foo.aln > foo.tbl first, and then use foo.tbl as input.

> **Comments:**
> **username_0**: That's great thanks. I have my tabular file containing the read counts and my aligned fasta file although I am only familar with using R not c++. Can this script be adapted/used in R?  
> **username_2**: That isn't C++! It's just a shell command using [`awk`](https://en.wikipedia.org/wiki/AWK). Just copy/paste the command above directly into a *nix terminal. You will have to also take the two scripts from the other answer I link to.  
> **username_0**: Ok sorry I did not recognise this. I have my two files: sequence.tbl (ID with aligned sequences already converted to tabular form) and counts.tbl (ID with read counts). I set up the directory to these files and entered your code into my ubuntu terminal and nothing happens. I doesn't print the output or come up with any errors? Do I need to add a command to print into an output tabular file?  
> **username_0**: I copy/paste the exact command above just changing the last line to,          }' counts.tbl < sequence.tbl  
> **username_2**: @clfrankling yeah, that's a big change :) The `<` is wrong. If you already have them in tbl format, you need the command to end with `}' counts.tbl  sequence.tbl`.  Ping me (`@terdon`) in [our chat room](https://chat.stackexchange.com/rooms/58858/bioinformatics) if you need more help. It will be easier to debug that way  

---

