A BLAST to SAM converter.
Some times ago, I've received a set of Ion-Torrent /mate-reads with a poor quality. I wasn't able to align much things using bwa. I've always wondered if I could get better alignments using NCBI-BLASTN (short answer: no) . That's why I asked guyduche, my intership student to write a C program to convert the output of blastn to SAM. His code is available on github at :
The input for blast2sam is- the XML output of NCBI blastn (or stdin)
- The single or pair of fastq file(s)
- The reference sequence indexed with picard
Example:
fastq2fasta in.R1.fq.gz in.R2.fq.gz |\ blastn -db REFERENCE -outfmt 5 | \ blast2bam -o result.bam -W 40 -R '@RG ID:foo SM:sample' - REFERENCE.dict in.R1.fq.gz in.R2.fq.gz
Output:
@SQ SN:gi|9629357|ref|NC_001802.1| LN:9181 @RG ID:foo SM:sample @PG ID:Blast2Bam PN:Blast2Bam VN:0.1 CL:../../bin/blast2bam -o results.sam -W 40 -R @RG ID:foo SM:sample - db.dict test_1.fastq.gz test_2.fastq.gz (...) ERR656485.2 83 gi|9629357|ref|NC_001802.1| 715 60 180S7=1X8=1X11=1X2=2X4=1X14=1X8=1X33=1X4=1X2=1X5=1X2=1X6=1S = 715 -119 CCTAGTGTTGCTTGCTTTTCTTCTTTTTTTTTTCAAGCAGAAGACGGCATACGAGATCCTCTATCGAGATCGGTCTCGGCATTCCTGCTGAACCGCTCTTCCGATCTAAGATAGAGGAAGAACAAAACAAATGTCAGCAAAGTCAGCAAAAGACACAGCAGGAAAAAGGGGCTGACGGGAAGGTCAGTCAAAATTATCCTATAGTGCAAAATCTCCAAGGGCAAATGGTACACCAGGCCATGTCACCTAGAACTTTAAATGCATGGGTAAAAGTAATAGAGGAAAAGGCCTTTAGCCCAN (),.((((,(((((,((.((.-(>69>20E>6/=>5EC@9-52?BEE::2951.)74B64=B==FFAF=A??59:>FFFDF:55GGFGF?DFGGFE868>GGGFGGGGED;FGFFGGGGGGGGGGGEFFGE9GGGGFGGGGGGGGDGECGGFGGGGGGGGGGFGGGGGEGGFGGGGGGFFGGGGGFF?EGGFFFEGGGGGGGGFEGGGEGGGFEGGGGGGGGGGDGFFCEGFGGGGGGGGGGGFFECFGGGGFGGGGGGGGGGGFCGGGGGGGGGGGGGGGGGGFGGGGGGGGF@CCA8! NM:i:13 RG:Z:foo AS:i:80 XB:f:148.852 XE:Z:4.07e-39 ERR656485.2 163 gi|9629357|ref|NC_001802.1| 715 60 73S7=1X8=1X11=1X2=2X4=1X14=1X8=1X33=1X4=1X2=1X5=1X2=1X8=106S = 715 119 NAGATAGAGGAAGAACAAAACAAATGTCAGCAAAGTCAGCAAAAGACACAGCAGGAAAAAGGGGCTGACGGGAAGGTCAGTCAAAATTATCCTATAGTGCAAAATCTCCAAGGGCAAATGGTACACCAGGCCATGTCACCTAGAACTTTAAATGCATGGGTAAAAGTAATAGAGGAAAAGGCCTTTAGCCCAGAGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTCTCATTACAAAAAAAACATACACAATAAATGATATAAGCGGAATCAACAGCATGA !8A@CGGEFGFGCDFGGGGGGGGGGGGGGFGGGGGFGFGGGGGGGGGGGGGGGGGGGGGGGGEGGGGGGGGGGGGGGFGGFGGGGGGGGGEGFGFGGGFFGGGGGGGGFGGGGGGGGGGGGFFFFGGGGGG=FFGGFFDGGGGGGGG8FGFGGGGGGGGGFGGGGGGGGGGFDGGFGGFGGGFFFGFF8DFDFDFFFFFFFFFBCDB<@EAFB@ABAC@CDFF?4>EEFE<*>BDAFB@FFBFF>((6<5CC.;C;=D9106(.))).)-46<<))))))))))((,(-)))()((())) NM:i:13 RG:Z:foo AS:i:82 XB:f:152.546 XE:Z:3.15e-40 (...)Now, I would be interested in finding another dataset where this tool could be successfully used.
That's it,
Pierre
8 comments:
Could it be used to align a de novo genome assembly and then get pseudochromosomes from it, rather than map reads and get consensus?
I suppose that such tool blast->assembly already exists (?)
this is not a XML document. You did it wrong, check your "output-format" blast parameter.
See also: https://github.com/guyduche/Blast2Bam
Post a Comment