Question on SAMtools variations calling of a BFAST-aligned SAM file?

[Ling]

Hi all,


I had problems in using SAMtools to call the variations from a BFAST- aligned SAM file.

The commands I used:

$ samtools-0.1.7/samtools view -bS -q 20 bfast.aln.sam > bfast.aln.bam
$ samtools-0.1.7/samtools sort bfast.aln.bam bfast-aln.sorted
$ samtools-0.1.7/samtools index bfast-aln.sorted.bam
$ samtools-0.1.7/samtools pileup -vcf TAIR9genomeSeq.txt bfast-aln.sorted.bam | tee bfast-raw.txt | samtools-0.1.7/samtools.pl varFilter -D100 > bfast-flt.txt

It had run for days in pileup step, with no error information and no output.

The SAMtools worked fine with the BWA-aligned sam file from the same data.

Is anything wrong with my Bfast-sam file? which look like:
@HD VN:0.1.2
@SQ SN:Chr1 LN:30427671
@SQ SN:Chr2 LN:19698289
@SQ SN:Chr3 LN:23459830
@SQ SN:Chr4 LN:18585056
@SQ SN:Chr5 LN:26975502
@SQ SN:chloroplast LN:154478
@SQ SN:mitochondria LN:366924
@PG ID:bfast VN:0.6.2a
HWI-EAS412:2:1:0:1469#0/1 4 * 0 0 * * nttcctatcaaggttaaggagtcaactgagcatctcagcggcgggattgaatacccggatcgaatcagagttcacg DOXVYWYWWWWVWWYYWWTWVWWWWWWWWYWYWWYWUSVYVWWWVNUUSUWWWQUGVRRRUPTUPWUUWTVSLRPU PG:Z:bfast AS:i:-2147483648 NM:i:0 NH:i:1 IH:i:1 HI:i:1 MD:Z:* XA:i:3
HWI-EAS412:2:1:1:1078#0/1 4 * 0 0 * * ncgccccggccgatatcctcaaccagctggccacgtagatcggaagagctcgtatgccgtcttcttctttaaaaaa DNVUWXXVWWSUVWWXUUWWVWWUWWVNSUURWQVWUSXWUUUUWWSSNUNSVWWFNJNVVVVQTKOVSGOVSUVS PG:Z:bfast AS:i:-2147483648 NM:i:0 NH:i:1 IH:i:1 HI:i:1 MD:Z:* XA:i:3
HWI-EAS412:2:1:1:660#0/1 16 Chr3 13591936 10 76M tgattaagtataagaacttaaaccgcaacatgatcttaaaggcgtaagaattgtatccttgttaaaagacacaaan BBXXXZ[[Z[Y[WWW[[[[[YVRVY[[[[YZY[ZYULV[[ZYY[[[ZZZ[[[[[[XVSX[[[[[[[[ZYYYVY[OD PG:Z:bfast AS:i:3000 NM:i:4 NH:i:1 IH:i:1 HI:i:1 MD:Z:5G23CC44G XA:i:3
HWI-EAS412:2:1:1:1493#0/1 0 Chr1 14102811 40 76M natcttcttgaagtgtgaaatacttcatcaaacggtggagaatacttctcccttggtcactcgatcaaacggtgag DOZZZZZZYYYYYXSXYYZZZZXZZZZZXZZZZXZZYYYXZZZZZZZZZYVXZZYWYWZZZXZZXXZZZZUURZXV PG:Z:bfast AS:i:3600 NM:i:1 NH:i:1 IH:i:1 HI:i:1 MD:Z:T75 XA:i:3
HWI-EAS412:2:1:1:1133#0/1 16 Chr4 3403799 50 76M * gtttgttctctcacttggcgagcaatacggtcgatgttatcgttgaatagaaggttctgaattccttgtgaccgtn BTWYYWYVYWWWTVWVYYWWYSSVYYVSTVWYWZZZXZZYWYZXZZZZZZZYYYZWSWYYZZWUUWZZYYVPLTOD PG:Z:bfast AS:i:3600 NM:i:1 NH:i:1 IH:i:1 HI:i:1 MD:Z:75G XA:i:3
HWI-EAS412:2:1:1:144#0/1 16 Chr2 15897161 255 76M agtaccgaatcgttggatgatgatggtagtgatgatcgtggtggtaatgagatccgctagaggaggaaccatttgn YYXWUWVWYWVWYVVUVYUYXWYXYXYYYYYYYYYXYYYYXYYYXYYYYYYYYYYYXYYXYYYYYYYYYYURUYOD PG:Z:bfast AS:i:3400 NM:i:2 NH:i:1 IH:i:1 HI:i:1 MD:Z:71A3A XA:i:3
HWI-EAS412:2:1:1:502#0/1 4 * 0 0 * * ntcttagatctcacaaaaaagaccaaaagtaattcaactgatgaatttagataaagtaatttgtcttctaaaccca DP[[[ZYZ[[[[[[[[[[[ZYZ[[YYYYXY[[[[[Z[[[[[[[[Z[[[YTYZ[[[[Y[Z[[XXXY[[W[YYXURVW PG:Z:bfast AS:i:-2147483648 NM:i:0 NH:i:1 IH:i:1 HI:i:1 MD:Z:* XA:i:3
HWI-EAS412:2:1:1:1922#0/1 0 mitochondria 210296 255 76M naagagcaaagcgtggctttttgccccggggcactagggcctgatcatctttatggctataaatctggcttccaag DOWWXWXYWYYYYWYYWYYYYWTUWXWXWWWWWWWWWXWWVWWYWWWWWWYWWWBBBBBBBBBBBBBBBBBBBBBB PG:Z:bfast AS:i:3000 NM:i:4 NH:i:1 IH:i:1 HI:i:1 MD:Z:C66C4AC2 XA:i:3
HWI-EAS412:2:1:1:1659#0/1 16 Chr2 18335755 255 76M ggcaggaacaaacagaaacagacagagtcgaaaacaaattcgtcgaatgataatagtagtaaacaagacacaggtn BTOTQOSYVYYXYXXXYYYXXXZZZZZZZXZZZZZZZZZZZZYXYZZZZZZZZZZZZZZXZZYYYZZZXZZZZZND PG:Z:bfast AS:i:3600 NM:i:1 NH:i:1 IH:i:1 HI:i:1 MD:Z:75A XA:i:3
HWI-EAS412:2:1:1:477#0/1 0 Chr3 18560494 255 76M natgtggttcagttttgtatggtttttttgtctgtgactatttgttaatggtacgtttttgaatgccccaaatcaa DM[[WSSLPWY[W[Y[[WZ[W[ZYVYYYY[[VZYSVYXZX[XXTYYVYXPXVSRT[YYYXXRUXXSYLQQUXXSUX PG:Z:bfast AS:i:3600 NM:i:1 NH:i:1 IH:i:1 HI:i:1 MD:Z:A75 XA:i:3
HWI-EAS412:2:1:1:1430#0/1 0 Chr4 4184778 28 76M * ntcgagtagtttcgtttcctcagtttacaagtgattgacatgctcaataaattaccctacacaacaagttcattaa DOVUXXXVXXXTTRVXYVTVYXXXZZZVTSYXXXYYYTTTXVVVXTXXYVXXYYVRTVVUWWSSOQVTXSOQYYTT PG:Z:bfast AS:i:3600 NM:i:1 NH:i:1 IH:i:1 HI:i:1 MD:Z:A75 XA:i:3


My reads are 76 bp and I used bfast-0.6.2a

Thank you very much for your attention and help!


Ling

[nilshomer]

Hi all,


I had problems in using SAMtools to call the variations from a BFAST- aligned SAM file.

The commands I used:

$ samtools-0.1.7/samtools view -bS -q 20 bfast.aln.sam > bfast.aln.bam
$ samtools-0.1.7/samtools sort bfast.aln.bam bfast-aln.sorted
$ samtools-0.1.7/samtools index bfast-aln.sorted.bam
$ samtools-0.1.7/samtools pileup -vcf TAIR9genomeSeq.txt bfast-aln.sorted.bam | tee bfast-raw.txt | samtools-0.1.7/samtools.pl varFilter -D100 > bfast-flt.txt

It had run for days in pileup step, with no error information and no output.

The SAMtools worked fine with the BWA-aligned sam file from the same data.

Is anything wrong with my Bfast-sam file? which look like:
@HD VN:0.1.2
@SQ SN:Chr1 LN:30427671
@SQ SN:Chr2 LN:19698289
@SQ SN:Chr3 LN:23459830
@SQ SN:Chr4 LN:18585056
@SQ SN:Chr5 LN:26975502
@SQ SN:chloroplast LN:154478
@SQ SN:mitochondria LN:366924
@PG ID:bfast VN:0.6.2a
HWI-EAS412:2:1:0:1469#0/1 4 * 0 0 * * nttcctatcaaggttaaggagtcaactgagcatctcagcggcgggattgaatacccggatcgaatcagagttcacg DOXVYWYWWWWVWWYYWWTWVWWWWWWWWYWYWWYWUSVYVWWWVNUUSUWWWQUGVRRRUPTUPWUUWTVSLRPU PG:Z:bfast AS:i:-2147483648 NM:i:0 NH:i:1 IH:i:1 HI:i:1 MD:Z:* XA:i:3
HWI-EAS412:2:1:1:1078#0/1 4 * 0 0 * * ncgccccggccgatatcctcaaccagctggccacgtagatcggaagagctcgtatgccgtcttcttctttaaaaaa DNVUWXXVWWSUVWWXUUWWVWWUWWVNSUURWQVWUSXWUUUUWWSSNUNSVWWFNJNVVVVQTKOVSGOVSUVS PG:Z:bfast AS:i:-2147483648 NM:i:0 NH:i:1 IH:i:1 HI:i:1 MD:Z:* XA:i:3
HWI-EAS412:2:1:1:660#0/1 16 Chr3 13591936 10 76M tgattaagtataagaacttaaaccgcaacatgatcttaaaggcgtaagaattgtatccttgttaaaagacacaaan BBXXXZ[[Z[Y[WWW[[[[[YVRVY[[[[YZY[ZYULV[[ZYY[[[ZZZ[[[[[[XVSX[[[[[[[[ZYYYVY[OD PG:Z:bfast AS:i:3000 NM:i:4 NH:i:1 IH:i:1 HI:i:1 MD:Z:5G23CC44G XA:i:3
HWI-EAS412:2:1:1:1493#0/1 0 Chr1 14102811 40 76M natcttcttgaagtgtgaaatacttcatcaaacggtggagaatacttctcccttggtcactcgatcaaacggtgag DOZZZZZZYYYYYXSXYYZZZZXZZZZZXZZZZXZZYYYXZZZZZZZZZYVXZZYWYWZZZXZZXXZZZZUURZXV PG:Z:bfast AS:i:3600 NM:i:1 NH:i:1 IH:i:1 HI:i:1 MD:Z:T75 XA:i:3
HWI-EAS412:2:1:1:1133#0/1 16 Chr4 3403799 50 76M * gtttgttctctcacttggcgagcaatacggtcgatgttatcgttgaatagaaggttctgaattccttgtgaccgtn BTWYYWYVYWWWTVWVYYWWYSSVYYVSTVWYWZZZXZZYWYZXZZZZZZZYYYZWSWYYZZWUUWZZYYVPLTOD PG:Z:bfast AS:i:3600 NM:i:1 NH:i:1 IH:i:1 HI:i:1 MD:Z:75G XA:i:3
HWI-EAS412:2:1:1:144#0/1 16 Chr2 15897161 255 76M agtaccgaatcgttggatgatgatggtagtgatgatcgtggtggtaatgagatccgctagaggaggaaccatttgn YYXWUWVWYWVWYVVUVYUYXWYXYXYYYYYYYYYXYYYYXYYYXYYYYYYYYYYYXYYXYYYYYYYYYYURUYOD PG:Z:bfast AS:i:3400 NM:i:2 NH:i:1 IH:i:1 HI:i:1 MD:Z:71A3A XA:i:3
HWI-EAS412:2:1:1:502#0/1 4 * 0 0 * * ntcttagatctcacaaaaaagaccaaaagtaattcaactgatgaatttagataaagtaatttgtcttctaaaccca DP[[[ZYZ[[[[[[[[[[[ZYZ[[YYYYXY[[[[[Z[[[[[[[[Z[[[YTYZ[[[[Y[Z[[XXXY[[W[YYXURVW PG:Z:bfast AS:i:-2147483648 NM:i:0 NH:i:1 IH:i:1 HI:i:1 MD:Z:* XA:i:3
HWI-EAS412:2:1:1:1922#0/1 0 mitochondria 210296 255 76M naagagcaaagcgtggctttttgccccggggcactagggcctgatcatctttatggctataaatctggcttccaag DOWWXWXYWYYYYWYYWYYYYWTUWXWXWWWWWWWWWXWWVWWYWWWWWWYWWWBBBBBBBBBBBBBBBBBBBBBB PG:Z:bfast AS:i:3000 NM:i:4 NH:i:1 IH:i:1 HI:i:1 MD:Z:C66C4AC2 XA:i:3
HWI-EAS412:2:1:1:1659#0/1 16 Chr2 18335755 255 76M ggcaggaacaaacagaaacagacagagtcgaaaacaaattcgtcgaatgataatagtagtaaacaagacacaggtn BTOTQOSYVYYXYXXXYYYXXXZZZZZZZXZZZZZZZZZZZZYXYZZZZZZZZZZZZZZXZZYYYZZZXZZZZZND PG:Z:bfast AS:i:3600 NM:i:1 NH:i:1 IH:i:1 HI:i:1 MD:Z:75A XA:i:3
HWI-EAS412:2:1:1:477#0/1 0 Chr3 18560494 255 76M natgtggttcagttttgtatggtttttttgtctgtgactatttgttaatggtacgtttttgaatgccccaaatcaa DM[[WSSLPWY[W[Y[[WZ[W[ZYVYYYY[[VZYSVYXZX[XXTYYVYXPXVSRT[YYYXXRUXXSYLQQUXXSUX PG:Z:bfast AS:i:3600 NM:i:1 NH:i:1 IH:i:1 HI:i:1 MD:Z:A75 XA:i:3
HWI-EAS412:2:1:1:1430#0/1 0 Chr4 4184778 28 76M * ntcgagtagtttcgtttcctcagtttacaagtgattgacatgctcaataaattaccctacacaacaagttcattaa DOVUXXXVXXXTTRVXYVTVYXXXZZZVTSYXXXYYYTTTXVVVXTXXYVXXYYVRTVVUWWSSOQVTXSOQYYTT PG:Z:bfast AS:i:3600 NM:i:1 NH:i:1 IH:i:1 HI:i:1 MD:Z:A75 XA:i:3


My reads are 76 bp and I used bfast-0.6.2a

Thank you very much for your attention and help!


Ling

What does the "samtools-0.1.7/samtools pileup -vcf TAIR9genomeSeq.txt" command give you? The other steps after pileup I have not seen before.

Nils

[Ling]

Hi Nils,

Thank you very much for your reply.

I always pipe the pileup results into filter step and thus never got the output from pileup.

I re-ran the pileup step after read your message. Here is the results:
Chr1 31 T T 33 0 47 2 .+3AAA. `b
Chr1 31 * */+AAA 32 32 47 2 * +AAA 0
Chr1 144 T W 0 0 60 2 ,a aB
Chr1 145 A A 33 0 60 2 ,,-1t aB
Chr1 145 * -T/-T 0 0 60 2 -T * 0
Chr1 1075 T C 30 30 60 1 c+3gag B
Chr1 1075 * +GAG/+GAG 41 35 60 1 +GAG 0
Chr1 1096 C C 30 0 60 1 ,+8cccgccct ]
Chr1 1096 * +CCCGCCCT/+CCCGCCCT 41 35 60 1 +CCCGCCCT * 1 0 0 0 0
Chr1 1099 A A 30 0 60 1 ,+1c M
Chr1 1099 * +C/+C 41 35 60 1 +C * 0
Chr1 1100 A A 30 0 60 1 ,+9gccgcccgt a
Chr1 1100 * +GCCGCCCGT/+GCCGCCCGT 41 35 60 1 +GCCGCCCGT * 1 0 0 0 0
Chr1 1102 T T 30 0 60 1 ,+5ccccc [
Chr1 1102 * +CCCCC/+CCCCC 41 35 60 1 +CCCCC 0
Chr1 1104 A C 30 30 60 1 c-1g X
Chr1 1104 * -G/-G 41 35 60 1 -G * 0
Chr1 1105 G A 0 0 0 1 * a
Chr1 1107 A A 30 0 60 1 ,+8ctttctct _
Chr1 1107 * +CTTTCTCT/+CTTTCTCT 41 35 60 1 +CTTTCTCT * 1 0 0 0 0
Chr1 1109 T T 30 0 60 1 ,+2gc a
Chr1 1109 * +GC/+GC 41 35 60 1 +GC * 0
Chr1 1579 T W 2 2 60 2 ,n [D
Chr1 1596 T Y 21 21 60 2 c, Wa
Chr1 2467 G R 19 19 60 5 ,.a., `_^a_
Chr1 3272 A W 28 28 60 3 t,, abb
Chr1 4018 A R 19 19 60 5 ,.g.. a`^``
Chr1 4030 G S 1 1 60 7 ,.,C... `]`Raa_
Chr1 4134 A W 7 7 60 10 ..,,,,,,T, ^[_a^`aaa]
Chr1 4264 A W 29 29 60 2 ,t a_
Chr1 4382 G R 2 2 60 2 N. Da
Chr1 4423 A R 22 29 60 2 ,G X_
Chr1 4635 T W 2 2 60 11 ,..A,,..... baa_aababab
Chr1 4654 A R 8 8 60 9 ..,,.G... __bbb_a``
Chr1 4870 C M 14 14 60 7 .,.,,a. [bBa`_a
Chr1 4881 T T 42 0 60 5 .+1G,,., BaYa_
Chr1 4881 * */+G 23 23 60 5 * +G 0
Chr1 4885 T T 42 0 60 5 .+1G,,., Ba`U`
Chr1 4885 * */+G 23 23 60 5 * +G 0


It looks good for me, and now I am runing the filtering step. I will let you know the results

Best,

Ling

[Ling]

Hi Nils,

It's working now!

But I do still not understand why it wasn't working by using:

$ samtools-0.1.7/samtools pileup -vcf TAIR9genomeSeq.txt bfast-aln.sorted.bam | tee bfast-raw.txt | samtools-0.1.7/samtools.pl varFilter -D100 > bfast-flt.txt

but working by:
$ samtools-0.1.7/samtools pileup -vcf TAIR9genomeSeq.txt bfast-aln.sorted.bam > bfast.pileup;

$ samtools-0.1.7/samtools.pl varFilter -D100 bfast.pileup > bwa-flt.txt

Thank you for you help!

Ling

[nilshomer]

Hi Nils,

It's working now!

But I do still not understand why it wasn't working by using:

$ samtools-0.1.7/samtools pileup -vcf TAIR9genomeSeq.txt bfast-aln.sorted.bam | tee bfast-raw.txt | samtools-0.1.7/samtools.pl varFilter -D100 > bfast-flt.txt

but working by:
$ samtools-0.1.7/samtools pileup -vcf TAIR9genomeSeq.txt bfast-aln.sorted.bam > bfast.pileup;

$ samtools-0.1.7/samtools.pl varFilter -D100 bfast.pileup > bwa-flt.txt

Thank you for you help!

Ling

Very odd indeed. I am glad it is working for you.