Parsing sam file

Question

[Sidney]

I have no idea on how to start the program.

Question is:

Write a Perl program that parses out the attached sequence read alignment file (6.SAM) to count how many reads a gene produced.

Output should be :

Gene ID:            Number of reads aligning:
gene29              3
gene27              6

[Sukhi Singh]

Homework???

[russhh]

all his other posts seem to have been

[Sidney]

Practice exercises, but I am trying to simply understand how to go about solving since I do not understand the question. So far I have :

foreach my $id (@{$ids}) {
    my $alignments = $sam -> get_features_by_location(-seq_id => $transcript_id, -iterator => 1);

    while (my $alignment = $alignments -> next_seq) {
        my $read_name = $alignment -> query -> name;

        if (exists($global_read_occurrences{$read_name})) {
            $global_read_occurrences{$read_name}++;
        }
        else {
            $global_read_occurrences{$read_name} = 1;
        }
    }
}

[Devon Ryan]

What do you interpret the question as meaning? We can start by simply telling you if your interpretation of the program you need to write is correct or not. Also, can you post a line or two of the SAM file? I'm guessing that it's aligned to the transcriptome, but if not that'll change a bit how the program would work.

[swbarnes2]

You don't need that 'if' statement.

$hash{whatever}++;

works fine whether or not 'whatever' is a key in the hash already.

For people not doing perl homework, samtools idxstats will count this up for you, you could also do it with unix tools: cut, sort, uniq -c, sort

[Zev.Kronenberg]

I love perl, but anyone assigning students perl homework is mad.