Create fasta of matched sequences with gaps
1
0
Entering edit mode
3.2 years ago
kingcohn ▴ 20

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.

alignment gene blast nucleotides fasta • 1.1k views
ADD COMMENT
0
Entering edit mode

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?

ADD REPLY
0
Entering edit mode
3.0 years ago

If you run your blast and request tab output format as such: -outfmt '6 std qseq sseq' you will get the sequence back as separate fields in the tabulated files, which then can be easily parsed out.

How do you get to the 1006 number btw?

ADD COMMENT

Login before adding your answer.

Traffic: 1062 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6