Поиск некодирующих последовательностей программой BLASTN

Определим, сколько гомологов каждой из тРНК найдет программа BLAST в геноме родственной бактерии Pasteurella multocida. Для этого выполним следующие действия: 
formatdb -i pm_genome.fasta -n pm -p F
blastall -p blastn -d pm -i trna_ecoli.fasta -m8 -o blastn.fasta (поиск по созданному нами банку)
grep ">" trna_ecoli.fasta (создание колонки из названий входных последовательностей)
С помощью скрипта script.scr было получено число находок для каждой последовательности.
Результат можно увидеть в файле trna.xls (1-2 колонки).

Далее был повторен поиск, с порогом на E-value равным 0.001. 
blastall -p blastn -d pm -i trna_ecoli.fasta -m8 -e 0.001 -o blastne.fasta
Результат можно увидеть в файле trna.xls (3 колонка).

Поиск некодирующих последовательностей программой megablast

Поиск гомологов тРНК программой megablast: 
megablast -d pm -i trna_ecoli.fasta -m8 > megablast.fasta
Поиск программой discontigous megablast: 
megablast -d pm -i trna_ecoli.fasta -m8 -D 2 -t 16 -W 11 -N 1 > dismegablast.fasta
где 
-D - тип выдачи (2 - стандартная)
-t - длина искомого слова с учетом разрывов (мною выбрана наименьшая)
-W - длина поискового слова без учета разрывов (мною также выбрано наименьшее значение)
-N - тип поисковых слов для поиска по некодирующим последовательностям
Результат можно увидеть в файле trna.xls (4-5 колонки).

Анализ результатов

При поиске гомологичных участков в геноме бактерии Pasteurella multocida для тРНК valV программа BLASTN находит последовательность AE006104, а программа MegaBLAST - нет.
Причина заключается в том, что BLASTN ищет по словам длины 11, а MegaBLAST - длины 28, поэтому при отсутствии 28 совпавших нуклеотидов, расположенных подряд без пропусков, MegaBLAST может пропустить данный участок. Программа discontigous megablast была запущена с длиной слова равной 16 (с учетом разрывов), но все равно не нашла данный участок.

Вырежу гомологичный участок в отдельный файл: ae006104.fasta 
Input (gapped) sequence(s): pm_genome.fasta:AE006104
     Begin at position [start]: 8727
       End at position [end]: 8747
        Reverse strand [N]: y
output sequence(s) [ae006104.fasta]:
Также выделю исходную последовательносмть в отдельный файл: valv.fasta
Выровняю эти последовательности программой needle: 
########################################
# Program: needle
# Rundate: Sun 22 Nov 2009 22:41:08
# Commandline: needle
#    [-asequence] valv.fasta
#    [-bsequence] ae006104.fasta
#    [-outfile] valae.needle
# Align_format: srspair
# Report_file: valae.needle
########################################

#=======================================
#
# Aligned_sequences: 2
# 1: valV
# 2: AE006104
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 77
# Identity:      21/77 (27.3%)
# Similarity:    21/77 (27.3%)
# Gaps:          56/77 (72.7%)
# Score: 105.0
# 
#
#=======================================

valV               1 gcgttcatagctcagttggttagagcaccaccttgacatggtgggggtcg     50
                            |||||||||||||||||||||                      
AE006104           1 -------tagctcagttggttagagcac----------------------     21

valV              51 ttggttcgagtccaattgaacgcacca     77
                                                
AE006104          21 ---------------------------     21
Программа MegaBLAST не нашла данный участок, т.к. в нем нет 28 подряд идущих совпавших нуклеотидов (только 21).

Аннотация участка записи AE006104 EMBL (замещена записью AE004439 19-FEB-2009, последовательность найдена в результате выравнивания blastn: 
>gb|AE004439.1|  Pasteurella multocida subsp. multocida str. Pm70, complete genome
Length=2257487


 Score = 42.1 bits (21),  Expect = 6e-06
 Identities = 21/21 (100%), Gaps = 0/21 (0%)
 Strand=Plus/Plus

Query  1       TAGCTCAGTTGGTTAGAGCAC  21
               |||||||||||||||||||||
Sbjct  757583  TAGCTCAGTTGGTTAGAGCAC  757603 ) :

FT   gene            complement(757549..757761)
FT                   /locus_tag="PM0654"
FT   CDS             complement(757549..757761)
FT                   /codon_start=1
FT                   /transl_table=11
FT                   /locus_tag="PM0654"
FT                   /product="unknown"
FT                   /db_xref="GOA:Q9CMZ6"
FT                   /db_xref="UniProtKB/Swiss-Prot:Q9CMZ6"
FT                   /protein_id="AAK02738.1"
FT                   /translation="MLDSIGMFWKDLSVIYGFLIDFIGDFGLFWIGLESEMVRLAGLEP
FT                   VTPTMSRWCSNQLSYSRILKCGAII"
Как мы видим, product="unknown" вряд ли гомологичен tRNA-Val, скорее всего кодирующие их нуклеотидные последовательности просто имеют похожий участок.

Назад