SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Question on calling SNPs using samtools/bcftools nkwuji Bioinformatics 6 02-19-2013 08:52 AM
From sam file, how to the mate aligned or not? fabrice Bioinformatics 2 08-09-2011 10:37 AM
SAMTOOLS SNP calling question harrb Bioinformatics 2 12-10-2010 06:37 AM
question on SAMtools consensus calling orionzhou Bioinformatics 9 11-16-2010 02:42 PM
Quality values for aligned reads in .sam file arnkas Bioinformatics 2 11-08-2009 03:41 PM

Reply
 
Thread Tools
Old 02-01-2010, 05:56 PM   #1
Ling
Junior Member
 
Location: Berkeley, CA

Join Date: Jan 2010
Posts: 8
Default Question on SAMtools variations calling of a BFAST-aligned SAM file?

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
Ling is offline   Reply With Quote
Old 02-01-2010, 07:59 PM   #2
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Quote:
Originally Posted by Ling View Post
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
nilshomer is offline   Reply With Quote
Old 02-02-2010, 11:17 AM   #3
Ling
Junior Member
 
Location: Berkeley, CA

Join Date: Jan 2010
Posts: 8
Default

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 is offline   Reply With Quote
Old 02-02-2010, 01:58 PM   #4
Ling
Junior Member
 
Location: Berkeley, CA

Join Date: Jan 2010
Posts: 8
Default

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
Ling is offline   Reply With Quote
Old 02-02-2010, 06:26 PM   #5
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Quote:
Originally Posted by Ling View Post
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.
nilshomer is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 01:31 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO