I have two sequences like those at the end. I want to find all local alignment larger than 10 and maybe shorter than some number (maybe I can filter that out later).
I was trying to do it with a smith-waterman, but I cannot find an implementation that retrieves me the from / to value to indentify where the results come from.
I want to allow mismatch and gaps, any ideas on what tool can I use?
It has to be fast since I will be splitting a full chromosome / genome and re doing this task many times. I'd also like to avoid reverse complement pairing.
>1
ATGGCTAGGAAACATAACCATCTTTGATCAACGAGCTAGTCAAGTAGAGGCATACTAGTG
ACACTCTGTTTGTCTATGTATTCACACATGTATCATGTTTCCGGTTAATACAATTCTAGC
ATGAATAATAAACATTTATCATGATATAAGGAAATAAATAATAACTTTATTGTTGCCTCT
AGGGCATATTTCCTTCAGAGTGCTAGTTGAAGGCACGTGATGTGCATCTTAGGTTAAACT
GTAAATCGAATCTTACTCTCGAGTAACTCCAAAAGCGGCCACCGAAAATCATTATTGGAT
TTGTTTTTCATTACCAGCACGTCGAGGTGAAAGATGTGCAGGTTCAGTGGGTCGACTTGC
GTTTTCACTGACATGTGGGACCGAATGTAAGCAAACAGACGATGGGCACACTCCGCGCGT
CTGCTGGCTGGCTGAGAAGAGAGAGGAGGGTTGAGCACGGACTGTAGGAAATTTCGTGAC
AATGTGAGTAGTACTAGTGGTGTACTGTACTCCTCGAGAGATAATGTGTTCGCTCTACGT
AACCAGCCCCTCGATATAGTGGGGTGGCATGGGCCATAACTAGCATTCACGTCTAACAGT
ACGGCACATTCCACCAATTTTCTTGGAATCCATGCTACAGCTATCTCTTCTCCCTCCTAT
TTTCTCTAACTCAGCCATCGCGCGGTGGGAGGGGTGGGGGTTTGCCGATGTCAGTTTGCT
CAACTTCGGTCCCACTTGTAAGTGAAAATGCAACACAGCACGCATCCCCCCTTGAGCGGC
GTGCCAGACCAACCGTGTGGAGGTGCGATGCGAAAATGGCAGGCTTTTCCTTTTGAACTT
GTAACATGTAAATTTTGGTTCTGCTCCTTGGCCCGACCGATAGATAGTAATGAGAAGTTG
ATGCCAAAAAATGAGAAGTTTGGTTCTGGCGCGGTTCATGCATGGTACGTACAAAAATTA
TGAAGTTGATGCAAAAGGCAAGATAATTACTAGTACAATATAGTCCCTCCATTCCTAAAT
ATTTGTATTTTTAGGAAAGATACATATGGATCCACCAACGACATCCTCACGCCCTAGTCC
GCCCGGTCACCAACCAACGGGTGAGGAGCAGCAGCGTTGGCGGCGAGCTCAACGTGCTTC
>1rc
CAATCGCAAGATCAGTTGACTAATGGAAGCTATACCCATGGAGGCGATGAGCGCGAGCAT
CAGCCGACTGTGCCAGATGGTCCATGACGCCGGCCTGCGGCCTGGCACCGAGGAACGTCT
CCAGGTCGTGCTTGAGGCCGCTAGGGCGAGAGGTCTCTTGGACGACAGCTTCGTCTCCTT
GTTCGACAAGGTCCTCGTTGGATTCCTCGACAAGTTCAACGTCGTGAAGAAGCTCACGGA
CGACCTTGACATACGCCTCCAGCCCACGCGCCCAGGCTCTACGATGCCCGCCACCCTCAA
CGACCTCTATGATGACAACCTCTTCGATGCACTGGTGGACCTGCGACTGCCCGTCGTCGT
GCCGGAGATTGTCCACCTCGAGGTCACGCTCGCCGCGCAGCGCCTGGCACAGTAAGACAC
CATCGACATAATCACCCACGTCTACGCGCAAATCGTCCACAAGGACTACTACATGCCAGA
GGAGGAGGATAGGACGCTGGCCTTCTTGGAACGCAGGGCAACCTTGGACGGCATTGTTTA
GAAGCACGTTGAGCTCGCCGCCAACGCTGCTGCTCCTCACCCGTTGGTTGGTGACCGGGC
GGACTAGGGCGTGAGGATGTCGTTGGTGGATCCATATGTATCTTTCCTAAAAATACAAAT
ATTTAGGAATGGAGGGACTATATTGTACTAGTAATTATCTTGCCTTTTGCATCAACTTCA
TAATTTTTGTACGTACCATGCATGAACCGCGCCAGAACCAAACTTCTCATTTTTTGGCAT
CAACTTCTCATTACTATCTATCGGTCGGGCCAAGGAGCAGAACCAAAATTTACATGTTAC
AAGTTCAAAAGGAAAAGCCTGCCATTTTCGCATCGCACCTCCACACGGTTGGTCTGGCAC
•
written
by
juan.crescente •
BLAST will give you all possible local alignments that satisify certain alignment parameters. Therefore, make sure you:
To run BLAST to align only two sequences (rather than search through database) just type:
blastn -query seq1.fa -subject seq2.fa -dust no -soft_masking false -evalue 1e+14 -word_size 4 -task blastn-short
You can add -outfmt 7 to save results in tab-separeted format.
You can also do this comparison on-line on the NCBI BLAST website by checking the option Align two or more sequences.
Sample output:
BLASTN 2.2.28+
Query= 1
Length=1140
Subject= 1rc
Length=900
Score = 714 bits (360), Expect = 0.0
Identities = 360/360 (100%), Gaps = 0/360 (0%)
Strand=Plus/Minus
Query 781 GTGCCAGACCAACCGTGTGGAGGTGCGATGCGAAAATGGCAGGCTTTTCCTTTTGAACTT 840
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 900 GTGCCAGACCAACCGTGTGGAGGTGCGATGCGAAAATGGCAGGCTTTTCCTTTTGAACTT 841
Query 841 GTAACATGTAAATTTTGGTTCTGCTCCTTGGCCCGACCGATAGATAGTAATGAGAAGTTG 900
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 840 GTAACATGTAAATTTTGGTTCTGCTCCTTGGCCCGACCGATAGATAGTAATGAGAAGTTG 781
Query 901 ATGCCAAAAAATGAGAAGTTTGGTTCTGGCGCGGTTCATGCATGGTACGTACAAAAATTA 960
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 780 ATGCCAAAAAATGAGAAGTTTGGTTCTGGCGCGGTTCATGCATGGTACGTACAAAAATTA 721
Query 961 TGAAGTTGATGCAAAAGGCAAGATAATTACTAGTACAATATAGTCCCTCCATTCCTAAAT 1020
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 720 TGAAGTTGATGCAAAAGGCAAGATAATTACTAGTACAATATAGTCCCTCCATTCCTAAAT 661
Query 1021 ATTTGTATTTTTAGGAAAGATACATATGGATCCACCAACGACATCCTCACGCCCTAGTCC 1080
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 660 ATTTGTATTTTTAGGAAAGATACATATGGATCCACCAACGACATCCTCACGCCCTAGTCC 601
Query 1081 GCCCGGTCACCAACCAACGGGTGAGGAGCAGCAGCGTTGGCGGCGAGCTCAACGTGCTTC 1140
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 600 GCCCGGTCACCAACCAACGGGTGAGGAGCAGCAGCGTTGGCGGCGAGCTCAACGTGCTTC 541
Score = 22.3 bits (11), Expect = 0.19
Identities = 11/11 (100%), Gaps = 0/11 (0%)
Strand=Plus/Minus
Query 894 GAAGTTGATGC 904
|||||||||||
Sbjct 719 GAAGTTGATGC 709
Score = 22.3 bits (11), Expect = 0.19
Identities = 11/11 (100%), Gaps = 0/11 (0%)
Strand=Plus/Minus
Query 962 GAAGTTGATGC 972
|||||||||||
Sbjct 787 GAAGTTGATGC 777
Score = 20.3 bits (10), Expect = 0.77
Identities = 10/10 (100%), Gaps = 0/10 (0%)
Strand=Plus/Minus
Query 833 TTGAACTTGT 842
||||||||||
Sbjct 219 TTGAACTTGT 210
Score = 20.3 bits (10), Expect = 0.77
Identities = 10/10 (100%), Gaps = 0/10 (0%)
Strand=Plus/Minus
Query 854 TTTGGTTCTG 863
||||||||||
Sbjct 762 TTTGGTTCTG 753
Score = 20.3 bits (10), Expect = 0.77
Identities = 10/10 (100%), Gaps = 0/10 (0%)
Strand=Plus/Minus
Query 919 TTTGGTTCTG 928
||||||||||
Sbjct 827 TTTGGTTCTG 818
and the list continues....
Sample tab-separeted output:
# BLASTN 2.2.28+
# Query: 1
# Subject: 1rc
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 5677 hits found
1 1rc 100.00 360 0 0 781 1140 900 541 0.0 714
1 1rc 100.00 11 0 0 894 904 719 709 0.19 22.3
1 1rc 100.00 11 0 0 962 972 787 777 0.19 22.3
1 1rc 100.00 10 0 0 833 842 219 210 0.77 20.3
1 1rc 100.00 10 0 0 854 863 762 753 0.77 20.3
1 1rc 100.00 10 0 0 919 928 827 818 0.77 20.3
1 1rc 92.31 13 1 0 1065 1077 153 141 3.0 18.3
1 1rc 100.00 9 0 0 1120 1128 396 388 3.0 18.3
1 1rc 100.00 9 0 0 620 628 503 511 3.0 18.3
1 1rc 100.00 8 0 0 781 788 70 77 12 16.4
1 1rc 100.00 8 0 0 625 632 206 199 12 16.4
1 1rc 100.00 8 0 0 532 539 276 283 12 16.4
1 1rc 100.00 8 0 0 322 329 381 374 12 16.4
1 1rc 91.67 12 1 0 1053 1064 418 429 12 16.4
1 1rc 100.00 8 0 0 588 595 437 444 12 16.4
1 1rc 100.00 8 0 0 442 449 478 485 12 16.4
1 1rc 100.00 8 0 0 861 868 570 577 12 16.4
1 1rc 100.00 8 0 0 610 617 630 623 12 16.4
and the list continues..
I think sequence masking is disabled by default? I cannot find an option to do so
By default, sequence masking is enabled in BLAST. To disable it, add these two arguments: -soft_masking false -dust no.
Does this answer your question?
Log in to answer this question.
can't you run BLAST with reduced word_length?
https://www.ncbi.nlm.nih.gov/books/NBK279684/
word_size blastn integer 11 Length of initial exact match.