Hello. I am attempting to isolate exons in thousands of genes using a highly fragmented draft assembly and a well curated official gene set of a closely related species. So far, I've gotten reasonable data from the program SPALN, which takes exonic regions and concatenates these in a closely related species. However, to continue with branch-site tests for selection using PAML I should have the codon files be of equal length. I am having trouble creating a fasta file with my query exons of the same length as my subject as the fragmented assemblies have large gaps/poor quality. Not sure how to phrase my problem adequately, so please look at my blast output below:
Query= scaffold5976_size9010.3 scaffold5976_size9010 + [1:9010] ( 1384 -
5266 ) Ldec\Orco + 1:1440 ( 1127 - 1440 ) N 627.00;C
join(1384..1490,5002..5157,5216..526)
Length=314
Subject= Ldec\Orco
Length=1440
Score = 569 bits (308), Expect = 1e-166
Identities = 312/314 (99%), Gaps = 0/314 (0%)
Strand=Plus/Plus
Query 1 AGATCAACGGAGTTACCGTGTATGCTGCTACTGTGATAGGTTACCTGGTGTATTCTTTGG 60
|||||| |||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 1127 AGATCAGCGGAGTTACCGTGTATGCTGCTACTGTGATAGGTTACCTGGTGTATTCTTTGG 1186
Query 61 CCCAGGTATTCCATTTCTGCATTTTTGGGAACAGGCTGATAGAGGAGAGTTCATCTGTTA 120
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 1187 CCCAGGTATTCCATTTCTGCATTTTTGGGAACAGGCTGATAGAGGAGAGTTCATCTGTTA 1246
Query 121 TGGAAGCAGCTTACAGCTGTCACTGGTATGATGGTTCAGAGGAAGCGAAAACATTCGTCC 180
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 1247 TGGAAGCAGCTTACAGCTGTCACTGGTATGATGGTTCAGAGGAAGCGAAAACATTCGTCC 1306
Query 181 AGATTGTATGTCAACAATGTCAAAAAGCCTTGTCGATATCTGGGGCGAAGTTTTTCACTA 240
||||||||||||||||||||||||||||||||||||||||||||||||||||||| ||||
Sbjct 1307 AGATTGTATGTCAACAATGTCAAAAAGCCTTGTCGATATCTGGGGCGAAGTTTTTTACTA 1366
Query 241 TTTCTCTAGATCTTTTTGCCTCGGTACTTGGTGCAGTAGTTACATATTTCATGGTACTGG 300
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 1367 TTTCTCTAGATCTTTTTGCCTCGGTACTTGGTGCAGTAGTTACATATTTCATGGTACTGG 1426
Query 301 TACAACTCAAATAA 314
||||||||||||||
Sbjct 1427 TACAACTCAAATAA 1440
Lambda K H
1.33 0.621 1.12
Gapped
Lambda K H
1.28 0.460 0.850
Effective search space used: 431256
Query= scaffold63454_size2903.1 scaffold63454_size2903 + [1:2903] ( 1043 -
1215 ) Ldec\Orco + 1:1440 ( 1 - 173 ) N 358.00;C join(1043..121)
Length=173
Subject= Ldec\Orco
Length=1440
Score = 320 bits (173), Expect = 7e-92
Identities = 173/173 (100%), Gaps = 0/173 (0%)
Strand=Plus/Plus
Query 1 ATGATGAAATTCAAGGTATCCGGCCTTGTGGCTGACCTTATGCCTAACATAAGGCTCATA 60
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 1 ATGATGAAATTCAAGGTATCCGGCCTTGTGGCTGACCTTATGCCTAACATAAGGCTCATA 60
Query 61 CAGGCATCGGGTCACTTCATGTTCAATTATCATGCTGATAATTCAGGAGCACTCCATGCC 120
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 61 CAGGCATCGGGTCACTTCATGTTCAATTATCATGCTGATAATTCAGGAGCACTCCATGCC 120
Query 121 TTACGACTGGGATATTCTTGCATGCATCTAGTTTTCTGTTTAGTGCAATACGG 173
|||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 121 TTACGACTGGGATATTCTTGCATGCATCTAGTTTTCTGTTTAGTGCAATACGG 173
Lambda K H
1.33 0.621 1.12
Gapped
Lambda K H
1.28 0.460 0.850
Effective search space used: 231498
As you can see the second match is the beginning of the subject gene, query and subject are both 1, this alignment continues for 173 bases then there is a huge gap of 1006 nucleotides, with respect to the subject, and the alignment finishes from 1127-1427. I would like to create a fasta file with the query alignment from 1-173 followed by 1006 N's, and ending again with the query that aligned to the reference from 1127-1427 and do this for all 20,000 genes in the official gene set.
doesn't the blast output already gives you the exons of equal length? Why would you then still want to fill in the 'gap' with Ns?