This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to force blastn to align full subject

My aim is to dispatch query sequences from fasta file using blastn according to 2 subject sequences.

Here, an example of a query sequence :

>M02945:227:000000000-BHPRH:1:2114:12657:28629 1:N:0:1

GAGAGCTCTATTACTGCACTGGAGGGTAGC

Here, my 2 subjects sequences, they are both 30bp long :

>ESR

GAGAGCTCTATTTCTTCGGTGTTACCAGCT

>Cercle

GAGAGCTCTATTACTGCACTGGAGGGAAGC

By look, the sequence M02945:227:000000000-BHPRH:1:2114:12657:28629 1:N:0:1 will go to Cercle.fa file (29 first bases match with 1 mismatch position 27).

And using blastn :

blastn -query test_sequence_for_blast.fa -subject ref.fa -culling_limit 1 -task blastn-short -outfmt 6

M02945:227:000000000-BHPRH:1:2114:12657:28629 ESR 100.00 8 0 0 2 9 9 2 0.008 16.4

M02945:227:000000000-BHPRH:1:2114:12657:28629 Cercle 100.00 26 0 0 1 26 1 26 1e-13 52.0

My problem is, I would like to know how many mismatch do blastn found in the full length subject sequence. So, first try I decrease the mismatch penalty to -1 (because here I just have the 1->26 information).

blastn -query test_sequence_for_blast.fa -subject ref.fa -culling_limit 1 -penalty -1 -task blastn-short -outfmt 6

M02945:227:000000000-BHPRH:1:2114:12657:28629 ESR 100.00 8 0 0 2 9 9 2 0.027 14.3

M02945:227:000000000-BHPRH:1:2114:12657:28629 Cercle 96.67 30 1 0 1 30 1 30 8e-12 46.0

Nice, I got the 1->30 information, ok that's great.

I inserted a second mismatch position 29

>M02945:227:000000000-BHPRH:1:2114:12657:28629 1:N:0:1

GAGAGCTCTATTACTGCACTGGAGGGTACC

blastn -query test_sequence_for_blast.fa -subject ref.fa -culling_limit 1 -penalty -1 -task blastn-short -outfmt 6

M02945:227:000000000-BHPRH:1:2114:12657:28629 ESR 100.00 8 0 0 2 9 9 2 0.008 16.4

M02945:227:000000000-BHPRH:1:2114:12657:28629 Cercle 100.00 26 0 0 1 26 1 26 1e-13 52.0

And again I lost the info from 26->30... Because at the end I want to totally send to trash sequences with more than 2 mismatches

How to force blastn to align 1->30 straigth and yield the number of mismatch ?

I also tried to modify gap open penalty and gap extend penalty but that didn't do the trick...

Thanks for helping !

blastn
• 3,499 views • link
updated by Biostar Biostar • written by Bastien Hervé Bastien Hervé

1 answer

BLAST, by definition, does only LOCAL alignments. Maybe you are looking for a global alignment tool https://en.wikipedia.org/wiki/List_of_sequence_alignment_software#Pairwise_alignment

• 0 views • link
written by Santosh Anand Santosh Anand

Thanks for this help, actually my query sequences are 100-150 bp long, I just trimmed the example sequence here for a better understanding. I think i'll substring the 30 first bases of each sequences and use bowtie2 --end-to-end on them.

• 0 views • link
written by Bastien Hervé Bastien Hervé

Why substring the first 30 bases? Bowtie2 or bwa mem should be supporting longer reads

• 0 views • link
written by Santosh Anand Santosh Anand

I'm just interested in the first 30 bases

• 0 views • link
written by Bastien Hervé Bastien Hervé

Nevertheless, don’t use a local aligner if you want a full sequence alignment.

You could also look at just calculating Hamming or Levenshtein distances if your sequences are pretty similar and don’t need aligning first

• 1 views • link
written by Joe Joe

Log in to answer this question.

More posts like this