{"action":"created","issue":{"url":"https://api.github.com/repos/FelixKrueger/Bismark/issues/170","repository_url":"https://api.github.com/repos/FelixKrueger/Bismark","labels_url":"https://api.github.com/repos/FelixKrueger/Bismark/issues/170/labels{/name}","comments_url":"https://api.github.com/repos/FelixKrueger/Bismark/issues/170/comments","events_url":"https://api.github.com/repos/FelixKrueger/Bismark/issues/170/events","html_url":"https://github.com/FelixKrueger/Bismark/issues/170","id":318068606,"number":170,"title":"Sorting bismark BAM files by read ID and problems with methylation extractor","user":{"login":"duncansproul","id":38759321,"avatar_url":"https://avatars2.githubusercontent.com/u/38759321?v=4","gravatar_id":"","url":"https://api.github.com/users/duncansproul","html_url":"https://github.com/duncansproul","followers_url":"https://api.github.com/users/duncansproul/followers","following_url":"https://api.github.com/users/duncansproul/following{/other_user}","gists_url":"https://api.github.com/users/duncansproul/gists{/gist_id}","starred_url":"https://api.github.com/users/duncansproul/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/duncansproul/subscriptions","organizations_url":"https://api.github.com/users/duncansproul/orgs","repos_url":"https://api.github.com/users/duncansproul/repos","events_url":"https://api.github.com/users/duncansproul/events{/privacy}","received_events_url":"https://api.github.com/users/duncansproul/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":1,"created_at":"2018-04-26T14:59:10Z","updated_at":"2018-04-27T14:07:05Z","closed_at":null,"author_association":"NONE","body":"Hi Felix,\r\n\r\nAs I mentioned by e-mail a few weeks ago, I think I’ve come across an issue with sorting Bismark BAM files that can cause problems with downstream analysis. Essentially, methylation calls for paired-end data can get misallocated from OB to CTOB when you sort a Bismark BAM file by position before resorting by read name.\r\n\r\nI’ve done a few more tests and read the methylation extractor code and confirmed that the problem does happen. Bizarrely, it only happens using Sambamba to sort and it is not seen using Samtools. I thought it was worth reporting it in case others come across it and also because it potentially has a fairly easy solution that could be incorporated into Bismark in future.\r\n\r\nAs I understand it, the methylation extractor module determines what strand data come from by comparing the XR and XG codes for the first read in each pair (read and genome conversion fields respectively). I think the following segment of code deals with this (this is code ~Line 2602):\r\n\r\n```perl\r\nif ($first_read_conversion eq 'CT' and $genome_conversion eq 'CT') {\r\n\t  $index = 0;\t\t## this is OT\r\n\t  $strand = '+';\r\n\t} elsif ($first_read_conversion eq 'GA' and $genome_conversion eq 'CT') {\r\n\t  $index = 1;\t\t## this is CTOT\r\n\t  $strand = '-';\r\n\t} elsif ($first_read_conversion eq 'GA' and $genome_conversion eq 'GA') {\r\n\t  $index = 2;\t\t## this is CTOB\r\n\t  $strand = '+';\r\n\t} elsif ($first_read_conversion eq 'CT' and $genome_conversion eq 'GA') {\r\n\t  $index = 3;\t\t## this is OB\r\n\t  $strand = '-';\r\n\t} else {\r\n\t  die \"Unexpected combination of read and genome conversion: $first_read_conversion / $genome_conversion\\n\";\r\n\t}\r\n```\r\n\r\nAfter bismark alignment all read pairs are arranged in the BAM with R1 followed by R2 and the methylation extractor assumes this. Both reads are also given exactly the same ID during alignment, ie the part indicating R1 or R2 is the same (they are both labelled R1). Here is some toy output from a short BAM (just showing the important fields with two pairs of reads):\r\n\r\n```\r\nE00397:87:HHMVMCCXY:6:1101:7720:1221_1:N:0:27\t99\tchr1\t106372402\tXR:Z:CT\tXG:Z:CT\r\nE00397:87:HHMVMCCXY:6:1101:7720:1221_1:N:0:27\t147\tchr1\t106372475\tXR:Z:GA\tXG:Z:CT\r\nE00397:87:HHMVMCCXY:6:1101:13930:1221_1:N:0:27\t83\tchr12\t45034123\tXR:Z:CT\tXG:Z:GA\r\nE00397:87:HHMVMCCXY:6:1101:13930:1221_1:N:0:27\t163\tchr12\t45033911\tXR:Z:GA\tXG:Z:GA\r\n```\r\n\r\nAfter positional sorting, pairs aligned to OB are reversed because R2 has a position lower than R1 (OT aligned reads remain the same):\r\n\r\n```\r\nE00397:87:HHMVMCCXY:6:1101:7720:1221_1:N:0:27\t99\tchr1\t106372402\tXR:Z:CT\tXG:Z:CT\r\nE00397:87:HHMVMCCXY:6:1101:7720:1221_1:N:0:27\t147\tchr1\t106372475\tXR:Z:GA\tXG:Z:CT\r\nE00397:87:HHMVMCCXY:6:1101:13930:1221_1:N:0:27\t163\tchr12\t45033911\tXR:Z:GA\tXG:Z:GA\r\nE00397:87:HHMVMCCXY:6:1101:13930:1221_1:N:0:27\t83\tchr12\t45034123\tXR:Z:CT\tXG:Z:GA\r\n```\r\n\r\nAfter sorting the position-sorted file by read name, this inverse ordering of OB aligned pairs persists. This is probably because both R1 and R2 are assigned exactly the same ID during bismark alignment.\r\n\r\n```\r\nE00397:87:HHMVMCCXY:6:1101:7720:1221_1:N:0:27\t99\tchr1\t106372402\tXR:Z:CT\tXG:Z:CT\r\nE00397:87:HHMVMCCXY:6:1101:7720:1221_1:N:0:27\t147\tchr1\t106372475\tXR:Z:GA\tXG:Z:CT\r\nE00397:87:HHMVMCCXY:6:1101:13930:1221_1:N:0:27\t163\tchr12\t45033911\tXR:Z:GA\tXG:Z:GA\r\nE00397:87:HHMVMCCXY:6:1101:13930:1221_1:N:0:27\t83\tchr12\t45034123\tXR:Z:CT\tXG:Z:GA\r\n```\r\n\r\nAs I mentioned, this appears to be specific to Sambamba. Repeating a similar double sort with Samtools reverts to the correct order. This implies that Samtools sort also looks at something other than the read name field when sorting by read ID.\r\n\r\nThe issue is not picked up by the “Now testing Bismark result file $filename for positional sorting (which would be bad...)” part of the program because it tests:\r\ni) that the BAM file does not contain the positional sorting flag.\r\nii) that the first 100,000 reads are paired by ID (which they are).\r\n\r\nThe problem can be recapitulated using the attached short BAM file and the following commands (assuming sambamba and samtools are installed):\r\n[TEST_bismark.bam.zip](https://github.com/FelixKrueger/Bismark/files/1951697/TEST_bismark.bam.zip)\r\n\r\n```bash\r\n# Original file\r\nsamtools view TEST_bismark.bam | cut -f 1,2,3,4,15,16\r\n\r\n# Sort by position then by Read ID with SAMBAMBA\r\nsambamba sort -o TEST_bismark_sortPos_SAMBAMBA.bam TEST_bismark.bam\r\nsambamba sort -n -o TEST_bismark_sortID_SAMBAMBA.bam TEST_bismark_sortPos_SAMBAMBA.bam\r\n\r\n# View results:\r\n# Sort by position - SAMBAMBA\r\nsamtools view TEST_bismark_sortPos_SAMBAMBA.bam | cut -f 1,2,3,4,15,16\r\n# Resort by read ID - SAMBAMBA\r\nsamtools view TEST_bismark_sortID_SAMBAMBA.bam | cut -f 1,2,3,4,15,16\r\n\r\nsamtools view TESTshort_bismark_sortPos.bam | cut -f 1,2,3,4,15,16\r\nsamtools view TESTshort_bismark_sortID.bam | cut -f 1,2,3,4,15,16\r\n\r\n# Now repeat with Samtools:\r\n# Sort by position then name\r\nsamtools sort TEST_bismark.bam TESTshort_bismark_sortPos_SAMTOOLS\r\nsamtools sort -n TESTshort_bismark_sortPos.bam TEST_bismark_sortID_SAMTOOLS\r\n\r\n# View results\r\n# Sort by position - Samtools\r\nsamtools view TESTshort_bismark_sortPos_SAMTOOLS.bam | cut -f 1,2,3,4,15,16\r\n# Resort by ID - Samtools\r\nsamtools view TESTshort_bismark_sortID_SAMTOOLs.bam | cut -f 1,2,3,4,15,16\r\n```\r\n\r\nThe two solutions are:\r\ni)\tDon’t use Sambamba for sorting if you are using Bismark.\r\nii)\tIn the future, R1/R2 IDs could be preserved after Bismark alignment. I did a quick test by manually editing the small bismark BAM file and it removed the problem. Appending a new _1 or _2 worked also. \r\n\r\nFor completeness, I should say I’ve seen this problem with Bismark v0.16.3 and v0.18.1 but have not yet tested on the latest version. I’m also using samtools v1.2 and sambamba v0.6.6.\r\n\r\nThanks\r\n\r\nDuncan"},"comment":{"url":"https://api.github.com/repos/FelixKrueger/Bismark/issues/comments/384980677","html_url":"https://github.com/FelixKrueger/Bismark/issues/170#issuecomment-384980677","issue_url":"https://api.github.com/repos/FelixKrueger/Bismark/issues/170","id":384980677,"user":{"login":"FelixKrueger","id":5870509,"avatar_url":"https://avatars2.githubusercontent.com/u/5870509?v=4","gravatar_id":"","url":"https://api.github.com/users/FelixKrueger","html_url":"https://github.com/FelixKrueger","followers_url":"https://api.github.com/users/FelixKrueger/followers","following_url":"https://api.github.com/users/FelixKrueger/following{/other_user}","gists_url":"https://api.github.com/users/FelixKrueger/gists{/gist_id}","starred_url":"https://api.github.com/users/FelixKrueger/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/FelixKrueger/subscriptions","organizations_url":"https://api.github.com/users/FelixKrueger/orgs","repos_url":"https://api.github.com/users/FelixKrueger/repos","events_url":"https://api.github.com/users/FelixKrueger/events{/privacy}","received_events_url":"https://api.github.com/users/FelixKrueger/received_events","type":"User","site_admin":false},"created_at":"2018-04-27T14:07:05Z","updated_at":"2018-04-27T14:07:05Z","author_association":"OWNER","body":"Hi Duncan,\r\n\r\nMany thanks again for providing such an excellent issue report, including test data and a detailed step-by-step procedure. I have now looked at the issue in more detail, and I am coming to the conclusion that the problem here is probably more an unwanted behaviour of `sambamba` rather than something I would like to address explicitly within Bismark.\r\n\r\nAs you stated, re-sorting a positionally sorted BAM file by name with `samtools` brings R1 and R2 into the correct order again. For the sorting-by-name procedure, I believe that Samtools does not only use the readID on its own, but uses the FLAG value in addition to that:\r\n\r\n>**SAM flag 83**:\r\n   read paired\r\n   read mapped in proper pair\r\n   read reverse strand\r\n   **first in pair**\r\n\r\n>**SAM flag 163**:\r\n   read paired\r\n   read mapped in proper pair\r\n   mate reverse strand\r\n   **second in pair**\r\n\r\nIn contrast, `sambamba` doesn't seem to be doing that. So we have 'encoded' the read identity into the SAM FLAG already if you will, albeit not explicitly via appending R1 or R2 in the read ID field. \r\n\r\nIn a nutshell I would like to settle for your first proposed solution:\r\ni)\tDon’t use `Sambamba` for sorting if you are using Bismark.\r\n\r\nor maybe: \"If you absolutely have to sort Bismark BAM files by chromosomal position between the mapping and deduplication and/or methylation extraction procedure (which we definitely do not advise you should do) you may use `samtools` or `sambamba` for the positional sorting; if you want to re-sort these files by read name you have to be using `samtools sort -n` but not `sambamba`.\r\n\r\nAnd as a final note: the corner-case issue of positionally sorting and then re-sorting by name that aforementioned positionally sorted file with `sambamba` discussed here may result in incorrectly assigning the strand from OB to CTOB. This should however not make any difference for the actual methylation call itself, because both OB and CTOB report on the very same C positions on the bottom strand. So it really is rather a cosmetic issue affecting read orientations anyway.\r\n\r\nDoes this sound like an acceptable solution? Cheers, Felix"}}
{"actor":{"display_login":"FelixKrueger"}}