Answers.Body: <p>For anyone interested in doing this from any sort of alignment file, I've implemented a <a href="https://github.com/biologyguy/BuddySuite/wiki/AB-Position-Frequency-Matrix" rel="nofollow noreferrer">position frequency matrix function in AlignBuddy</a>.</p>
  
  <p>input:</p>
  
  <pre><code> 15 9
  seq_1      CAGGTAGCC
  seq_2      CCGGTCAGA
  seq_3      AGGGTTTGA
  seq_4      TTGGTGAGG
  seq_5      CAAGTATGA
  seq_6      ACTGTATGC
  seq_7      CTGGTAACC
  seq_8      TATGTACTG
  seq_9      GCTGTGAGA
  seq_10     CAGGTGGGC
  seq_11     TCAGTGAGA
  seq_12     GGGGTGAGT
  seq_13     TGGGTATGT
  seq_14     GAGGTGAGA
  seq_15     CAGGTGGAG
  </code></pre>
  
  <p>Command:</p>
  
  <pre><code>:$ alignbuddy input.phy --pos_freq_mat
  </code></pre>
  
  <p>Output:</p>
  
  <pre><code>### Alignment 1 ###
  A       0.133   0.400   0.133   0.000   0.000   0.400   0.467   0.067   0.400
  C       0.400   0.267   0.000   0.000   0.000   0.067   0.067   0.133   0.267
  G       0.200   0.200   0.667   1.000   0.000   0.467   0.200   0.733   0.200
  T       0.267   0.133   0.200   0.000   1.000   0.067   0.267   0.067   0.133
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-16T17:08:29.763Z
Answers.LastActivityDate: 2017-12-21T15:15:59.897Z
Answers.LastEditDate: 2017-12-21T15:15:59.897Z
Answers.UserDisplayName: Steve Bond
Answers.Body: <p>awk</p>
  
  <pre><code>awk '{L=length($1);for(i=1;i&lt;=L;i++) {B=substr($1,i,1);T[i][B]++;}} END{for(BI=0;BI&lt;4;BI++) {B=(BI==0?"A":(BI==1?"C":(BI==2?"G":"T")));printf("%s",B); for(i in T) {tot=0.0;for(B2 in T[i]){tot+=T[i][B2];}printf("\t%0.2f",(T[i][B]/tot));} printf("\n");}}' input.txt
  
  A   0.13    0.40    0.13    0.00    0.00    0.40    0.47    0.07    0.40
  C   0.40    0.27    0.00    0.00    0.00    0.07    0.07    0.13    0.27
  G   0.20    0.20    0.67    1.00    0.00    0.47    0.20    0.73    0.20
  T   0.27    0.13    0.20    0.00    1.00    0.07    0.27    0.07    0.13
  </code></pre>
  
  <p>Or, expanded for clarity:</p>
  
  <pre><code>awk '{
          L=length($1);
          for(i=1;i&lt;=L;i++) {
              B=substr($1,i,1);
              T[i][B]++;
          }
        } 
        END{
              for(BI=0;BI&lt;4;BI++) {
                  B=(BI==0?"A":(BI==1?"C":(BI==2?"G":"T")));
                  printf("%s",B); 
                  for(i in T) {
                      tot=0.0;
                      for(B2 in T[i]){
                          tot+=T[i][B2];
                      }
                  printf("\t%0.2f",(T[i][B]/tot));
                  } 
              printf("\n");
              }
  }' input.txt
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-08-04T07:32:18.193Z
Answers.Text: Can we do same for taking dinucleotide at time (eg; AA, TT, CG etc); the output matrix would reduce to half?
Answers.UserDisplayName: kashiff007
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-15T11:21:09.063Z
Answers.LastActivityDate: 2017-12-16T16:34:41.507Z
Answers.LastEditDate: 2017-12-16T16:34:41.507Z
Answers.UserDisplayName: Pierre
Answers.Body: <p>Here is one example of how to do this with a bit of python. Alternatively one could create strings of each column and using <code>letterFrequency()</code> from the <code>Biostrings</code> package.</p>
  
  <pre class="lang-python prettyprint-override"><code>#Make a list of hashes
  hl = []
  for i in range(9):
      hl.append({'A': 0, 'C': 0, 'G': 0, 'T': 0})
  
  f = open("foo.txt")  # CHANGE ME
  nLines = 0
  for line in f:
      for idx, c in enumerate(line.strip()):
          hl[idx][c] += 1
      nLines += 1
  f.close()
  
  nLines = float(nLines)
  for char in ['A', 'C', 'G', 'T']:
      print("{}\t{}".format(char, "\t".join(["{:0.2f}".format(x[char]/nLines) for x in hl])))
  </code></pre>
  
  <p>The output of your example is then:</p>
  
  <pre><code>A   0.13    0.40    0.13    0.00    0.00    0.40    0.47    0.07    0.40
  C   0.40    0.27    0.00    0.00    0.00    0.07    0.07    0.13    0.27
  G   0.20    0.20    0.67    1.00    0.00    0.47    0.20    0.73    0.20
  T   0.27    0.13    0.20    0.00    1.00    0.07    0.27    0.07    0.13
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-08-04T07:37:02.727Z
Answers.Text: how to change enumerate if I want to take dinucleotides at time (eg; AA, TT, CG etc); the output matrix would reduce to half?
Answers.UserDisplayName: kashiff007
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T21:01:48.197Z
Answers.LastActivityDate: 2017-12-15T10:22:28.050Z
Answers.LastEditDate: 2017-12-15T10:22:28.050Z
Answers.UserDisplayName: Devon Ryan
Answers.Body: <p>Here's a Perl approach:</p>
  
  <pre class="lang-perl prettyprint-override"><code>#!/usr/bin/env perl
  use strict;
  my %counts;
  
  ## Read the input file line by line
  while (my $line = &lt;&gt;) {
    print;
    ## remove trailing '\n' characters
    chomp $line;
    ## split the line into an array at every character
    my @columns=split(/./, $line); 
    ## iterate over the array from the first to the last position
    for my $i (0..$#columns){
      ## The nucleotide at this position
      my $nt = $columns[$i];    
      ## Save the count in this hash of arrays. The keys are the
      ## nucleotides, and the value at each position is increased
      ## when that nucleotide is found in that position.
      $counts{$nt}[$i]++;
    }
  }
  ## Iterate over the counts hash
  for my $nt (sort keys(%counts)){
    print "$nt\t";
    ## dereference the array stored in the hash
    my @countsForThisNt = @{$counts{$nt}};
    ## iterate over the counts for each position for this nt
    for (my $l=0;$l&lt;=$#countsForThisNt; $l++) {#
      ## If the value for this position isn't defined,
      ## set it to 0.
      $countsForThisNt[$l]||=0;
      ## Print all the things
      printf "%.2f\t", $countsForThisNt[$l]/$., $l;
    }
    print "\n";
  }
  </code></pre>
  
  <p>Save the script somewhere in your PATH, make it executable and run:</p>
  
  <pre><code>$ foo.pl file  
  A   0.13    0.40    0.13    0.00    0.00    0.40    0.47    0.07    0.40    
  C   0.40    0.27    0.00    0.00    0.00    0.07    0.07    0.13    0.27    
  G   0.20    0.20    0.67    1.00    0.00    0.47    0.20    0.73    0.20    
  T   0.27    0.13    0.20    0.00    1.00    0.07    0.27    0.07    0.13    
  </code></pre>
  
  <p>Alternatively, <a href="https://www.youtube.com/watch?v=Be7Og9Gc_KY&amp;feature=youtu.be&amp;t=18s" rel="nofollow noreferrer">if you're into the whole brevity thing</a>, and enjoy some golfing, here's the same thing as a one-liner:</p>
  
  <pre><code> perl -ne 'chomp;@F=split(//);for$i(0..$#F){$k{$F[$i]}[$i]++}}{for $nt(sort keys(%k)){print"$nt\t";for$i(0..$#{$k{$nt}}){$g=$k{$nt}[$i]||0; printf"%.2f\t",$g/$.}print"\n";}' file
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-15T10:01:03.897Z
Answers.LastActivityDate: 2017-12-15T10:12:18.873Z
Answers.LastEditDate: 2017-12-15T10:12:18.873Z
Answers.UserDisplayName: terdon
Answers.Body: <p>benn’s answer works but is very un-R-like: you wouldn’t normally iteratively assign to items in a matrix: R has powerful vectorised operations which make this manual work unnecessary. Here’s a more idiomatic solution:</p>
  <pre class="lang-r prettyprint-override"><code>sequences = readLines('sequences.txt')
  bases_matrix = do.call(rbind, strsplit(sequences, ''))
  
  apply(bases_matrix, 2L, function (col) {
      str = DNAString(paste(col, collapse = ''))
      letterFrequency(str, letters = 'ACGT', OR = 0L, as.prob = TRUE)
  })
  </code></pre>
  <p>This uses the same Bioconductor packages. Since this is such a simple problem, it can also be written without Bioconductor:</p>
  <pre class="lang-r prettyprint-override"><code>bases = strsplit(sequences, '')
  # Use a data.frame here so we can use factors in the next step:
  # R does not support matrices of factors. Ugh.
  bases_by_column = setNames(do.call(rbind.data.frame, bases), seq_along(bases[[1L]]))
  # Ensure that every column will be a complete set of ACGT frequencies
  bases_by_column = lapply(bases_by_column, factor, c('A', 'C', 'G', 'T'))
  sapply(lapply(bases_by_columns, table), prop.table)
  </code></pre>
  <p>Using modern R idioms from the ‘magrittr’ package, I’d write this as a pipeline; this very directly shows the sequence of transformations.</p>
  <pre class="lang-r prettyprint-override"><code>do.call(rbind.data.frame, bases) %&gt;%
      setNames(seq_along(bases[[1L]])) %&gt;%
      lapply(factor, c('A', 'C', 'G', 'T')) %&gt;%
      lapply(table) %&gt;%
      sapply(prop.table)
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-21T15:53:36.027Z
Answers.Text: Double for loop is what I learned from java. In my opinion bioinformatics is about that it works, if you are only interested good code practices you should restrict to [SO](https://stackoverflow.com/) I guess.
Answers.UserDisplayName: benn
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-21T16:04:49.997Z
Answers.Text: @b.nota Nothing wrong with “it works”. But this being a Q&A site, I think we should definitely also teach best practices. in fact, there’s a problem in bioinformatics with this “it works” attitude, which leads to a proliferation of a lot of truly atrocious code. “But it works”. Yeah, just about … as long as you don’t nudge it.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-21T16:09:19.520Z
Answers.Text: I don't completely agree, no offense. With many bioinformaticians the BIO part is much more important than the best script practices. At the company where I work, they don't care about these best script practises, but do want someone that understands e.g., what a tissue resident memory cell is. That's just my point of view.
Answers.UserDisplayName: benn
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-21T16:13:05.920Z
Answers.Text: @b.nota It’s akin to lab (or generally science) best practices, to be honest: Many people use the same argument to justify sloppy lab work or dodgy use of statistics (p-hacking etc). At the moment bad software engineering is still somewhat accepted but I’m confident that this acceptance will lessen as the field matures. Just to clarify: this doesn’t so much concern your code, which works fine. But there’s a lot of genuinely bad, unmaintainable code in circulation in bioinfomatics *tools* that costs money and time to fix because of bad practices.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-21T16:19:55.093Z
Answers.Text: Yeah I agree with that part. I am PhD in biology, and learned to code with a java course here and a course in R there. I work with R almost daily, but unfortunately I still try to avoid `apply` functions (I do try them more often, but still the old school java double for loop is more intuitive for me).
Answers.UserDisplayName: benn
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2017-12-21T15:08:09.717Z
Answers.LastActivityDate: 2020-07-22T13:58:33.167Z
Answers.LastEditDate: 2020-07-22T13:58:33.167Z
Answers.UserDisplayName: Konrad Rudolph
Answers.Body: <p>Here an example in R using <code>Biostrings</code> and <code>letterFrequency</code> (as suggested by <a href="https://bioinformatics.stackexchange.com/users/77/devon-ryan">Devon Ryan</a>).</p>
  
  <pre class="lang-r prettyprint-override"><code>library(Biostrings)
  
  data &lt;- read.table("DNA.txt", stringsAsFactors = F)
  
  new &lt;- matrix(nrow = 9, ncol = 15)
  
  for(i in 1:9){
    for(j in 1:15){
      new[i,j] &lt;- substring(data[j,], i, i)
    }
  }
  
  countTable &lt;- matrix(nrow = 9, ncol = 4)
  for(i in 1:9){
    columnSeq &lt;- DNAStringSet(paste0(new[i,], collapse = ""))
    columnCounts &lt;- letterFrequency(columnSeq, letters = "ACGT", OR = 0)
    countTable[i,] &lt;- columnCounts
  }
  colnames(countTable) &lt;- c("A", "C", "G", "T")
  
  freqTable &lt;- countTable/15
  round(t(freqTable), digit = 2)
    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
  A 0.13 0.40 0.13    0    0 0.40 0.47 0.07 0.40
  C 0.40 0.27 0.00    0    0 0.07 0.07 0.13 0.27
  G 0.20 0.20 0.67    1    0 0.47 0.20 0.73 0.20
  T 0.27 0.13 0.20    0    1 0.07 0.27 0.07 0.13
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-21T15:10:38.493Z
Answers.Text: The whole `new` matrix generation can be written in a single line as `do.call(rbind, strsplit(data, ''))`. The subsequent `for` loop can likewise be replaced.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-21T15:55:59.830Z
Answers.Text: Yes, there are many roads that lead to Rome.
Answers.UserDisplayName: benn
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-15T09:07:09.480Z
Answers.LastActivityDate: 2017-12-15T19:59:40.387Z
Answers.LastEditDate: 2017-12-15T19:59:40.387Z
Answers.UserDisplayName: benn
Answers.Body: <p>Here is how you get the distributions by column using the <code>TraMineR</code> R package.</p>
  
  <pre><code>library(TraMineR)
  
  sts.data &lt;- c(
    "CAGGTAGCC",
    "CCGGTCAGA",
    "AGGGTTTGA",
    "TTGGTGAGG",
    "CAAGTATGA",
    "ACTGTATGC",
    "CTGGTAACC",
    "TATGTACTG",
    "GCTGTGAGA",
    "CAGGTGGGC",
    "TCAGTGAGA",
    "GGGGTGAGT",
    "TGGGTATGT",
    "GAGGTGAGA",
    "CAGGTGGAG"
  )
  seq &lt;- seqdef(seqdecomp(sts.data,sep=''))
  
  seqstatd(seq)
  </code></pre>
  
  <p>and here is the outcome of the <code>seqstatd</code> function</p>
  
  <pre><code>  [State frequencies]
     [1]  [2]  [3] [4] [5]   [6]   [7]   [8]  [9]
  A 0.13 0.40 0.13   0   0 0.400 0.467 0.067 0.40
  C 0.40 0.27 0.00   0   0 0.067 0.067 0.133 0.27
  G 0.20 0.20 0.67   1   0 0.467 0.200 0.733 0.20
  T 0.27 0.13 0.20   0   1 0.067 0.267 0.067 0.13
  
     [Valid states]
    [1] [2] [3] [4] [5] [6] [7] [8] [9]
  N  15  15  15  15  15  15  15  15  15
  
     [Entropy index]
     [1]  [2]  [3] [4] [5]  [6]  [7]  [8]  [9]
  H 0.94 0.94 0.62   0   0 0.78 0.87 0.62 0.94
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-21T10:51:33.163Z
Answers.LastActivityDate: 2017-12-21T10:51:33.163Z
Answers.UserDisplayName: Gilbert
Body: <p>I have some sequences shown below</p>
  <pre><code>CAGGTAGCC
  CCGGTCAGA
  AGGGTTTGA
  TTGGTGAGG
  CAAGTATGA
  ACTGTATGC
  CTGGTAACC
  TATGTACTG
  GCTGTGAGA
  CAGGTGGGC
  TCAGTGAGA
  GGGGTGAGT
  TGGGTATGT
  GAGGTGAGA
  CAGGTGGAG
  </code></pre>
  <p>Each line has 9 nucleotides.
  Consider it to be 9 columns.I want to calculate the nucleotide frequency of each nucleotide for each of the 9 columns.
  For example 1st column will have these bases C,C,A,T,C,A etc
  Out put should be something like this</p>
  <pre><code>A   0.25    0.34    0.56    0.43    0.00    0.90    0.45    0.34    0.31
  C   0.45    0.40    0.90    0.00    0.40    0.90    0.30    0.25    0.2
  G   0.00    0.00    0.34    1.00    0.30    0.30    0.35    0.90    0.1
  T   0.24    0.56    0.00    0.00    1.00    0.34    0.45    0.35    0.36
  </code></pre>
  <hr />
  <p><strong>Note</strong>, I just made up the numbers to show you the output file format</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-12-14T16:57:48.377Z
Comments.Text: I think a package in Bioconductor can do this, I don't recall exactly one. Has you searched there? Or what have you tried to calculate this? The original file is a plain text or a fasta file ( I don't recall if it is only for fasta or for other file types)?
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-12-14T17:02:58.210Z
Comments.Text: Original file is a text file only. I am thinking of using awk(associative arrays) to calculate it, still writing the script though
Comments.UserDisplayName: user3138373
Tags: <sequence-analysis><nucleotide-models>
Title: calculating nucleotide frequency per column
ThreadText: # calculating nucleotide frequency per column
  
  (Asked by: username_0 on 2017-12-14)
  
  I have some sequences shown below
  
  ```
  CAGGTAGCC
  CCGGTCAGA
  AGGGTTTGA
  TTGGTGAGG
  CAAGTATGA
  ACTGTATGC
  CTGGTAACC
  TATGTACTG
  GCTGTGAGA
  CAGGTGGGC
  TCAGTGAGA
  GGGGTGAGT
  TGGGTATGT
  GAGGTGAGA
  CAGGTGGAG
  
  ```
  
  Each line has 9 nucleotides.
  Consider it to be 9 columns.I want to calculate the nucleotide frequency of each nucleotide for each of the 9 columns.
  For example 1st column will have these bases C,C,A,T,C,A etc
  Out put should be something like this
  
  ```
  A   0.25    0.34    0.56    0.43    0.00    0.90    0.45    0.34    0.31
  C   0.45    0.40    0.90    0.00    0.40    0.90    0.30    0.25    0.2
  G   0.00    0.00    0.34    1.00    0.30    0.30    0.35    0.90    0.1
  T   0.24    0.56    0.00    0.00    1.00    0.34    0.45    0.35    0.36
  
  ```
  
  ---
  
  **Note**, I just made up the numbers to show you the output file format
  
  **Tags:** sequence-analysis, nucleotide-models
  
  > **Comments:**
  > **username_1**: I think a package in Bioconductor can do this, I don't recall exactly one. Has you searched there? Or what have you tried to calculate this? The original file is a plain text or a fasta file ( I don't recall if it is only for fasta or for other file types)?  
  > **username_0**: Original file is a text file only. I am thinking of using awk(associative arrays) to calculate it, still writing the script though  
  
  ---
  
  ## 7 answers
  
  ### Answer by username_2 on 2017-12-16. Score: 3 | Accepted answer
  
  For anyone interested in doing this from any sort of alignment file, I've implemented a [position frequency matrix function in AlignBuddy](https://github.com/biologyguy/BuddySuite/wiki/AB-Position-Frequency-Matrix).
  
  input:
  
  ```
   15 9
  seq_1      CAGGTAGCC
  seq_2      CCGGTCAGA
  seq_3      AGGGTTTGA
  seq_4      TTGGTGAGG
  seq_5      CAAGTATGA
  seq_6      ACTGTATGC
  seq_7      CTGGTAACC
  seq_8      TATGTACTG
  seq_9      GCTGTGAGA
  seq_10     CAGGTGGGC
  seq_11     TCAGTGAGA
  seq_12     GGGGTGAGT
  seq_13     TGGGTATGT
  seq_14     GAGGTGAGA
  seq_15     CAGGTGGAG
  
  ```
  
  Command:
  
  ```
  :$ alignbuddy input.phy --pos_freq_mat
  
  ```
  
  Output:
  
  ```
  ### Alignment 1 ###
  A       0.133   0.400   0.133   0.000   0.000   0.400   0.467   0.067   0.400
  C       0.400   0.267   0.000   0.000   0.000   0.067   0.067   0.133   0.267
  G       0.200   0.200   0.667   1.000   0.000   0.467   0.200   0.733   0.200
  T       0.267   0.133   0.200   0.000   1.000   0.067   0.267   0.067   0.133
  
  ```
  
  
  ---
  
  ### Answer by username_3 on 2017-12-15. Score: 5
  
  awk
  
  ```
  awk '{L=length($1);for(i=1;i<=L;i++) {B=substr($1,i,1);T[i][B]++;}} END{for(BI=0;BI<4;BI++) {B=(BI==0?"A":(BI==1?"C":(BI==2?"G":"T")));printf("%s",B); for(i in T) {tot=0.0;for(B2 in T[i]){tot+=T[i][B2];}printf("\t%0.2f",(T[i][B]/tot));} printf("\n");}}' input.txt
  
  A   0.13    0.40    0.13    0.00    0.00    0.40    0.47    0.07    0.40
  C   0.40    0.27    0.00    0.00    0.00    0.07    0.07    0.13    0.27
  G   0.20    0.20    0.67    1.00    0.00    0.47    0.20    0.73    0.20
  T   0.27    0.13    0.20    0.00    1.00    0.07    0.27    0.07    0.13
  
  ```
  
  Or, expanded for clarity:
  
  ```
  awk '{
          L=length($1);
          for(i=1;i<=L;i++) {
              B=substr($1,i,1);
              T[i][B]++;
          }
        } 
        END{
              for(BI=0;BI<4;BI++) {
                  B=(BI==0?"A":(BI==1?"C":(BI==2?"G":"T")));
                  printf("%s",B); 
                  for(i in T) {
                      tot=0.0;
                      for(B2 in T[i]){
                          tot+=T[i][B2];
                      }
                  printf("\t%0.2f",(T[i][B]/tot));
                  } 
              printf("\n");
              }
  }' input.txt
  
  ```
  
  > **Comments:**
  > **username_4**: Can we do same for taking dinucleotide at time (eg; AA, TT, CG etc); the output matrix would reduce to half?  
  
  ---
  
  ### Answer by username_5 on 2017-12-14. Score: 4
  
  Here is one example of how to do this with a bit of python. Alternatively one could create strings of each column and using `letterFrequency()` from the `Biostrings` package.
  
  ```python
  #Make a list of hashes
  hl = []
  for i in range(9):
      hl.append({'A': 0, 'C': 0, 'G': 0, 'T': 0})
  
  f = open("foo.txt")  # CHANGE ME
  nLines = 0
  for line in f:
      for idx, c in enumerate(line.strip()):
          hl[idx][c] += 1
      nLines += 1
  f.close()
  
  nLines = float(nLines)
  for char in ['A', 'C', 'G', 'T']:
      print("{}\t{}".format(char, "\t".join(["{:0.2f}".format(x[char]/nLines) for x in hl])))
  
  ```
  
  The output of your example is then:
  
  ```
  A   0.13    0.40    0.13    0.00    0.00    0.40    0.47    0.07    0.40
  C   0.40    0.27    0.00    0.00    0.00    0.07    0.07    0.13    0.27
  G   0.20    0.20    0.67    1.00    0.00    0.47    0.20    0.73    0.20
  T   0.27    0.13    0.20    0.00    1.00    0.07    0.27    0.07    0.13
  
  ```
  
  > **Comments:**
  > **username_4**: how to change enumerate if I want to take dinucleotides at time (eg; AA, TT, CG etc); the output matrix would reduce to half?  
  
  ---
  
  ### Answer by username_6 on 2017-12-15. Score: 3
  
  Here's a Perl approach:
  
  ```perl
  #!/usr/bin/env perl
  use strict;
  my %counts;
  
  ## Read the input file line by line
  while (my $line = <>) {
    print;
    ## remove trailing '\n' characters
    chomp $line;
    ## split the line into an array at every character
    my @columns=split(/./, $line); 
    ## iterate over the array from the first to the last position
    for my $i (0..$#columns){
      ## The nucleotide at this position
      my $nt = $columns[$i];    
      ## Save the count in this hash of arrays. The keys are the
      ## nucleotides, and the value at each position is increased
      ## when that nucleotide is found in that position.
      $counts{$nt}[$i]++;
    }
  }
  ## Iterate over the counts hash
  for my $nt (sort keys(%counts)){
    print "$nt\t";
    ## dereference the array stored in the hash
    my @countsForThisNt = @{$counts{$nt}};
    ## iterate over the counts for each position for this nt
    for (my $l=0;$l<=$#countsForThisNt; $l++) {#
      ## If the value for this position isn't defined,
      ## set it to 0.
      $countsForThisNt[$l]||=0;
      ## Print all the things
      printf "%.2f\t", $countsForThisNt[$l]/$., $l;
    }
    print "\n";
  }
  
  ```
  
  Save the script somewhere in your PATH, make it executable and run:
  
  ```
  $ foo.pl file  
  A   0.13    0.40    0.13    0.00    0.00    0.40    0.47    0.07    0.40    
  C   0.40    0.27    0.00    0.00    0.00    0.07    0.07    0.13    0.27    
  G   0.20    0.20    0.67    1.00    0.00    0.47    0.20    0.73    0.20    
  T   0.27    0.13    0.20    0.00    1.00    0.07    0.27    0.07    0.13    
  
  ```
  
  Alternatively, [if you're into the whole brevity thing](https://www.youtube.com/watch?v=Be7Og9Gc_KY&feature=youtu.be&t=18s), and enjoy some golfing, here's the same thing as a one-liner:
  
  ```
   perl -ne 'chomp;@F=split(//);for$i(0..$#F){$k{$F[$i]}[$i]++}}{for $nt(sort keys(%k)){print"$nt\t";for$i(0..$#{$k{$nt}}){$g=$k{$nt}[$i]||0; printf"%.2f\t",$g/$.}print"\n";}' file
  
  ```
  
  
  ---
  
  ### Answer by username_7 on 2017-12-21. Score: 2
  
  username_8’s answer works but is very un-R-like: you wouldn’t normally iteratively assign to items in a matrix: R has powerful vectorised operations which make this manual work unnecessary. Here’s a more idiomatic solution:
  
  ```r
  sequences = readLines('sequences.txt')
  bases_matrix = do.call(rbind, strsplit(sequences, ''))
  
  apply(bases_matrix, 2L, function (col) {
      str = DNAString(paste(col, collapse = ''))
      letterFrequency(str, letters = 'ACGT', OR = 0L, as.prob = TRUE)
  })
  
  ```
  
  This uses the same Bioconductor packages. Since this is such a simple problem, it can also be written without Bioconductor:
  
  ```r
  bases = strsplit(sequences, '')
  # Use a data.frame here so we can use factors in the next step:
  # R does not support matrices of factors. Ugh.
  bases_by_column = setNames(do.call(rbind.data.frame, bases), seq_along(bases[[1L]]))
  # Ensure that every column will be a complete set of ACGT frequencies
  bases_by_column = lapply(bases_by_column, factor, c('A', 'C', 'G', 'T'))
  sapply(lapply(bases_by_columns, table), prop.table)
  
  ```
  
  Using modern R idioms from the ‘magrittr’ package, I’d write this as a pipeline; this very directly shows the sequence of transformations.
  
  ```r
  do.call(rbind.data.frame, bases) %>%
      setNames(seq_along(bases[[1L]])) %>%
      lapply(factor, c('A', 'C', 'G', 'T')) %>%
      lapply(table) %>%
      sapply(prop.table)
  
  ```
  
  > **Comments:**
  > **username_8**: Double for loop is what I learned from java. In my opinion bioinformatics is about that it works, if you are only interested good code practices you should restrict to [SO](https://stackoverflow.com/) I guess.  
  > **username_7**: @b.nota Nothing wrong with “it works”. But this being a Q&A site, I think we should definitely also teach best practices. in fact, there’s a problem in bioinformatics with this “it works” attitude, which leads to a proliferation of a lot of truly atrocious code. “But it works”. Yeah, just about … as long as you don’t nudge it.  
  > **username_8**: I don't completely agree, no offense. With many bioinformaticians the BIO part is much more important than the best script practices. At the company where I work, they don't care about these best script practises, but do want someone that understands e.g., what a tissue resident memory cell is. That's just my point of view.  
  > **username_7**: @b.nota It’s akin to lab (or generally science) best practices, to be honest: Many people use the same argument to justify sloppy lab work or dodgy use of statistics (p-hacking etc). At the moment bad software engineering is still somewhat accepted but I’m confident that this acceptance will lessen as the field matures. Just to clarify: this doesn’t so much concern your code, which works fine. But there’s a lot of genuinely bad, unmaintainable code in circulation in bioinfomatics *tools* that costs money and time to fix because of bad practices.  
  > **username_8**: Yeah I agree with that part. I am PhD in biology, and learned to code with a java course here and a course in R there. I work with R almost daily, but unfortunately I still try to avoid `apply` functions (I do try them more often, but still the old school java double for loop is more intuitive for me).  
  
  ---
  
  ### Answer by username_8 on 2017-12-15. Score: 2
  
  Here an example in R using `Biostrings` and `letterFrequency` (as suggested by [username_5 username_5](https://bioinformatics.stackexchange.com/users/77/username_5-username_5)).
  
  ```r
  library(Biostrings)
  
  data <- read.table("DNA.txt", stringsAsFactors = F)
  
  new <- matrix(nrow = 9, ncol = 15)
  
  for(i in 1:9){
    for(j in 1:15){
      new[i,j] <- substring(data[j,], i, i)
    }
  }
  
  countTable <- matrix(nrow = 9, ncol = 4)
  for(i in 1:9){
    columnSeq <- DNAStringSet(paste0(new[i,], collapse = ""))
    columnCounts <- letterFrequency(columnSeq, letters = "ACGT", OR = 0)
    countTable[i,] <- columnCounts
  }
  colnames(countTable) <- c("A", "C", "G", "T")
  
  freqTable <- countTable/15
  round(t(freqTable), digit = 2)
    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
  A 0.13 0.40 0.13    0    0 0.40 0.47 0.07 0.40
  C 0.40 0.27 0.00    0    0 0.07 0.07 0.13 0.27
  G 0.20 0.20 0.67    1    0 0.47 0.20 0.73 0.20
  T 0.27 0.13 0.20    0    1 0.07 0.27 0.07 0.13
  
  ```
  
  > **Comments:**
  > **username_7**: The whole `new` matrix generation can be written in a single line as `do.call(rbind, strsplit(data, ''))`. The subsequent `for` loop can likewise be replaced.  
  > **username_8**: Yes, there are many roads that lead to Rome.  
  
  ---
  
  ### Answer by username_9 on 2017-12-21. Score: 1
  
  Here is how you get the distributions by column using the `TraMineR` R package.
  
  ```
  library(TraMineR)
  
  sts.data <- c(
    "CAGGTAGCC",
    "CCGGTCAGA",
    "AGGGTTTGA",
    "TTGGTGAGG",
    "CAAGTATGA",
    "ACTGTATGC",
    "CTGGTAACC",
    "TATGTACTG",
    "GCTGTGAGA",
    "CAGGTGGGC",
    "TCAGTGAGA",
    "GGGGTGAGT",
    "TGGGTATGT",
    "GAGGTGAGA",
    "CAGGTGGAG"
  )
  seq <- seqdef(seqdecomp(sts.data,sep=''))
  
  seqstatd(seq)
  
  ```
  
  and here is the outcome of the `seqstatd` function
  
  ```
    [State frequencies]
     [1]  [2]  [3] [4] [5]   [6]   [7]   [8]  [9]
  A 0.13 0.40 0.13   0   0 0.400 0.467 0.067 0.40
  C 0.40 0.27 0.00   0   0 0.067 0.067 0.133 0.27
  G 0.20 0.20 0.67   1   0 0.467 0.200 0.733 0.20
  T 0.27 0.13 0.20   0   1 0.067 0.267 0.067 0.13
  
     [Valid states]
    [1] [2] [3] [4] [5] [6] [7] [8] [9]
  N  15  15  15  15  15  15  15  15  15
  
     [Entropy index]
     [1]  [2]  [3] [4] [5]  [6]  [7]  [8]  [9]
  H 0.94 0.94 0.62   0   0 0.78 0.87 0.62 0.94
  
  ```
  
  
  ---
  
  