2019 blast Practice example
The location where gene prediction of pearl Exoc7 seems to be wrong
First of all file preparation
echo ">exoc7_pearl atgactctgaaggatgccctcaacaaaagtcacacaaatacaggaaacatgctcacaata cttcaaagctttgaaaatcgtttaaagaagttagagggaacagttgagcctgtttacaat gagacagaaatgctgcggcgcagacaagaaaatatagagaaaactatgacaacactggac aatgtgctgggttactaccatattgctaaagatgtacaagatttgattaaagaaggtcca gtagtttgtggtctggagaagtacctgtctactatggaccggctgctccaagcactgaac tactttaataaacataacccaaccagtctggaagtgacagatgtcatcaaagtatatgat gatggtaaagatacattgaatgcagagttccgtagtttacttggtcgtcactgtcgtccg gtgccggctgttactatactggatttactaggaccagatgaagagttacaaacaatggaa aatgatgcacccatagaacatctgcctgagaaaattgtgaatgatttaaccctcatcgca aagtggctatacaccaatggtaaagctacagagtatatgaaagattacaccaaagtcagg tcccaaatgctcctctactctctgcaggggaactcaataaagcggaaggctaccacggcc ttgatgcagtccccttttgatccaggtcatagaagacaaggctcttataacgaattgaca aaagaggaaagttttgatgttgaaattgatatctacataacagaactaacagcattgctg aaacttattcagaatgaccctgagagatcttcgatgccccgagacggtacagttcatgaa ctgacaaaccataccattatagtactggagcccctgttagattatgctgagacagctggg gccatgttactcacccatggtgaacatgcagttccatctgatgctgtggatgtcaagaaa agtaaactcaagttggctgactatatcactaaggttttgtcagcattaggattaaactta agtaacaaggcagaaacttacagtgatccaatactcagacatgtgttcatgcttaataac tatcactacatactcaagtctttaaaaaggtctggggtattagaattaattcacacatgg aataaagatgtaggacagttttatgaggaccagatacatgaacaaaaaagactttattcc cagagctggagtaaagttctacattttgtactggaaatgaatgagccaatatcccaacaa agaatccagcaaatggagacatcaaagataaaggacaaagaaaagcagaatataaaagac aagttctctggattcaacaaagagttggaagaaatctcacgtgttcagaaagcatacgcc attcctgatccagaactgagggacaatatcaagaaagacaataaagaatatattgtgccg cgatacaagcttttcttagaaaaatttcaacggctgaacttcacaaagaattcagaaaaa tatatgaaatacactgtaaaggatgtggaagaaacacttgataaatttttcgatacttca gcttaa" > exoc7_pearl.fna
echo ">EKC30356.1 Exocyst complex component 7 [Crassostrea gigas] MLTILQSFENRLRKLENTVEPVYNETEMLRRRQENIEKTMVTLDNVLGYYHVGKEVEEFIKEGPHNCGLE KYLSIMDRLVQAHNYFNKHNPTSLELTDVIRVYDDGKEALVIEFRTLLGRHCRPVPPVMVLDMISTDEEL QGSDDIQLEHLPEKILTELSLISTWLFNNTKNTEYMKDYTRSRSSMLIKSLQGHSFKRRAVITLMQSPFD PGNKRQGSHAELPKEENLDVEVDIYITELSALLKLIQSEAQLMSGIIADKHHRSVFDNIIQEGLDSVIKN GELLAVNAKKSIAKHDFINVLSVFPVLKHLRSIKPEFDLTLEGCATPTRAKLTSLLSTLGSTAAKALEEF ALSIKTDPEKASMPKDGTVHELTNRTIIFLEPLQDYADTAGAMLLLHGEQAAPSEAVDPKKSKMRLADYI TKTLSALGLNLTIKAETYSDPTLRPVFMLNNYHYILKSLKRSGLLDLIHTWNKDVGQFYEDRINEQKKLY SESWSRVMHYITEVHEPISQQRIQAMENSKLKDKEKQNIKDKFSGFNKELEDILKIQKGYAIPDPELREQ MKKDNKDFIIPAFRMFLDKFKRLNFTKNPEKYIKYSVQDVAEVVDKLFDMSA" > exoc7_c.gigas.faa
Download the genome of the pearl oyster which you will use later.
wget https://marinegenomics.oist.jp/pearl/download/pfu_genome1.0.fasta.gz gzip -d pfu_genome1.0.fasta.gz
When you're ready, start docker to set up an environment that uses blast.
docker run -it --rm -v $PWD:$PWD -w $PWD quay.io/biocontainers/blast:2.7.1--boost1.64_1 bash
Entering the above commands will enter the virtual environment available to blast. The amino acid sequence of exoc7_c.gigas is used as a DB, and the nucleotide sequence of exoc7_pearl is used as a query. First, create a database of exoc7_c.gigas.
makeblastdb -in exoc7_c.gigas.faa -dbtype prot
Next, blastx will be performed with the amino acid sequence of exoc7_c.gigas as DB, and the DNA sequence of exoc7_pearl as a query.
blastx -db exoc7_c.gigas.faa -query exoc7_pearl.fna -num_threads 4
As a result, the following output is obtained.
Database: exoc7_c.gigas.faa
1 sequences; 612 total letters
Query= exoc7_pearl
Length=1566
Score E
Sequences producing significant alignments: (Bits) Value
EKC30356.1 Exocyst complex component 7 [Crassostrea gigas] 426 6e-148
> EKC30356.1 Exocyst complex component 7 [Crassostrea gigas]
Length=612
Score = 426 bits (1096), Expect = 6e-148, Method: Compositional matrix adjust.
Identities = 196/259 (76%), Positives = 235/259 (91%), Gaps = 0/259 (0%)
Frame = +1
Query 787 IQNDPERSSMPRDGTVHELTNHTIIVLEPLLDYAETAGAMLLTHGEHAVPSDAVDVKKSK 966
I+ DPE++SMP+DGTVHELTN TII LEPL DYA+TAGAMLL HGE A PS+AVD KKSK
Sbjct 354 IKTDPEKASMPKDGTVHELTNRTIIFLEPLQDYADTAGAMLLLHGEQAAPSEAVDPKKSK 413
Query 967 LKLADYITKVLSALGLNLSNKAETYSDPILRHVFMLNNYHYILKSLKRSGVLELIHTWNK 1146
++LADYITK LSALGLNL+ KAETYSDP LR VFMLNNYHYILKSLKRSG+L+LIHTWNK
Sbjct 414 MRLADYITKTLSALGLNLTIKAETYSDPTLRPVFMLNNYHYILKSLKRSGLLDLIHTWNK 473
Query 1147 DVGQFYEDQIHEQKRLYSQSWSKVLHFVLEMNEPISQQRIQQMETSKIKDKEKQNIKDKF 1326
DVGQFYED+I+EQK+LYS+SWS+V+H++ E++EPISQQRIQ ME SK+KDKEKQNIKDKF
Sbjct 474 DVGQFYEDRINEQKKLYSESWSRVMHYITEVHEPISQQRIQAMENSKLKDKEKQNIKDKF 533
Query 1327 SGFNKELEEISRVQKAYAIPDPELRDNIKKDNKEYIVPRYKLFLEKFQRLNFTKNSEKYM 1506
SGFNKELE+I ++QK YAIPDPELR+ +KKDNK++I+P +++FL+KF+RLNFTKN EKY+
Sbjct 534 SGFNKELEDILKIQKGYAIPDPELREQMKKDNKDFIIPAFRMFLDKFKRLNFTKNPEKYI 593
Query 1507 KYTVKDVEETLDKFFDTSA 1563
KY+V+DV E +DK FD SA
Sbjct 594 KYSVQDVAEVVDKLFDMSA 612
Score = 389 bits (998), Expect = 2e-133, Method: Compositional matrix adjust.
Identities = 203/341 (60%), Positives = 250/341 (73%), Gaps = 33/341 (10%)
Frame = +1
Query 49 MLTILQSFENRLKKLEGTVEPVYNETEMLRRRQENIEKTMTTLDNVLGYYHIAKDVQDLI 228
MLTILQSFENRL+KLE TVEPVYNETEMLRRRQENIEKTM TLDNVLGYYH+ K+V++ I
Sbjct 1 MLTILQSFENRLRKLENTVEPVYNETEMLRRRQENIEKTMVTLDNVLGYYHVGKEVEEFI 60
Query 229 KEGPVVCGLEKYLSTMDRLLQALNYFNKHNPTSLEVTDVIKVYDDGKDTLNAEFRSLLGR 408
KEGP CGLEKYLS MDRL+QA NYFNKHNPTSLE+TDVI+VYDDGK+ L EFR+LLGR
Sbjct 61 KEGPHNCGLEKYLSIMDRLVQAHNYFNKHNPTSLELTDVIRVYDDGKEALVIEFRTLLGR 120
Query 409 HCRPVPAVTILDLLGPDEELQTMENDAPIEHLPEKIVNDLTLIAKWLYTNGKATEYMKDY 588
HCRPVP V +LD++ DEELQ +D +EHLPEKI+ +L+LI+ WL+ N K TEYMKDY
Sbjct 121 HCRPVPPVMVLDMISTDEELQG-SDDIQLEHLPEKILTELSLISTWLFNNTKNTEYMKDY 179
Query 589 TKVRSQMLLYSLQGNSIKRKATTALMQSPFDPGHRRQGSYNELTKEESFDVEIDIYITEL 768
T+ RS ML+ SLQG+S KR+A LMQSPFDPG++RQGS+ EL KEE+ DVE+DIYITEL
Sbjct 180 TRSRSSMLIKSLQGHSFKRRAVITLMQSPFDPGNKRQGSHAELPKEENLDVEVDIYITEL 239
Query 769 TALLKLIQNDPERSSMPRDGTVHELTNHTII--VLEPLLDYAETAGAMLLTHGEHAVPSD 942
+ALLKLIQ++ + S G + + + ++ +++ LD G +L
Sbjct 240 SALLKLIQSEAQLMS----GIIADKHHRSVFDNIIQEGLDSVIKNGELL----------- 284
Query 943 AVDVKKSKLKLADYITKVLSALGLNLSNKAETYSDPILRHV 1065
AV+ KKS K D+I VLS P+L+H+
Sbjct 285 AVNAKKSIAK-HDFIN-VLSVF-------------PVLKHL 310
There are two aligned sequences, and looking at the second alignment, “Sbjct” is the sequence of exoc7_c.gigas used as DB and is homologous to exoc7_pearl, but the homology after 248 amino acids is low. Since exoc7_pearl of the query is a DNA sequence, it increases by 3 bp per character, and the region of 787 bp or later corresponds to a region of low homology. Then, when looking at the first alignment, the alignment of exoc7_pearl starts from 787 bp, while the exoc7_c.gigas starts from 354 amino acids. Therefore, it is understood that the 248 to 353 amino acids of exoc7_c.gigas are deleted in predicted exoc7_pearl.
Investigate which scaffold of the pearl oyster genome has the deleted sequences in gene prediction.
First, the deleted sequence is stored in FASTA format.
echo ">exoc7_deleted SEAQLMSGIIADKHHRSVFDNIIQEGLDSVIKNGELLAVNAKKSIAKHDFINVLSVFPVLKHLRSIKPEFDLTLEGCATPTRAKLTSLLSTLGSTAAKALEEFALS" > exoc7_deleted.faa
Next, create a blast index of the pearl oyster genome.
makeblastdb -in pfu_genome1.0.fasta -dbtype nucl
Search the pearl oyster genome with tblastn, using the deleted sequence as a query.
tblastn -db pfu_genome1.0.fasta -query exoc7_deleted.faa -num_threads 4
As a result,
Query= exoc7_deleted
Length=106
Score E
Sequences producing significant alignments: (Bits) Value
scaffold294819.1|size544 118 2e-33
> scaffold294819.1|size544
Length=544
Score = 118 bits (295), Expect = 2e-33, Method: Compositional matrix adjust.
Identities = 55/75 (73%), Positives = 67/75 (89%), Gaps = 0/75 (0%)
Frame = -3
Query 1 SEAQLMSGIIADKHHRSVFDNIIQEGLDSVIKNGELLAVNAKKSIAKHDFINVLSVFPVL 60
SEAQLMSGII +KHHRSVF++II+ LD V+K GE LA NAKKSI+KHDF++VLSVFPV+
Sbjct 227 SEAQLMSGIIPEKHHRSVFESIIEGSLDMVVKGGETLASNAKKSISKHDFLSVLSVFPVV 48
Query 61 KHLRSIKPEFDLTLE 75
+HLR++KPEFDL LE
Sbjct 47 RHLRTVKPEFDLALE 3
is obtained, and it can be seen that there is one Exoc7 exon in the scaffold 29491. 1 (Exon of Exoc7 other than this exon is in scaffold 1200).