Please help me with bowtie2 output and samtools report

Question

[seta]

Hi friends,

I have made a de novo transcriptome assembly, followed by mapping back read to the assembly using bowtie2 with this command:

bowtie2 -x assembly_4 -1 file_1.fq -2 file_2.fq --local --no-unal -p 8 -S assembly_4.sam

The output of bowtie2 is here:

182786653 reads; of these:
  182786653 (100.00%) were paired; of these:
    6674426 (3.65%) aligned concordantly 0 times
    125631460 (68.73%) aligned concordantly exactly 1 time
    50480767 (27.62%) aligned concordantly >1 times
    ----
    6674426 pairs aligned concordantly 0 times; of these:
      978106 (14.65%) aligned discordantly 1 time
    ----
    5696320 pairs aligned 0 times concordantly or discordantly; of these:
      11392640 mates make up the pairs; of these:
        8790860 (77.16%) aligned 0 times
        1476379 (12.96%) aligned exactly 1 time
        1125401 (9.88%) aligned >1 times
97.60% overall alignment rate

When I used the flagstat with samtools (samtools flagstat assembly_4.bam > report.stats), the report was:

356782446 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
356782446 + 0 mapped (100.00%:-nan%)
356782446 + 0 paired in sequencing
178733950 + 0 read1
178048496 + 0 read2
352224454 + 0 properly paired (98.72%:-nan%)
355623980 + 0 with itself and mate mapped
1158466 + 0 singletons (0.32%:-nan%)
2758204 + 0 with mate mapped to a different chr
2576765 + 0 with mate mapped to a different chr (mapQ>=5)

I want to know how the properly paired is 98.72% while based on bowtie2 the percentage of aligned concordantly is 96.35%. Why there is such a difference?

Answer 1

[Devon Ryan]

You appear to have filtered the BAM or SAM file at some point (or are looking at the wrong file). You started with ~170 million aligned reads but there are only a ~36 million in the BAM file, which obviously can't happen.

[seta]

Thanks, Devon. I used the command

./samtools view -S -bu assembly_4.sam | samtools sort -@ 3 - assembly_4.sorted

for converting the BAM to SAM. Also, I checked the number of alignment in bam file that returned "356,782,446". So, I thought everythings is OK. Sorry, I didn't undestand your mean from ~36 milion read in bam file!. Please correct me if I'm wrong as this is my first experience.

[Devon Ryan]

Ah, I misread ~356 million as 36 million, mea culpa!

Anyway, the difference in percentages is due to you telling bowtie2 to not print unmapped reads to the SAM file. You'll note that bowtie2 says there are 125631460+50480767 pairs of reads that should be properly paired. That equals 352224454 total reads, which is exactly what samtools reports too. The only difference is the total number of reads in the file, due to you having used --no-unal with bowtie2.

[seta]

Thanks for clearing me.