Question: How to correctly generate alignment Stat in RNA-Seq
0
gravatar for rohitsatyam102
4 months ago by
rohitsatyam102120 wrote:

Hi All!!

While aligning RNA-Seq data and quantifying the actual number of reads mapped to my transcriptome, I realised that samtools flagstat does not report the actual alignied reads stats. I used RSEM for my alignment and expression calculation. I am aware the samtools counts secondary alignments thereby increasing total number of reads that pass QC and therefore influence the alignment percentage that it throws as the output. Source here

Is there any other tool that can precisely give the actual number of reads correctly mapping in pairs to the transcriptome? Somewhat like the log.out file in STAR does.

sequencing rna-seq alignment • 297 views
ADD COMMENTlink modified 4 months ago by brianj.park50 • written 4 months ago by rohitsatyam102120
0
gravatar for Devon Ryan
4 months ago by
Devon Ryan96k
Freiburg, Germany
Devon Ryan96k wrote:

Just exclude secondary alignments then: samtools view -hF 256 foo.bam | samtools flagstat -

ADD COMMENTlink written 4 months ago by Devon Ryan96k

Hi Devon. Used the above command and I am getting the following stats

102767366 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
97087946 + 0 mapped (94.47% : N/A)
102767366 + 0 paired in sequencing
51383683 + 0 read1
51383683 + 0 read2
97087946 + 0 properly paired (94.47% : N/A)
97087946 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

However, from the FastQC report that I have, I know my sample has 88.5 million reads (each in R1 and R2). This 2 x 88.5 million should be equal to the number of reads I get in QC Passed reads ~177 million. How should I account for the rest of the ~74 million reads?

ADD REPLYlink modified 4 months ago • written 4 months ago by rohitsatyam102120

They're not in the BAM file, so likely you excluded them elsewhere.

ADD REPLYlink written 4 months ago by Devon Ryan96k
1

Thanks @Devon

Found a way around. I revisited the RSEM source code and read it. I found that the Log.final.out file is generated in the sample_temp folder which is deleted at the end. Using --keep-intermediate-files we can preserve the alignment stats and other logs too.

ADD REPLYlink written 4 months ago by rohitsatyam102120
0
gravatar for brianj.park
4 months ago by
brianj.park50
Montréal, Canada
brianj.park50 wrote:

Check RSeQC: in particular bam_stat.py

ADD COMMENTlink written 4 months ago by brianj.park50
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1482 users visited in the last hour