Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

Close

#279 Bowtie2 2.1.0 missing unmapped reads for --un-conc

v0.9.0
closed
Val
None
5
2014-02-28
2013-02-26
skatzman
No

bowtie2 version 2.1.0 pairedEnd mappings with params:

--end-to-end --no-mixed --no-dovetail --no-contain --no-overlap --no-unal

I specify --un-conc to get unmapped reads, but the resulting fastq files are nearly empty.

Using version 2.0.5 all the expected unmapped reads are present. I have not tried intermediate versions.

I will send Ben examples by private email.

Discussion


  • Anonymous
    2013-02-27

    I got the same error as you noted above using 2.1.0. My analysis also worked fine with v2.0.5. I tested version 2.0.6 and it had the same error as 2.1.0.

     
    Last edit: Anonymous 2013-11-26
  • David Huen
    David Huen
    2013-03-15

    I have also observed what may be similar but not identical behaviour.

    WIth2.1.0, this command:-
    bowtie2 -x DATABASE -1 SRR527876_1.paired.3070.fasta.gz -2 SRR527876_2.paired.3070.fasta.gz -f -p 3 --un-conc one_hit.fa --al-conc paired_hits.fa --no-unal --local --dovetail -S SRR527876.sam

    yields --un-conc files that have the read that MAPPED rather than the one that did not. For example,

    SRR527876.21418 89 sequence1 446 44 11S88M = 446 0 AAACTTTCAATTTAGCACTTAACCTGCCATTTTGCCAATGCGATGTTATGTGAAGGAATTTTTAGTATCTTTGTATTCAAATAATTCACATTTTATGTA IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII AS:i:176 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:88 YT:Z:UP

    shows the first read of a pair (x1) to have mapped (x40) to reverse strand (x10) with mate unmapped (x8).

    In one_hit.1.fa, I see

    >SRR527876.21418 21418/1
    TACATAAAATGTGAATTATTTGAATACAAAGATACTAAAAATTCCTTCACATAACATCGCATTGGCAAAATGGCAGGTTAAGTGCTAAATTGAAAGTTT

    which is the successfully mapped read. There is no corresponding unmapped read in one_hit.2.fa. Could your test result be inverted?

    With 2.0.5, I get a lot of reads back in both files including pairs where both mates were unmapped which doesn't seem right either.

    Have I misunderstood the semantics of --un-conc?

     
  • skatzman
    skatzman
    2013-03-15

    To David Huen -- for pairedEnd mapping, the 2 unmapped fastq output files (_1 and _2) should always have the identical set of IDs. Otherwise they would be useless for the next steps. I think you should double check your output.

    To clarify more precisely the problem with 2.1.0, it seems to be writing to the unmapped fastq files, the set of paired reads that are "aligned discordantly 1 time", rather than (as in 2.0.5) the set of paired reads that are "aligned concordantly 0 times".

    Here is an example. The stats from either 2.0.5 or 2.1.0 are the same:

    72400 reads; of these:
    72400 (100.00%) were paired; of these:
    69262 (95.67%) aligned concordantly 0 times
    2533 (3.50%) aligned concordantly exactly 1 time
    605 (0.84%) aligned concordantly >1 times
    ----
    69262 pairs aligned concordantly 0 times; of these:
    55 (0.08%) aligned discordantly 1 time
    4.41% overall alignment rate

    Here are the counts in the unmapped fastq files for 2.0.5:

    $ gunzip -c hek_hna1_jh01_00of1_trim_unmapped_1.fastq.gz | gawk '(NR % 4 == 2)' | wc -l
    69262

    $ gunzip -c hek_hna1_jh01_00of1_trim_unmapped_2.fastq.gz | gawk '(NR % 4 == 2)' | wc -l
    69262

    Here are the counts in the unmapped fastq files for 2.1.0

    $ gunzip -c hek_hna1_jh01_00of1_trim_unmapped_1.fastq.gz | gawk '(NR % 4 == 2)' | wc -l
    55

    $ gunzip -c hek_hna1_jh01_00of1_trim_unmapped_2.fastq.gz | gawk '(NR % 4 == 2)' | wc -l
    55

     
  • Brian Walenz
    Brian Walenz
    2013-05-07

    The problem is that --no-unal is overriding --un-conc.

    Line 426 in bowtie2 (from 2.0.6) is the change that broke it. I kind-of see a fix -- test if ($filt) && !("the test on line 448") -- but I haven't quite figured out all the cases in this area.

     
  • Val
    Val
    2014-02-21

    • assigned_to: Val
    • Group: --> v0.9.0
     
  • Val
    Val
    2014-02-28

    • status: open --> closed
     
  • Val
    Val
    2014-02-28

    Fixed in 2.2.1