{"action":"created","issue":{"url":"https://api.github.com/repos/samtools/bcftools/issues/306","labels_url":"https://api.github.com/repos/samtools/bcftools/issues/306/labels{/name}","comments_url":"https://api.github.com/repos/samtools/bcftools/issues/306/comments","events_url":"https://api.github.com/repos/samtools/bcftools/issues/306/events","html_url":"https://github.com/samtools/bcftools/issues/306","id":102039731,"number":306,"title":"bcftools consensus bug","user":{"login":"bioprog","id":13045632,"avatar_url":"https://avatars.githubusercontent.com/u/13045632?v=3","gravatar_id":"","url":"https://api.github.com/users/bioprog","html_url":"https://github.com/bioprog","followers_url":"https://api.github.com/users/bioprog/followers","following_url":"https://api.github.com/users/bioprog/following{/other_user}","gists_url":"https://api.github.com/users/bioprog/gists{/gist_id}","starred_url":"https://api.github.com/users/bioprog/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/bioprog/subscriptions","organizations_url":"https://api.github.com/users/bioprog/orgs","repos_url":"https://api.github.com/users/bioprog/repos","events_url":"https://api.github.com/users/bioprog/events{/privacy}","received_events_url":"https://api.github.com/users/bioprog/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":2,"created_at":"2015-08-20T01:58:51Z","updated_at":"2015-08-26T00:33:51Z","closed_at":null,"body":"Hi,\r\nI'm pretty sure this is a bug resulting from something like a for loop in consensus.c where something is not checked/done after the last iteration.\r\n\r\nI have a reference sequence that ends with a couple of \"N\"s. Exactly 1 read covers ALL sites of that N region (and some sites before the N region). No other read reaches into that N region.\r\n\r\nNote: The situation was set up by me with test data on purpose to do some checks on bcftools' consensus behaviour.\r\n\r\nWhen you call the consensus, bcftools fills the sites within the N region with the the nucleotide of that 1 read which covers it.\r\n\r\nHowever at the LAST position it puts an N into the consensus instead of the nucleotide from the read (a T in my case).\r\nWhen I check the output of 'bcftools call' via 'bcftools view', I see that the scores are all the same in that N region sites (as expected because they are all subject to the same conditions [i.e. 1 read covering the site ... oh and the quality of every base is the same in my test data]). It's only that last position which gets assigned an N in the consensus instead of the read base."},"comment":{"url":"https://api.github.com/repos/samtools/bcftools/issues/comments/134778718","html_url":"https://github.com/samtools/bcftools/issues/306#issuecomment-134778718","issue_url":"https://api.github.com/repos/samtools/bcftools/issues/306","id":134778718,"user":{"login":"bioprog","id":13045632,"avatar_url":"https://avatars.githubusercontent.com/u/13045632?v=3","gravatar_id":"","url":"https://api.github.com/users/bioprog","html_url":"https://github.com/bioprog","followers_url":"https://api.github.com/users/bioprog/followers","following_url":"https://api.github.com/users/bioprog/following{/other_user}","gists_url":"https://api.github.com/users/bioprog/gists{/gist_id}","starred_url":"https://api.github.com/users/bioprog/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/bioprog/subscriptions","organizations_url":"https://api.github.com/users/bioprog/orgs","repos_url":"https://api.github.com/users/bioprog/repos","events_url":"https://api.github.com/users/bioprog/events{/privacy}","received_events_url":"https://api.github.com/users/bioprog/received_events","type":"User","site_admin":false},"created_at":"2015-08-26T00:33:51Z","updated_at":"2015-08-26T00:33:51Z","body":"Alright here's an example:\r\n\r\n$ cat tmp_ref.fasta\r\n>myseq_contig1\r\nNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN\r\nNNNNNNNNNNCCTCGGAAACTCTTAATGGCACATAAACGACCGAATCAATGACTCACCTT\r\nTGTCCGAGTCTCTAACCGTTTAAGGTGAATTAACATACCGGGGAGGTGAATCCAAGGTCT\r\nCAGACCACTTCGGACCGAGATGCTTCGCATNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN\r\nNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN\r\n\r\n$ cat reads.fastq\r\n@left_overhang\r\nACGTACGACTGACTGAGCGGGCGCATCTATTGCATTCGATCGGAACGATCGATCGACAGCGACTAGCATCCCTCGGAAACTCTTAATGGCACATAAACGA\r\n+\r\nJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ\r\n@read71_171\r\nCCTCGGAAACTCTTAATGGCACATAAACGACCGAATCAATGACTCACCTTTGTCCGAGTCTCTAACCGTTTAAGGTGAATTAACATACCGGGGAGGTGAA\r\n+\r\nJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ\r\n@read91_191\r\nACATAAACGACCGAATCAATGACTCACCTTTGTCCGAGTCTCTAACCGTTTAAGGTGAATTAACATACCGGGGAGGTGAATCCAAGGTCTCAGACCACTT\r\n+\r\nJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ\r\n@read111_211\r\nGACTCACCTTTGTCCGAGTCTCTAACCGTTTAAGGTGAATTAACATACCGGGGAGGTGAATCCAAGGTCTCAGACCACTTCGGACCGAGATGCTTCGCAT\r\n+\r\nJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ\r\n@right_overhang\r\nTCCAAGGTCTCAGACCACTTCGGACCGAGATGCTTCGCATCCGAGTATCTCGGCAAAGTTATAGGCGTGGGGGCTACGCTCCAAGGCGTTATTTACTTAT\r\n+\r\nJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ\r\n\r\nThe 3 reads in the middle match the reference with 100% identity at the positions indicated (the read \"@read71_171\" matches from 71 to 171 ... etc ...).\r\n71 is where the non-N-region starts. 211 (name of the 3rd read: \"@read111_211\") is where the non-N-region ends.\r\nThus those 3 reads span the non-N-region.\r\n\r\nNow to the 2 \"overhang\" reads:\r\n\r\nThe LAST 30 bases of the read \"@left_overhang\" match the FIRST 30 bases of \"@read71_171\" (thus the first 30 bases of the non-N-region). So the first 70 bases of that read overlap the N-region at the start of the reference.\r\n\r\nThe FIRST 40 bases of \"@right_overhang\" match the LAST 40 bases of the non-N-region. So the last 60 bases of that read overlap the N-region at the end.\r\n\r\nNow the commands (You can just copy them into a script and execute it, it will show you the consensus in the terminal [assuming you have copied the reads and the reference into the according filenames beforehand]):\r\nbwa index tmp_ref.fasta\r\nbwa mem tmp_ref.fasta reads.fastq > tmp.sam\r\nsamtools view -h tmp.sam | samtools mpileup -B -vf tmp_ref.fasta - | bcftools call -Mm -Oz -o calls.vcf.gz -\r\nbcftools tabix -f calls.vcf.gz\r\nbcftools consensus -f tmp_ref.fasta calls.vcf.gz\r\n\r\nThis results in a consensus sequence that shows you \"N\"s in the regions, where I want to have the bases of the overhangs.\r\nBut the overhangs get SOFT-CLIPPED (you can see it in the CIGAR string in the tmp.sam file).\r\n\r\nSo I just modify that tmp.sam file via a script and do the following:\r\nI replace the CIGAR strings from 70S30M (\"@left_overhang\") and 40M60S (\"@right_overhang\") to 100M. Furthermore I change the start position (column 4 in tmp.sam) of \"@left_overhang\" from 71 to 1 (The start of \"@right_overhang\" stays the same).\r\n\r\nThe resulting tmp_edited.sam file looks like this:\r\n\r\n@SQ\tSN:myseq_contig1\tLN:270\r\n@PG\tID:bwa\tPN:bwa\tVN:0.7.12-r1039\tCL:bwa mem tmp_ref.fasta reads.fastq\r\nleft_overhang\t0\tmyseq_contig1\t1\t60\t100M\t*\t0\t0\tACGTACGACTGACTGAGCGGGCGCATCTATTGCATTCGATCGGAACGATCGATCGACAGCGACTAGCATCCCTCGGAAACTCTTAATGGCACATAAACGA\tJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ\tNM:i:0\tMD:Z:30\tAS:i:30\tXS:i:0\r\nread1_101\t0\tmyseq_contig1\t71\t60\t100M\t*\t0\t0\tCCTCGGAAACTCTTAATGGCACATAAACGACCGAATCAATGACTCACCTTTGTCCGAGTCTCTAACCGTTTAAGGTGAATTAACATACCGGGGAGGTGAA\tJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ\tNM:i:0\tMD:Z:100\tAS:i:100\tXS:i:0\r\nread21_121\t0\tmyseq_contig1\t91\t60\t100M\t*\t0\t0\tACATAAACGACCGAATCAATGACTCACCTTTGTCCGAGTCTCTAACCGTTTAAGGTGAATTAACATACCGGGGAGGTGAATCCAAGGTCTCAGACCACTT\tJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ\tNM:i:0\tMD:Z:100\tAS:i:100\tXS:i:0\r\nread41_141\t0\tmyseq_contig1\t111\t60\t100M\t*\t0\t0\tGACTCACCTTTGTCCGAGTCTCTAACCGTTTAAGGTGAATTAACATACCGGGGAGGTGAATCCAAGGTCTCAGACCACTTCGGACCGAGATGCTTCGCAT\tJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ\tNM:i:0\tMD:Z:100\tAS:i:100\tXS:i:0\r\nright_overhang\t0\tmyseq_contig1\t171\t60\t100M\t*\t0\t0\tTCCAAGGTCTCAGACCACTTCGGACCGAGATGCTTCGCATCCGAGTATCTCGGCAAAGTTATAGGCGTGGGGGCTACGCTCCAAGGCGTTATTTACTTAT\tJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ\tNM:i:0\tMD:Z:40\tAS:i:40\tXS:i:0\r\n\r\nIf I call the consensus pipeline as above (only replacing \"samtools view -h tmp.sam\" by \"samtools view -h tmp_edited.sam\") I get the desired result: The bases in the N-region are the bases of the read. EXCEPT for the last position."}}