Question: flagstats results show 0.0 % properly paired reads
0
gravatar for SOHAIL
4.2 years ago by
SOHAIL290
Beijing Institute of Genomics, CAS.
SOHAIL290 wrote:

Hi Everyone!

I used samtools Flagstats walker. I am getting 0.0% properly paired reads in the BAM file. I sorted my BAM file and marked the duplicates as well. is there any BAM FLAG is missing in my file or mapping results are bad quality. Note that my average mapping quality (MAPQ) is more than 30.. i used bwa-mem for mapping to reference genome. i did not mark any duplicates in this step.

can you please help help me tracking this issue?

My command arguments were:

samtools flagstat /CHG000691/CHG000691_PE_02AL3.sam > /Flagstat/CHG000691/lanewise/CHG000691_PE_02AL3_flagstat.txt

Output:
75031887 + 0 in total (QC-passed reads + QC-failed reads)
733759 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
74460957 + 0 mapped (99.24%:-nan%)
74298128 + 0 paired in sequencing
37149064 + 0 read1
37149064 + 0 read2
0 + 0 properly paired (0.00%:-nan%)
73397388 + 0 with itself and mate mapped
329810 + 0 singletons (0.44%:-nan%)
2837608 + 0 with mate mapped to a different chr
982994 + 0 with mate mapped to a different chr (mapQ>=5)

Thank you very much!

alignment next-gen • 2.4k views
ADD COMMENTlink modified 19 months ago • written 4.2 years ago by SOHAIL290
1

It would be helpful if you could answer these questions -

1) What was your bwa mem command?

2) What was the library type? (e.g. fragment, LMP [plus CLRS, Nextera, etc.])

3) What was the expected insert size and orientation?

4) Are you sure the sequencing was paired?

5) If you used 2 files (rather than interleaved), are you sure the two files go together?

6) What kind of preprocessing did you do? (For example, fastx toolkit is notorious for destroying pairing order and should never be used with paired reads, or ever, IMO).

7) What kind of experiment is it?  DNA, RNA, etc.

8) What version of bwa are you using?

You can validate that the names appear to be correctly paired with the BBMap package:

reformat.sh in1=r1.fq in2=r2.fq verifypaired

If the reads are interleaved in a single file, use the "verifyinterleaved" flag instead.

In any case, it's odd that 0 reads are paired.  I would expect there to be a few just by random chance out of 75 million (though that depends on the genome size).  I'm guessing that bwa treated these as single-ended.

ADD REPLYlink written 4.2 years ago by Brian Bushnell16k

1) What was your bwa mem command?

    bwa mem -t 5 -R '@RG\tID:02AL1\tSM:HMN15-1\tLB:hs\tPL:illumina' -P -M ./human_g1k_v37.fasta ./S17_02A_CHG000691_L001_R1.fastq.gz ./S17_02A_CHG000691_L001_R2.fastq.gz > ./CHG000691/CHG000691_PE_02AL1.sam
I think it might be due to -P option, but i am not sure, because it is written in BWA manual:
        "In the paired-end mode, perform SW to rescue missing hits only but do not try to find hits that fit a proper pair. " 

2) What was the library type? (e.g. fragment, LMP [plus CLRS, Nextera, etc.])

    Fragment-Sonicator
3) What was the expected insert size and orientation?
    Forward-Reverse (FR)   (-----> <----)  
    insert size 300-500

4) Are you sure the sequencing was paired?
    
    yes, I am sure the sequencing was paired-end. as you can see the flagstat output:
        
        74298128 + 0 paired in sequencing
        37149064 + 0 read1
        37149064 + 0 read2
        

5) If you used 2 files (rather than interleaved), are you sure the two files go together?
    yes.. :) 

6) What kind of preprocessing did you do? (For example, fastx toolkit is notorious for destroying pairing order and should never be used with paired reads, or ever, IMO).
    No preprocessing was performned.    

7) What kind of experiment is it?  DNA, RNA, etc.
        Whole Genome DNA Sequencing

8) What version of bwa are you using?

    BWA-Version: 0.7.12-r1039

ADD REPLYlink written 4.2 years ago by SOHAIL290

You could easily check if -P is the culprit, by repeating the alignment omitting it.

ADD REPLYlink written 4.2 years ago by h.mon27k

Can you please post the output of head for the two FASTQ files? If they're gzipped, you can do this by typing zcat [fastq_file] | head

ADD REPLYlink modified 4.2 years ago • written 4.2 years ago by Dan D6.8k

@Dan D 

here is the output:

First File:

[sohail@datatrans CHG000691]$ cat S17_02A_CHG000691_L001_R1.fastq | head

@ST-E00169:61:H2CCHCCXX:1:1101:7475:1713 1:N:0:TCCGCGAA
CCACAAAAGGAGTATGTCAAAACTGGTCTATCAAAAGAAAGGTTCATCTCTTGGAGGTGAATGCACATATCACAAAAAAGTTTCCAAGATGGCTCCTGTCTAGATTATATGTGAAGATATATCATTTTCTCCCAAAGGCTGCATAGCGCT
+
AAFFFAFKKKKKAAF,AFK7,K,<FKF<FA7FF,AKFKKKKK7FFKK<,7F<FAKFFAFKKKK,<,<<<A<KA,A,7,K<,,,,,,,A<A,,77,<7A,,,,F,,7,A,AF7FA7A77F<,,<,,,AFAKK,A<F7,<7A,(,7,,<,7A
@ST-E00169:61:H2CCHCCXX:1:1101:8815:1713 1:N:0:TCCGCGAA
TGATAGAGCAGTATTGAAACACTCTTTTTGTGGAAAATGCAAGTGGATATTTGGATATCTTGGAGGATTTAGTTGGAAGTGGGAATTCACATAAAAGGTATACAGCAGCATTCTCAGAAATTTCTTTCTGATGTCTGCATTCAACTAATA
+
AAFFFKKKFKKFKAA7FKK,F7<AA,<<FAFK7KKFK7F7AF7,AA,77FKKKKKKKKKKKKK,FK,7FA,AK<AAA,<,7,7F7FFKKA,,AK7AFAKK7K,FAK,AAKKKKFFKFFKKKAKA7AFKKK7AKKKKK<,A77<AKK,<AF
@ST-E00169:61:H2CCHCCXX:1:1101:12835:1713 1:N:0:TCCGCGAA
CATGTGAATNGAAAATTATAAGAGAAGTCATTATATGAACAAGACATATGCACACACATGCTTCTAGCAGAACAAGTAAAAATTGCAAAGATATAGACCTGACGTAAGTGCCCATCAAACAGTGGGTGGATAAAGAAAATGTGGTATATA

Second file:
[sohail@datatrans CHG000691]$ cat S17_02A_CHG000691_L001_R2.fastq | head
@ST-E00169:61:H2CCHCCXX:1:1101:7475:1713 2:N:0:TCCGCGAA
ATCCTGAGGAACTTATTTGTGATGTGTGCATTCATCTCAAAGAGTTGAACTTTTCGTTTGATTGCATAGTTTTAAAATACTCTTTTTGTAGAATTTGTATGTGGAAATTTGGGGCGCTTTGCGGCCTTTGGGTGGAATTGAAATATGTTC
+
AA,7FKKKKAKAFF,7FKKKFF7FFFFF,F<F,A<7AF7AFF,F7FKKFFFFFF7,7FFFKK7F,7<,A7K7,,A7<,A<FK<F<K,FKFFA,A,<<,,,,7,7,,,7,,,7,7((((,7<,,,,(,,,,,,,,(,,,,,,,,,,,,,,,
@ST-E00169:61:H2CCHCCXX:1:1101:8815:1713 2:N:0:TCCGCGAA
AAATAGAGTTTCAAAGCTGCTCTGTAAAAAGAAAGGTTCCACTCTGTTAGCTGAGTACACACATCACAAACTTGTTTCTGAGAATCCTTCTGTCTCGTTTTTATGGGAAGATATTCCCTTGTTCACCGTAGGCATCAAAGCGCTTCAAAT
+
AA,7,,A,,7<7AA,7,7FA77FK7FFFFFF,K,<A,,,,F,,AFFAKA,,A,7FA,,A<,FAKF,,,A,F7AFKKK7FAAA,7FK<FF,7<FAF7<AAFFK,<<AA,A<FFKK77,<AF,A,7<F,7A77A,AFA,,A7,,,,,,,,,,
@ST-E00169:61:H2CCHCCXX:1:1101:12835:1713 2:N:0:TCCGCGAA
CTCATAGCTTAGCTCCCACTTATAAGTGAGAACATGTGATATTTGGTTTTCCATTCCTGAGTTACTTCATTTAGAATAATGGTCTCCCATCTCCATTCAAGTTACTGCAAAAGACATTACGTAATTCCTTTCTATGGCAGAGAAATCTTC
[sohail@datatrans CHG000691]$ 

 

ADD REPLYlink written 4.2 years ago by SOHAIL290

Did you ever figure out why this was happening? I'm having the same issue and can't figure out why.

ADD REPLYlink written 19 months ago by joshuau4900
2
gravatar for SOHAIL
19 months ago by
SOHAIL290
Beijing Institute of Genomics, CAS.
SOHAIL290 wrote:

@joshuau490

Hi, Please find the answer in my posted thread at:

https://gatkforums.broadinstitute.org/gatk/discussion/5786/gatk-flagstats-show-0-0-properly-paired-reads

hope it helps.

-sohail

ADD COMMENTlink written 19 months ago by SOHAIL290
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: 906 users visited in the last hour