How to calculate proportion of reads with each base at every reference position

[lindseyjane]

I have a need to calculate the following, and wonder if anyone else has done something similar or could advise me on how to do it. The question is:

For each base position on a reference sequence, what is the proportion of mapped reads with an A/G/C/T at this position?
e.g. for reference position x, 50% of reads show an A/50% of reads show a T.

At the moment I have Illumina paired end 51bp data and some 76bp data and have mapped the reads to a reference using the maq software. Any advice oon how to do this type of downstream analysis would be much appreciated, thanks.

Lindsey

[spenthil]

Are you familiar with python or perl? This can be be done pretty easily with BioPython/BioPerl.

Hopefully someone else can chime in with an existing solution.

[lindseyjane]

Yes, I am learning perl. Does anyone have an existing solution?

[spenthil]

AFAIK, nothing existing but you could use bioperl to accomplish it without too much work (my guess is a few hours tops).

[simonandrews]

I had a quick go since this seemed a fairly simple task. The documentation for mapview output isn't very clear so I may not be handling the sequence the right way when it's a reverse hit, but you can at least use this to get an idea for how this could be done.

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

my @mapview_files = @ARGV;

foreach my $file (@mapview_files) {
  process_file ($file);
}

sub process_file {
  my ($file) = @_;
  open (IN,$file) or die "$file: $!";

  # Data Structure is a hash of chromosomes
  # containing arrays of positions where each
  # position is a hash of GATC mapping to a count
  my %chrs;

  while (<IN>) {
    chomp;
    # Ignore headers
    next if (/^#/);

    # Extract data
    my ($chr,$pos,$seq) = (split(/\t/))[1,2,14];

    # Split sequence into bases
    my @seq = split(//,$seq);

    # Add each base to the data structure
    for my $offset (0..$#seq) {
      ++$chrs{$chr}->[$pos+$offset]->{uc($seq[$offset])};
    }
  }

  # Print a header
  print join("\t",qw(File Chr Pos G A T C)),"\n";

  # Go through each chromosome
  foreach my $chr (sort keys %chrs) {

    # Go through the positions on that chromosome
    my @positions = @{$chrs{$chr}};
    for my $position (1..$#positions) {
      my @line = ($file,$chr,$position);

      # Get the counts for each base
      foreach my $base qw(G A T C) {
	if (exists $positions[$position]->{$base}) {
	  push @line, $positions[$position]->{$base};
	}
	else {
	  push @line, 0;
	}
      }

      # Print the result
      print join("\t",@line),"\n";
    }
  }
}

[lindseyjane]

Thank you so very much for taking the time to produce this code. It is an excellent starting point for me and I appreciate the time you have spent.

I hope soon I will be able to produce perl code this quickly myself!

[krobison]

Hmm, perhaps the site should have a wiki for contributed code (I once posted a Perl program as well)

[lindseyjane]

Hmm, perhaps the site should have a wiki for contributed code (I once posted a Perl program as well)

Excellent idea, yes please! That would be so useful.