Hi Jiri,
That looks okay to me.
Samtools does have an option
-x, --ignore-overlaps disable read-pair overlap detection
perhaps you need to set that.
Colin
On 21 February 2017 at 12:32, Jiri Nehyba <ji...@utexas.edu> wrote:
> Hi Colin,
>
> here are couple of alignments from sam file. I used name sorting for this
> demonstration to show reads of both orientations together. Otherwise for
> mpile I am using the same bam file sorted by plain sort command.
>
> Jiri
>
> M03964:25:000000000-AKT18:1:1101:1792:12037 83 chr16 3249344
> 42 173M = 3249342 -175 ACAACCCAGAGTTGTTGGGAAAATGAAGTA
> AGGCCCAGTGTGTCCAAGTGCCTGGCAGAGAAGAGCCCACAGGCAGGGAGTGCCTACCTT
> GTGTTCCAGGGCGACCTCCTCAATGGGGCGCACCCGGTGGCCTTGGTGCTCCTGACTCAGACTGCAGATGAGGCAGATGGGCT
> GGGFGGGGFGGGGGGGGGFGGGGGGGFGFFGGFGGFCDGGFGGGGGGGGGGGGGGGFGGG
> GGGGGGGFFGGGFFGGCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
> GGGGGGGGGGGGGGGGGGGGGGGGEGFGGGGGGGGGGGGGGGGEGGGGGFGGG AS:i:0 XN:i:0
> XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:173 YS:i:0 YT:Z:CP
> M03964:25:000000000-AKT18:1:1101:1792:12037 163 chr16 3249342
> 42 173M = 3249344 175 CAACAACCCAGAGTTGTTGGGAAAATGAAG
> TAAGGCCCAGTGTGTCCAAGTGCCTGGCAGAGAAGAGCCCACAGGCAGGGAGTGCCTACC
> TTGTGTTCCAGGGCGACCTCCTCAATGGGGCGCACCCGGTGGCCTTGGTGCTCCTGACTCAGACTGCAGATGAGGCAGATGGG
> GGGGGGGGGDGGGGGGGGGCDGGF?FFFGAFGGGGGGGGFGGGGGGGGGGGGGGG
> GGGGEGGGGGGGGGFFECFF@EG7FCGGGGFGGGDFFGFFGGGGGG?CCCEFE>
> FGGGFFFCD8CEECGGGGGGGGGGGGGGGFCFCCGGFGF>F6F>FGGFGFDFFFAD>CAFFDAF
> AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:173 YS:i:0
> YT:Z:CP
> M03964:25:000000000-AKT18:1:1101:1856:12139 83 chr16 3248843
> 42 198M = 3248841 -200 ACCCCTGCTCACTCTTCCCACCTTCCTCCC
> AGGGACGGATGGGCCATCAGCCACCTCTGACCTTACCAGAAAGCTCACTGCCTTCTCCTC
> CCCATAGGATCGCTGCTCCTCCCCTGATTTTCTCAGCTTCTTCAGATGCTCCAGCTGCTT
> CTGAATTTTCTTCTGGAAAAACAGCACTTGTTGAAAAGCTTGAATTTG FDAFFFFGGGFGGDF@5C77GGC>
> GEGEGGFEEGGGGGGGGFGGGGFDEF6E8FGGGGFGGGGGFFGGGEA9GGGGFGGGGGGG
> GGGGGGGDGGFGGGFGGGGGGGEGGGGGFGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGG
> GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGG AS:i:0 XN:i:0
> XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:198 YS:i:-3 YT:Z:CP
> M03964:25:000000000-AKT18:1:1101:1856:12139 163 chr16 3248841
> 42 198M = 3248843 200 GGACCCCTGCTCACTCTTCCCACCTTCCTC
> CCAGGGACGGATGGGCCATCAGCCACCTCTGACCTTACCAGAAAGCTCACTGCCTTCTCC
> TCCCCATAGGATCGCTGCTCCTCCCCTGATTTTCTCAGCTTCTTCAGATGCTCCAGCTGC
> TTCTGACTTTTCTTCTGGAAAAACAGCACTTGTTGAAAAGCTTGAATT
> GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGG
> GGGGGGGGGGDEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGFG
> GGGGGFGGGGEGGGCFGGGGGGGGGGGGFGGGGDGG+<D7FFGGFGFFFGFGFGF?5?FFFCFFF<@C5==@5@=CFF
> AS:i:-3 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:156A41 YS:i:0
> YT:Z:CP
> M03964:25:000000000-AKT18:1:1101:1936:12864 83 chr16 3254291
> 42 159M = 3254289 -161 AATTTCTGGATTTGCGGGCGCCTTCTCCCC
> TGTAGAAATGGTGACCTCAAGGCTTCTAGGTCGCATCTTTCCCGAGGGCAGGTACACTTC
> GAAGGGCCTGCACTCCTTCTGCCCCGGGGCGCCCCCCGCCAGCCCCTGCAGCCTCCCCGCGGAGCTGGC
> ?CGC:<FF<7:E>GGC>>DF88ECEGGEC<<+7@F<CC?CFGGGDFGGGGGGGF9DGEE>
> E:C7FGGFEGGDEE<EGFGGFGGGGGEFFGE6G?;5@CGGGGGGGGGGCEEE:GGGGGCCGGGGGCGGEDGF@7FGGGEGDGGGGGEEGCEE<GGGGDF
> AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:159 YS:i:0
> YT:Z:CP
> M03964:25:000000000-AKT18:1:1101:1936:12864 163 chr16 3254289
> 42 159M = 3254291 161 AGAATTTCTGGATTTGCGGGCGCCTTCTCC
> CCTGTAGAAATGGTGACCTCAAGGCTTCTAGGTCGCATCTTTCCCGAGGGCAGGTACACT
> TCGAAGGGCCTGCACTCCTTCTGCCCCGGGGCGCCCCCCGCCAGCCCCTGCAGCCTCCCCGCGGAGCTG
> FCFGAEFCF9CEGGFFFGGGGEEGG:<@<:@B<FGGGGGGAF,C<EGFGD@FCA@FFFFGC,9C:@C
> @:FFG,EFGGGCB@F?@EFGFGCFDEG,BCGGGCEDCF>DFG9E9@DDFGC=6+@EEGG:@:EGGG*=:8**?F4D=<49EE7DDG557C46
> AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:159 YS:i:0
> YT:Z:CP
> M03964:25:000000000-AKT18:1:1101:1992:12226 83 chr16 3243623
> 42 218M = 3243621 -220 GTCTAACACTCTTCAGATCATCAGAGAAGA
> TGAGGTTGGGGTAAGCGGTTTCTGCATCCAGAATCACATTAACTGCAAAGAAAATTTGAA
> TACCTAGGTAGGGGTCCATGGGCAACATCCCTACAGGGTTCTCCCCACCTGCAGGAAACA
> GGGACAGGGTAGTTCTTCTGGAACGTGGTAGGGGAGAGCACAGGGATCCAGCAGGCCAGGGCCACTTG
> AFFFAFB:;+AFFA6FFFGFFGFGFFFGGGFFC>GGD?C6FEGGGED@CFGGGFDCFDECE,
> 38FFEAGGGFGGGGGGGGGGECGGGGGCEFF@F9A=GGECFGGGGGDGEFGGGGGGGGGGGGE@
> CGGGGGGGGGGGGGGGGGGFGGGGFFGGGFGGFGGFGGGGGFGGFCGGGGGGGGDGGFCG
> FGGGFCFE@EFFDGFGFDCGGFF?GFFFFGFF AS:i:0 XN:i:0 XM:i:0 XO:i:0
> XG:i:0 NM:i:0 MD:Z:218 YS:i:-8 YT:Z:CP
> M03964:25:000000000-AKT18:1:1101:1992:12226 163 chr16 3243621
> 42 218M = 3243623 220 AAGTCTAACACTCTTCAGATCATCAGAGAA
> GATGAGGTTGGGGTAAGCGGTTTCTGCATCCAGAATCACATTAACTGCAAAGAAAATTTG
> AATACCGAGGTAGGGGTCCATGGGCAACATCCCTACAGGGTTATCCCCACCTGCAGGAAA
> CAGGGACAGGGTAGTTCTTCTGGAACGTGGTAGGGGAGAGGACAGGGATCCAGCAGGCCAGGGCCACT
> EEGCFDGFFEFFFD9CDEEFGGGFGGFGE,@@<AEF?C@E@FCFDGFFFE,@:EEFGFAE,B@FF
> ,?AAFA9FFEFG,EFEFFF,,AFC,A<EFG,,@>CC4C@BC=?F?9EDCGGGF=,@DFGGGEC88D,+6@=2@
> D>E6=EFGGCG61=DDD?=*=?8*8=?F7CA+3<8@?5<F?>5@A5?)*:8))0::@BAD94>0*7((/58(/885>?04
> AS:i:-8 XN:i:0 XM:i:3 XO:i:0 XG:i:0 NM:i:3 MD:Z:96T35C57C27
> YS:i:0 YT:Z:CP
> M03964:25:000000000-AKT18:1:1101:2092:11235 99 chr16 3254128
> 42 189M = 3254130 191 CGATATAAAGTAGGAAAGAACACAATTTAC
> CGGTGACCGAATTTTCTGGATTTCCAGGGCCTTCCTTCAGGTCCGCAGATGCCCCTCCAT
> CCGGCGTGGGCCTTGCCCGGGGTTCTGTTGCCGAGTCCAGATTCGCAGCTGTCTTTTCTT
> CTAGAGTCAGGAGAGTTTCTGGATTTGCGGGCGCCTTCT @+F@:FF@EFFG9F9FE<88FEF8,;
> C6CFG@E:FGGFG@@C,6CFGFFFF9FFGGGDGGGGFFDFDF9FE,
> 9EEE7++:,5,BFFGDG<EFGE+FEGG+:AFCF@,==FEF+8A,<FEG,?F7F@,3,8733EEG@
> +>EGFFCCFEC;,@7BDC,3@9C8<*,*45,2,2,2,,29:19?CDD5::2; AS:i:-12
> XN:i:0 XM:i:4 XO:i:0 XG:i:0 NM:i:4 MD:Z:42G51A53C15A24
> YS:i:-19 YT:Z:CP
> M03964:25:000000000-AKT18:1:1101:2092:11235 147 chr16 3254130
> 42 189M = 3254128 -191 ATATAAAGGAGGAAAGAACTCAATTTACCG
> GTGACCGAATGTTCTGGAGTTCCAGGGCCTTCCTTCAGGTCCGCAGAGGCCCCTCCATCC
> GGAGTGGGCCTTGCCCGGGGTTCTGTTGCCGAGGCCAGATTCGCAGCTGTCTTTTCCTCT
> AGAGTCAGGAGAATTTCTGGATTTGCGGGCGCCTGCTCC ))5A5655:8+8;1+1+4*+>:=+:81D8*
> C8@/?7GC9A70CFD;2**CFCC75;,@,EFE9C@53++++8D8>D8+ECB83,+@
> EGFGGGCF?@44+++F=+C:@+<GF8F=,+7CDFEAA,GEFC@+8:F8FC8FGFF8F@
> GCDC,9GFFF<<E<C,FFEDF6,GGGF6:::F7F7@@6CF7F@@, AS:i:-19 XN:i:0
> XM:i:6 XO:i:0 XG:i:0 NM:i:6 MD:Z:8T10A28T28T45T60T4 YS:i:-12
> YT:Z:CP
>
>
>
> On 2/20/2017 9:01 PM, Colin Hercus wrote:
>
> Hi Jiri,
>
> Could you paste a few alignments from the SAM file.
>
> Colin
>
> On 21 February 2017 at 00:31, Jiri Nehyba <ji...@utexas.edu> wrote:
>
>> Hi Tom,
>>
>> Thanks for your help.
>>
>> I think Bowtie2 is not changing base qualities. To verify that I
>> extracted from Bowtie2 bam file back the fastqs (bedtools bamtofastq)
>> and looked base scores in FastQC - they look the same as before
>> alignment. Bowtie2 is assigning alignment quality scores (bam fifth
>> column) and those are mostly perfect (151364 from 158220 bam/sam lines
>> have score 42).
>>
>> mpileup of samtools is changing the qualities. Below is one line from
>> mpileup file (samtools mpileup -Q 0 -d 1000000 -f hg38.fa PREFIX.bam >
>> PREFIX.pileup) . I folded long fifth and sixth column for better
>> readibility. This is the position where the reads have C instead of T
>> that is in reference. You might see about half of the bases have
>> qualities higher than 0 (!), most frequent quality is m that would
>> correspond to 109-33=76 on Illumina 1.8+ scale. Now original qualities
>> both in fastqs and in bam/sam file don't have zeroes and the most common
>> quality is G that is 38 on Illumina 1.8+ scale (lowest I can see just by
>> looking at random places in bam file is + at the end of reads that is
>> quality 10).
>>
>>
>> Jiri Nehyba
>>
>> chr16 3243407 T 928
>> CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
>> CCCCCCCCCCCCCCCCCCCC
>> CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
>> CCCCCCC.CCCCCCCCCCCC
>> CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
>> CCCCCCCCCCCCCCCCCCCC
>> CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
>> CCCCCCCCCCCCCCCCCCCC
>> CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
>> CCCCCCCCCCCCCCCCACCC
>> CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC.CCCCCCCCCCC
>> CCCCcccccccccccccccc
>> cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
>> cccccccccccccccccccc
>> cccccccccccccccccccccccccccccccccccccccccccccccccccc,ccccccc
>> cccccccccccccccccccc
>> cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
>> cccccccccccccccccccc
>> cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
>> cccccccccccccccccccc
>> cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
>> cacccccccccccccccccc
>> cccccccccccccccccccccccccccccccccccccccccccccccc
>> mhmmmmklgmlmmlmmmmmlmmmmmmmimlmmmljmml_mmmmm\llmmRkmmkjm_mmi
>> mmmlmmkmmlmmkmmfllmm
>> m`mmlmmmmmmm[m_mmmbkkmimmmmemmmRmmmkiml`glmlmmmmhi`_mm]
>> gmmmbmkll]emLlkmmmmmkchgi
>> jmlmlkm_ilmmmmmmmmmikmXjmhjmm\jjjm`llmmimZmmh`mmmmmllRmjkmkm
>> mkmllmmmmmmmmmmmmm[m
>> l^mmmmlmmlmhmmmmmmmmmlmmmmlmmYmmiYmmmjmmmmkmmjImmmkmmmmmgmmm
>> klmbimmlmkmmmmmlmmim
>> emmmmmlmmmlmmmmmmmmmmmmmmkimmmkhmmjlGlhmllmmmmllmm[mmmgmeimi
>> mm_bmmkmQmamjmmlmmmf
>> mmhimmmmmlmmmmimmmml_mmjmmmmmlmmlmmmmimhemhR_^mm)mmjlmmmmkmm
>> mm`m!!!!!!!!!!!!!!!!
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>> !!!!!!!!!!!!!!!!!!!!
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>> !!!!!!!!!!!!!!!!!!!!
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>> !!!!!!!!!!!!!!!!!!!!
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>> !!!!!!!!!!!!!!!!!!!!
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>> !!!!!!!!!!!!!!!!!!!!
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>>
>>
>> On 2/20/2017 6:38 AM, Thomas W. Blackwell wrote:
>> >
>> > I think this is a bowtie2 question, not a samtools question. Why is
>> > bowtie2 setting those base call qualities to zero ?
>> >
>> > - tom blackwell -
>> >
>> > On Sun, 19 Feb 2017, Jiri Nehyba wrote:
>> >
>> >> Hi,
>> >>
>> >> I would like to ask for help with mutation calling using
>> >> samtools-bcftools. Specifically my problem is that I am loosing
>> >> (almost) all my reverse bases and mutations are called only on
>> >> forward bases (reads).
>> >>
>> >> I have installed last version of samtools and bcftools (1.3.1).
>> >>
>> >> I have paired fastqs from amplicon panel libraries sequenced on
>> >> Illumina MiSeq. I cut off the primer sequences using Cutadapt,
>> >> aligned reads to hg38 genome with Bowtie2 (paired alignment), sorted
>> >> and indexed bam file with samtools and then I used following command:
>> >>
>> >> samtools mpileup -u -v -d 1000000 -f hg38.fa PREFIX.bam | bcftools
>> >> call -c -v
>> >>> PREFIX.vcf
>> >>
>> >> Resulting vcf file looks like this (showing just three rows and only
>> >> some of the fields)
>> >>
>> >> chr16 3243407 . T C 221.999 . DP=928; . . .
>> >> ;DP4=1,0,462,0;MQ=42;FQ=-281.989;PV4=1,1,0.454462,1 GT:PL
>> >> 1/1:255,255,0
>> >> chr16 3243888 . C T 221.999 . DP=2982; . . .
>> >> ;DP4=2,0,1486,2;MQ=42;FQ=-281.989;PV4=1,0.458362,0.4298,1 GT:PL
>> >> 1/1:255,255,0
>> >> chr16 3243922 . A T 221.999 . DP=2982; . . .
>> >> ;DP4=1,0,1486,4;MQ=42;FQ=-281.989;PV4=1,0.338485,0.450054,1 GT:PL
>> >> 1/1:255,255,0
>> >>
>> >> What I don't like is the DP4 field that suggests extreme strand bias.
>> >>
>> >> I tried to find out why are those reverse bases discarded - there is
>> >> no reason for generally lower quality score - amplicons are ligated
>> >> randomly in both orientations and I am getting both sides of each
>> >> amplicon in R1 and R2 reads.
>> >>
>> >> The only thing that prevents the "loss" of the bases is setting the Q
>> >> to 0, like that:
>> >>
>> >> samtools mpileup -u -v -Q 0 -d 1000000 -f hg38.fa PREFIX.bam |
>> >> bcftools call -c -v > PREFIX.vcf
>> >>
>> >> Result:
>> >>
>> >> chr16 3243407 . T C 221.999 . DP=928; . . .
>> >> ;DP4=2,1,462,463;MQ=42;FQ=-281.989;PV4=1,1,0.419714,1 GT:PL
>> >> 1/1:255,255,0
>> >> chr16 3243888 . C T 221.999 . DP=2982; . . .
>> >> ;DP4=2,2,1489,1489;MQ=42;FQ=-281.989;PV4=1,0.493227,0.40081,1 GT:PL
>> >> 1/1:255,255,0
>> >> chr16 3243922 . A T 221.999 . DP=2982; . . .
>> >> ;DP4=4,3,1487,1488;MQ=42;FQ=-281.989;PV4=1,1,1,1 GT:PL
>> 1/1:255,255,0
>> >>
>> >> If I set Q to more than 0, for example 1 I will again loose almost
>> >> all reverse bases.
>> >>
>> >> Finally, a look at mpileup file obtained by this command: samtools
>> >> mpileup -Q 0 -d 1000000 -f hg38.fa PREFIX.bam > PREFIX.pileup
>> >>
>> >> I see that half of the bases is getting score ! that is 0 :
>> >>
>> >>
>> >>
>> >> Thank you for your help,
>> >>
>> >> Jiri
>> >>
>>
>>
>> ------------------------------------------------------------
>> ------------------
>> Check out the vibrant tech community on one of the world's most
>> engaging tech sites, SlashDot.org! <http://sdm.link/slashdot>
>> http://sdm.link/slashdot
>> _______________________________________________
>> Samtools-help mailing list
>> Samtools-help@lists.sourceforge.net
>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>>
>
>
>
------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, SlashDot.org! http://sdm.link/slashdot
_______________________________________________
Samtools-help mailing list
Samtools-help@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/samtools-help
