count mapped reads

[Laurent Manchon]

--Hi,

i need to count the mapped reads from RSEM output file '*.genome.bam' to compare with
the accepted.results.bam output generated by tophat.
So, i have tried to use this command: samtools view -c -F 4 my.bam
For tophat output it runs fine but for RSEM i obtain a big value, i don't know if this command is correct for RSEM outputs.

thank you --

[b...@cs.wisc.edu]

Hi Laurent,

RSEM output every alignment a read has and also their associated posterior
probabilities. Thus what you have with your command will be the total
number of alignments instead of reads. If you want to know the mapped
reads, you can go to sample_name.stat/sample_name.cnt. The first line
contains N0, N1, N2, N, which are bowtie statistics. They represent number
of reads failed to align, number of reads alignable, number of reads that
aligning to more than 200 places and therefore filtered and the total
number of reads.

Best,
Bo

> --
> RSEM website: http://deweylab.biostat.wisc.edu/rsem/
> ---
> You received this message because you are subscribed to the Google Groups
> "RSEM Users" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to rsem-users+...@googlegroups.com.
> To post to this group, send email to rsem-...@googlegroups.com.
> Visit this group at http://groups.google.com/group/rsem-users?hl=en-US.
>
>
>

[Shawn Driscoll]

Bo-

Is it also reliable to gain this count by adding in the --sampling-for-bam option and then using samtools view to filter out alignments with MAPQ == 0?  It seems like that option makes it so the alignments selected based on probability have MAPQ == 100 while everything else is 0. 

[Bo Li]

Hi Shawn,

Not exactly, since there is a chance to sample the read from noise although this read has alignments. The most reliable way is to write a script by yourself, which count number of reads do have an alignment in the bam/sam file.

Best,
Bo

[Colin Dewey]

Probably a simpler way to get at the count of aligned reads (or fragments) is to sum the "expected count" column of the gene.results or isoform.results files.  This gives you the expectation of the number of aligned reads (fragments).  Recall that RSEM treats the true alignment of a read as being unknown, so it can only provide an estimate of the number of reads that truly align.

Colin

[c...@ncgr.org]

Hey All,

I was curious if you might comment on why so few of my aligned reads are counted in the sum of my expected counts file.  The total read count is ~44 M.  In my expected counts file I sum to ~22 M reads.  I wrote a quick one liner to count the number of reads where the bitflag (field 2) was not set for unaligned using a bitwise and with 0x04 (can be done for any bam file) and then sorting and unique counting the first field.  This produced a very appropriate ~40 M reads aligned, which is approximately what I receive using bwa as well which does not use multi-record bam files and can be counted quickly using flagstat.  Using an independent (and much simpler) counting procedure based on unique read counting (parsed from bwa alignments) and iterative apportioning, based on unique ratios, I obtain a count that is ~40 M reads.

The context of this is a denovo assembly using an in-house platform based on ABySS and Cap3.  Much care was taken to avoid unnecessary redundancy in the assembly process, however, it still exists due to large homologous regions flanked by differential sequence (definitely not confident enough to assume these are isoforms, and no reference genome is available).

Thoughts?

Thank you, I really like the method but lack the understanding to completely rationalize this using the noise transcript,

Connor Cameron 

[Colin Dewey]

Hi Connor,

Are your data paired-end by any chance?  If so, then it's important to note that the expected counts are for *fragments*, not reads, so one pair of reads counts as one fragment.

Colin

[c...@ncgr.org]

Excellent, I had thought about this, but didn't really revisit it.

Thank you for your time,

Connor Cameron

[Laurent Manchon]

okay good,
thank you for this help,
but be careful when i sum the "expected count" column of the gene.results
i don't get exactly the same number recorded by bowtie in sample.cnt file.

[Colin Dewey]

Hi Laurent,

The numbers in the sample.cnt file are summaries of the raw alignment generated by bowtie using liberal parameters (assuming you let RSEM do the alignment for you).  The number of reads found to be aligned in the raw alignment will generally be higher than the expected count because RSEM will assign very small probabilities to alignments with many mismatches (at high quality score positions).  The expected count should generally be a more accurate estimate of how many reads are truly aligned.

Best,

Colin

[Shawn Driscoll]

Interesting information. It's such a bizzarr concept that aligned reads may not be truly aligned but from what I can tell RSEM does work quite well. I was using the genome bam file generated by RSEM to compare back to a control BAM to sort of benchmark the RSEM selected alignments and to do that I need the bam to contain only the selected alignment and no duplicates. It would be great if RSEM could produce that output...maybe I'm bit seeing the option? Either that or if the genome bam is annotated in a way that only those alignments could be extracted with a quick samtools view call. I was thinking that the mapq field could be used for this but I guess not?

[b...@cs.wisc.edu]

Hi Shawn,

To achieve this purpose, you can use '--sampling-for-bam' option, which
randomly samples an alignment based on the posterior weights. If you use
this option, all alignment's MAPQ field is 0 execpt those selected.

If you want exactly the best weighted alignments are selected, you need to
write your own script.

Hope it helps,
Bo

[Shawn Driscoll]

I see - and if I want to go with my own script then I need to leave off the '--sampling-for-bam' option, correct? otherwise the weights are all set to 1/0, right?

[archana bhardwaj]

Hello everyone 

I am new to NGS . I have paired end sequence data for differential gene expression analysis. First i did the prepossing of my data such as adapter removal, ambiguous character and quality filtration . Now i want to assemble my data in form of contigs. I run trininty.pl to get the assembled contigs. There are so many output files . What is the exact contig files??? Inchworm contigs or bundled.fasta ??? How can i get the read count of each contig ??? 

waiting for reply 

[Mafalda Ferreira]

Hello,

I'm sorry if this is out of the blue since this was answered more than a year ago...
I was struggling with the exact same thing. I want to know the number of reads that mapped so I can have an idea of the mapping % and success. 
I used samtools flagstats first but I got a huge number! However, now I understand that it's due to the multiple alignments that are allowed.

Next, I checked the library sizes in edgeR and they seem about right, giving me aproximately 460M reads in total for all samples (I doubled the expected counts, since they are counts for fragments, right?), when the input is 650M, which means that is some amount of failure, which is ok. Now I was looking to the my_sample.rsem.out.cnt file in the my_sample.rsem.out.stat/ folder and I can see that the #total reads is exactly half of what it should be (counted by the input fastqcs). Am I reading this wrong or this is a count for fragments and I should double it too? Here is an example for 1 sample

                Library size in EdgeR: 16061092 (the double would be the number of reads right?)

     # of alignable reads in .cnt file: 16279180

         # of total "reads" in .cnt file: 22002705 (which is exactly half of what it should be by the .cnt)

# of input reads counted in fastqc: 44005410 (which is exactly double of what it should be by the .cnt)

I'm pretty sure that I didn't count the input wrong because I've did it with a custom script, FastQC and with wc -l and they all converge in the same number.

This is the head of my_sample.rsem.out.cnt:

5723525 16279180 0 22002705

16222463 56717 6069312

37528223 3

0       5723525

1       10209868

2       2878169

3       645611

4       970148

5       263091

6       386832

[Mafalda Ferreira]

I'm sorry, I wrote something wrong there. 

In:

"# of total "reads" in .cnt file: 22002705 (which is exactly half of what it should be by the .cnt)"

I meant to say:   

# of total "reads" in .cnt file: 22002705 (which is exactly half of what it should be by the fastqcs)

Thanks

[Bo Li]

Hi Mafalda,

If you have paired-end reads, the two mates of a pair are counted as one
read.

Hope it helps,
Bo

>>> RSEM website: http://deweylab.biostat.wisc.edu/rsem/ [1]

>>> ---
>>> You received this message because you are subscribed to the Google
>> Groups
>>> "RSEM Users" group.
>>> To unsubscribe from this group and stop receiving emails from it,
>> send an
>>> email to rsem-users+...@googlegroups.com.
>>> To post to this group, send email to rsem-...@googlegroups.com.
>>> Visit this group at

>> http://groups.google.com/group/rsem-users?hl=en-US [2].
>>>
>>>
>>>
>
> --
> RSEM website: http://deweylab.biostat.wisc.edu/rsem/ [3]

> ---
> You received this message because you are subscribed to the Google
> Groups "RSEM Users" group.
> To unsubscribe from this group and stop receiving emails from it,
> send an email to rsem-users+...@googlegroups.com.
> To post to this group, send email to rsem-...@googlegroups.com.

> Visit this group at http://groups.google.com/group/rsem-users [4].
>
>
> Links:
> ------
> [1]
> http://www.google.com/url?q75http%3A%2F%2Fdeweylab.biostat.wisc.edu%2Frsem%2F46sa75D46sntz75146usg75AFQjCNEU-wRL_aNCO7ziCWa-Y12BbbySXw
> [2] http://groups.google.com/group/rsem-users?hl75en-US
> [3] http://deweylab.biostat.wisc.edu/rsem/
> [4] http://groups.google.com/group/rsem-users