#242 reference sequence dependent mis-assembly

open
nobody
None
5
2012-12-14
2012-09-17
No

I am getting a strange assembly artifact with bowtie2 beta 6.
I am assembling HiSeq and MiSeq reads of avian influenza A virus samples.
bird flu consists of 8 segments that total 13Kb. My reference sequence consists of :
segment 1: 2 sequences.
segment 2: 2 sequences.
segment 3: 3 sequences.
segment 5: 3 sequences.
segment 7: 2 sequences.
segment 8: 3 sequences.
segment 4: 20 sequences.
segment 6: 9 sequences.

The last two segments, segment 4 and segment 6, are highly divergent, and thus require large number of reference sequences. The sequences are present in the reference file in that order.

When I assemble with this set, I get correct assembly for the first 6 segments, 1,2,3,5,7,8, and sometimes correct segment4 and most often errounous assembly of segment 6. It maps significant number of reads to one of the segment 6 reference. This however is not the correct one.
If I assemble just to the segment 6 references, the same reads assemble correctly. The sequences in segments 4 and segment 6 are highly divergent (40% amino acids identity, no nucleotide level similarity) and reads should not get sucked into a similar looking reference.
At first I thought the program was running out of a physical resources but this strange behaviour happens for Miseq data as well.

Discussion

  • Ben Langmead

    Ben Langmead - 2012-12-14

    Dear Shin,

    Can you give a specific example where you think Bowtie 2 is aligning a read incorrectly? It's hard to tell from your describing whether Bowtie 2 is doing something wrong, or whether it's just doing something other than what you expected.

    Best,
    Ben

     
  • Ben Langmead

    Ben Langmead - 2012-12-14
    • status: open --> pending
     
  • Shin Enomoto

    Shin Enomoto - 2012-12-14
    • status: pending --> open
     
  • Shin Enomoto

    Shin Enomoto - 2012-12-14

    Initially when I assemble against the entire set of references I get mapping results like below saying most of the sequences mapped to seg6N6

    gi|190352319|gb|CY032666|seg6N1 1458 6 4
    gi|241828973|gb|GQ257383|seg6N2 1442 33 1
    gi|153805739|gb|EU030978|seg6N3 1410 0 0
    gi|190194859|gb|EU557562|seg6N4 1413 0 0
    gi|110733396|gb|CY012842|seg6N5 1435 0 0
    gi|229598099|gb|CY039772.1|seg6N6 1464 149503 144
    gi|241828923|gb|GQ257390|seg6N7 1436 0 0
    gi|190682238|gb|CY032894|seg6N8 1460 50 3
    gi|241828297|gb|GQ257481|seg6N9 1435 0 0
    gi|343971488|gb|CY097195.1|seg4H1 1746 0 0
    H2consensus|seg4H2 1788 64 12
    gi|390517097|emb|HE802063.1|seg4H5 1767 0 0
    H5usaConsens|seg4H5 1800 10 4
    H9consensus|seg4H9 1742 0 0
    H14consensus|seg4H14 1749 0 0
    H15consensus|seg4H15 1764 165 29

    However the consensus sequence extracted looks horrible and sometimes incorrect.

    @gi|229598099|gb|CY039772.1|seg6N6
    AGCRAAAGCAGGGTGAAAATGAATCCAAATCAGAAGATAATATGCATCTCAGCAACAGGA
    ATGACACTATCCGTAGTAAGTCTGCTAATAGGATTGGCCAACTTGGGGCTGAACATTGGG
    CTTCATTTCAAGGTAGGAGACACACCGGAAAggaccagTGCTAGCACCAATGAGACctgt
    ggtgttaattctgataccacaggttggtCatggccCAACTTCRCTcTgtTGACCTttGAC
    atagacaagtaangnnccAATATAATGTTcnacgnngnnngnnnttaaaaatgaatccaa
    atcagaaaataataacaattggttctgtgtcattggcactagttgtattcaacatACtgc
    ttcatattgtatcAatagtcataggaATaatatcagtgacaaaagaaagcaGTatatCat
    caacctGCAACCCCACTACGGTGTacaatgaaactgtaaggctggaaactataacaattc
    ctatcaaCaacactgtttacatAGAaagagagtcacatcaagaacctgAGtttTtaAaca
    atacagaacctcTctgcaatgtatcCGggttcgcaatagtttCCAAGGAcaatggaatca
    GAaTtGGGACAaGgTAACATGTCTTCGGCAGTGGAgaaccattcgtggcGtgTGGTccca
    cagAaTGTaggacttttttcctaacgcaaggtgCCTCActgaatgataaacactccaaca
    atacagtgaaagataGGAGTGATGACGGtgcattgatgagtgtgccattagGATcttcac
    cCAATGcttaTCAGgccaagtttgaGACAGttgcatGGtcggccacagcatgtcatgatg
    gcaAagagtggctggcaGTagggataAGCGGtgcAGatgAcgatgcttATgctgtAAGCC
    ATTatgggggggtgccaacagatgtggtgaggtcatggagaaAgcaaattctaagaaCac
    aagaatcatcatgtgtgtgtatgaaaggaaactgtTAGTCGGtaatgacggatggtcctg
    cgaacaatcaggctAgttACAAaatTTTCAAGTCtcataagggaatggtgacaaacgaaA
    gagaagtgtcatttcagggaggtcacattgaaGAgtgctcttgctacCTCAAGGTgggta
    aagtGGAgtgcgtttGCCCAGAtaattGGAACRGAAtgaatagaccaattttgacttttg
    atGATAActtgaactatgaggtgggttatttatgtgcCggaataccgacagacactccac
    gggttCaagaCAgtggtttcactggttccTgtactaatgCTGATCgagggagtggaacAA
    acAAGTAGGGagTTATGGGtTttggcttcagacaaggtaatagtgtgtgggCAGGAGGAA
    CTGTTAGCATttcatcccgAaGtgggtttgagatcctattgaTagAAgacggttggacca
    aaacaagcaaaaCTgtcgTCaaaa
    +

    ~~~~~~~~~~~~~~~~~~~~~~~~~~HHHHB?????????E~~~~~~~~~~~~~~H????
    ????????????????????????????~??????~~~NEEf$W?????K~~~~??????
    ??????????????????~~T:?HKQ]NH???????????????????????????????
    ???????????????????????????????????????????????????????~~???
    ?????????????B????????????B?????????????????????????B???????
    ??????`B0c`KEc~~~:NZZBHW????????????????????????????????????
    ???????B????????????????????????????????????????????????????
    ????????????B??????????????????????????????~~~~(N???????????
    ?????>~~~~????BB~~~B?BB~~R~~~~H=WWB??????????????????BQ~????
    ?????~~???????????????????????????~~Hc??????????????????????
    ???????????????K?B]~~R~~~~i~???????????????????????QN???????
    ?B~~~?????~~{????????????BK~~~??????BB??????????????????????
    ?????????????????BB???????W~>~Q?????????????????KE??????HN,~
    ~?c?????????????????????????????????????????????????????????
    ????????????????????????????????????~~~~~~??????????????????
    ??????????????????EBEE?????B?B~~?f?????????????????????????K
    ?????????????????????????????????H??????????????K~~~~F~?????
    ????BBf????????T?NB?Q???????BQNH*~~?????????????????????????
    ??~xK~~?????????????????????????????????????????????????????
    ???????????N?????????????????f??????????~~~~H??????????????B
    ??B~~~~WTZ??KEBEQT>?????????????????????????????????~~~~C?BB
    ?EH??NB]?E???????????B????????????????????B??BB?????????????
    ?????????????~??????????

    If instead I assemble to just the segN1-segN9 references, the mapping looks ok.

    gi|190352319|gb|CY032666|seg6N1 1458 6 4
    gi|241828973|gb|GQ257383|seg6N2 1442 33 1
    gi|153805739|gb|EU030978|seg6N3 1410 0 0
    gi|190194859|gb|EU557562|seg6N4 1413 0 0
    gi|110733396|gb|CY012842|seg6N5 1435 0 0
    gi|229598099|gb|CY039772.1|seg6N6 1464 149719 391
    gi|241828923|gb|GQ257390|seg6N7 1436 0 0
    gi|190682238|gb|CY032894|seg6N8 1460 50 4
    gi|241828297|gb|GQ257481|seg6N9 1435 0 0

    And the consensus sequence also appears to be high quality.

    @gi|229598099|gb|CY039772.1|seg6N6
    AGCRAAAGCAGGGTGAAAATGAATCCAAATCAGAAGATAATATGCATCTCAGCAACAGGA
    ATGACACTATCCGTAGTAAGTCTGCTAATAGGATTGGCCAACTTGGGGCTGAACATTGGG
    CTTCATTTCAAGGTAGGAGACACACCGGAAACAGGAACCCCTAGCACCAATGAGACAAAC
    TCCACAACCACAATAATTAACTACAACACCCAAAACAACTTCACAAATGTGACCAACATC
    GTGTTAGTTAAAGAAGAAAATAAAATGTTTACAAACCTTTCTAAGCCCTTGTGTGAAGTA
    AACTCATGGCATATTCTATCCAAGGACAATGCGATTAGAATAGGAGAGGATGCCCACATC
    CTTGTCACAAGAGAACCATATCTCTCATGTGGACCACATGAGTGCAGAATGTTTGCCCTC
    AGCCAAGGTACCACACTACGGGGTCGGCATGCAAATGGGACTATACATGACAGAAGCCAG
    TTTAGGGCATTAATAAGCTGGGAAATGGGGCAAGCGCCGAGTCCGTACAATACCAGAGTA
    GAATGTGTGGGATGGTCCAGCACTTCATGCCATGACGGCATCTCAAGAATGTCAATCTGC
    ATGTCAGGACCTAATAACAATGCTTCGGCAGTGGTCTGGTACAATGGAAGGCCAGTAACC
    GAAATTGCTTCATGGACAGGAAATATATTAAGGACTCAGGAATCAGAATGTGTATGCCAT
    AATGGAATATGCCCTGTGGTGATGACGGATGGCCCAGCTAATAACAGAGCAGAGACAAAA
    ATAATTTACTTCAAAGAAGGAAAAATACAGAAAATAGGGGAATTGACAGGAAGTGCACAG
    CATATAGAAGAATGCTCGTGCTATGGAGCGGAAGAAGTAATTAAATGCATTTGCAGAGAC
    AATTGGAAAGGTGCAAATAGACCAGTAATTACCATAAACCCAACAGCAATGACCCATACG
    AGCAAATACTTGTGTTCAAAGATTCTAACTGACACAAGTCGGCCTAATGACCCCGGAAGC
    GGGAACTGTGATGCACCAATAACCGGAGGGAGCCCAGATCCTGGCGTAAAAGGATTTGCA
    TTCTTAGATGGGGGGAATTCCTGGTTGGGAAGGACCATAAGCAAAGATTCAAGGTCAGGG
    TATGAGATGCTAAAAGTCCCGAATGCGGAAACAGACAATCAGTCCGCTCCAGTTGCACAT
    CAGATAATAGTGAACAACCAAAACTGGTCAGGATACTCAGGAGCGTTCATCGATTATTGG
    GCTGATAGAGAATGCTTCAATCCTTGCTTTTATGTGGAATTGATCAGAGGCAGGCCAAAA
    GAGAGTAGTGTATTGTGGACATCTAATAGTATAGTAGCGCTCTGTGGATCCAAGGAGCGA
    TTGGGATCTTGGTCATGGCATGATGGGGCTGAAATCATCTACTTCAAGTAGACGAGATTT
    TGAAAAAACACCCTTGTTTCTACT
    +

    ~~~~~~~~~~~~~~~~~~~~~~~~

     
  • Shin Enomoto

    Shin Enomoto - 2013-02-03

    I looked at the alignment with samtools tview. Two bam files, first one (sort0.bam) was created using 43 reference sequences and the second one (H1-16.sort0.bam) was created by 21 reference sequences. Below is a tview of sort0.bam that mapped strange:
    61 71 81 91 101 111 121 131 141 151 161 171 181 191
    GCAGGGGTCACAATGGAGAAGTTTATCGCAA*TAGCAATGCTCTTG*****G*******TGAGCACAAATGCATACGATAGGATATGCATTGGTTACCAATCGAACAACTCCACAGATACAGTGAACACTCTTATAGAACAGAATGTGCCGGTCACTCAA
    ...T..TCTTG.CAAA.T.T.CC.TGGA..C C.TGC.GT.C.CAA T G.CA.CATCG.AA.G..TC.CAC.A.C.A.AA..AAGAGGTGACC.ATG..A.TG.A.CGGT.GA..GCAAAAGCCT.G...A.C.T...AT.AA.AA.CG
    C.CTC.TCTTGTCAAATTCTCCC.TGGC.TC*CCTGCCCT.C.CCT*****T*******C.CC.C.TCG.T..G..TC.CCCCACC aa,,aagaggtgacc,atg,,a,tg,a,cggt,ga,,gcaaaagcct,g,,,a,c,t,,,at,aa,aa,cg
    ...T..TCTTG.CAAA.T.T.CC.TGGA..C*C.TGC.GT.C.CAA*****T*******G.CA.CATCG.AA.G..TC.CAC.A.C.A.AA..A gaggtgacc,atg,,a,tg,a,cggt,ga,,gcaaaagcct,g,,,a,c,t,,,at,aa,aa,cg
    ...T..TCTTGGCAAA.T.T.CC.TGGA..C*C.TGC.GT.C.CAA*****T*******G.CA.CATCG.AA.G..CC.CAC.A.C.A.AA..A gaggtgacc,atg,,a,tg,a,cggt,ga,,gcaaaagcct,g,,,a,c,t,,,ataaa,aa,cg
    ...T..TCTTG.CAAA.T.T.CC.TGGA..C*C.TGC.GT.C.CAA*****T*******G.CA.CATCG.AA.G..TC.CAC.A.C.A.AA..A gaggtgacc,atg,,a,tg,a,cggt,ga,,gcaaaagcct,g,,,a,c,t,,,at,aa,aa,cg
    ...T..TCTTG.CAAA.T.T.CC.TGGA..C*C.TGC.GT.C.CAA*****T*******G.CA.CATCG.AA.G..TC.CAC.A.C.A.AA..A gaggtgacc,atg,,a,tg,a,cggt,ga,,gcaaaagcct,g,,,a,c,t,,,at,aa,aa,cg
    ......CCTTG.CAAA.T.T.CC.TGGA..C*C.TGC.GT.C.CAA*****T*******G.CC.CATCG.AA.G..TC.CAC.A.C.A.AA..A gaggtgacc,atg,,a,tg,a,cggt,ga,,gcaaaagcct,g,,,a,c,t,,,at,aa,aa,cg
    ...T..TCTTG.CAAA.T.T.CC.TGGA..C*C.TGC.GT.C.CAA*****T*******G.CA.CATCG.AA.G..TC.CAC.A.C.A.AA..A gaggtgacc,atg,,a,tg,a,cggt,ga,,gcaaaagcct,g,,,a,c,t,,,at,aa,aa,cg

    The reference sequence is not the correct one.

    More details here: http://www.sriflu.net/shin2/

     
  • Shin Enomoto

    Shin Enomoto - 2013-03-05

    The same problem happens on bam files created by Mosaik. I don't think this is a bowtie2 problem.

     

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks