procedure for extracting raw read counts

[NitaC]

Good morning everyone.

I am a newb and this is my first time touching any sequencing data. It just so happens that the best way to learn is to dive right in. Incidentally I have hit a bit of a stumbling block and was hoping someone here could shed some light on it.

I am trying to extract the raw read counts for transcripts mapped only to the forward strand. The mapping was done using Gsnap (I did not do the mapping.) We did paired-end sequencing using Illumina's HiSeq platform. Again, this was my first time working with sam/bam files and it took some time to understand the formats. Here is what I did:

To get only the forward read counts I took the sam file I was supplied and converted it back to a bam file. Then using the 0x0040 flag I used a command to get just the forward reads and redirected that output to another sam file. Then, looking at column 3 in the sam file, I saw either the name of the transcript or a *. It is my understanding that the "*" is for those reads that did not align. So I wrote a little bash script to filter those out. I was told by a co-worker that using column 3 I could get the raw counts of the mapped transcripts. So, I wrote another little script to do that. I know my script is working correctly, however, my numbers do not match the numbers my co-worker has. In fact, I have counts for some transcripts that my co-worker reported as not detected at all. I know there are tools to automatically get counts, but I want to actually understand these file formats before using some black box of tools.

This is how one line looks in the sam file:
[anitabow@kirby2 spikeIn]$ head -n 35 DR38spikeInForward.sam|tail -n 1
FCC1WG7ACXX:5:1101:3115:1979#GGACCCAT 99 ERCC-00130 851 40 100M = 958 207 CTTGAAGCAGATGATATCCCGATGCTGACGGCTTACAATAAACGTGATCAAAAACTGCCTGATTTTATACCGACCGCCGGAAGGGATCACATTATGGTCA CCCFFFFFHHHHHJJJJJJJJJJJJJJJJJJJJIFFHHDGIJIJDEHJIIJJJJIHHIJJHHHHHHHFFFFFDDDDDDDDDDDDDDDDDDDDDDEDDCCD MD:Z:100 NH:i:1 HI:i:1 NM:i:0 SM:i:40 XQ:i:40 X2:i:0

As you can see, the third column has the name of one transcript. This is what I have been using to extract my counts. So basically, my script will count every time it sees ERCC-00130 for example.

Thanks in advance for any help.

[maubp]

You need to read more about the FLAG field in SAM/BAM. Here bits 0x4 and 0x10 will tell you if a read is mapped and to which strand.

You can't just look at column 3 (RNAME), as that can be a reference name if the read's mate was mapped even if the read itself was not.

[NitaC]

Hi maubp!

It is my understanding that 0x4 tells you that the read was unmapped and that 0x10 indicates that the read is on the reverse strand. I don't think either of those flags really address my issue. I've already filtered for unaligned reads.

And I'm not sure I understand your second response. If I've filtered out unaligned reads and filtered for only the forward strand, ultimately I still need some sort of reference name or gene id to count right? Are you saying that after filtering for forward reads, there may still be reads in my results that have not actually mapped to the forward strand but are listed because their mate did?

As it turns out, my little pipeline actually matches the results using htseq-count so I'm a little more confident in what I did. And my co-worker's counts were from one sample whereas mine were from all the samples combined. That little tidbit would've saved me from a lot of headaches.

[maubp]

Are you saying that after filtering for forward reads, there may still be reads in my results that have not actually mapped to the forward strand but are listed because their mate did?

Yes in general. This is deliberate so that when coordinate sorted the unmapped read is next to its mapped partner in the BAM file. See "Recommended Practise for the SAM format" in the specification. Also "Bit 0x4 is the only reliable place to tell whether the segment is unmapped".

[NitaC]

Ohhh, good to know. Thank you.

[Simon Anders]

Also, remember that a read that aligns to multiple places will be listed in multiple lines in the SAM file.