Поиск в геноме участков, кодирующих белки, похожие на заданный
Используя Standalone BLAST, попытаемся найти гомологи нашего белка YPWA_BACSU в геноме Geobacillus thermodenitrificans. Для этого выполним следующую команду в рабочей директории с последовательностью генома в fasta-формате и последовательностью нашего белка:makeblastdb -in gt_genome.fasta -dbtype nucl
Будем использовать алгоритм tblastn - подходящий для наших целей - с параметром e-value 0,001:tblastn -query YPWA_BACSU.fasta -db gt_genome.fasta -out homologyes.out -evalue 0.001
Результатом стал файл homologyes.out, в котором содержалась одна запись.| Число находок с E-value < 0,001 | 1 |
| E-value лучшей находки | 0.0 |
| Название последовательности с лучшей находкой | CP000557 CP000557.1 Geobacillus thermodenitrificans NG80-2, complete genome. |
| Координаты лучшей находки (от-до) | 1483914-1485374 |
| Доля последовательности вашего белка, вошедшая в выравнивание с лучшей находкой |
(488-3+1)/501=0,97 |
Поиск гомологов некодирующих последовательностей программой BLASTN
Определим, сколько гомологов каждой тРНК (из файла trna_bacsu.fasta) содержится в геноме бактерии Geobacillus thermodenitrificans. Для этого запустим BLSTN со следующими парамтрами:blastn -task blastn -query trna_bavsu.fasta -db gt_genome.fasta -out trna.out outfmt 6 -evalue 0.01
Для того,чтобы узнать, сколько находок для конкретной тРНК, выполним например, команду:grep -c 'BSn5_t20894' trna.out
Теперь сделаем то же, но не для выбранной последовательности, а для всех. Создадим колонку из названий входных последовательностей командой и импортируем ее в Excel
grep ">" trna_bacsu.fasta > trna.xls
Функцией "СЦЕПИТЬ" в соседнем столбце сцепим команду grep -c, значение столбца слева (с названиями тРНК), файл для поиска и конечный файл, распространим на соседние ячейки столбца. Получим много строк следущего типа: grep -c >BSn5_t20894 tRNA-Gln trna.out >> trnacount.txt, скопируем их в файл scriptcount.sh, сохраним в формате Unix. Далее в FARe удалим столбцы ">" и " tRNA-Gln". Конечный вариант скрипта запускаем командами:chmod +x scriptcount.sh
./scriptcount.sh
Импортируем результаты работы скрипта в Excel, чтобы получить файл, содержащий колонку с названиями и колонку с количеством находок.Поиск гомологов при изменённых параметрах программы BLASTN
Повторим предыдущее задание с измененными параметрами:blastn -task blastn -query trna_bacsu.fasta -db gt_genome.fasta -out trna2.out -outfmt 6 -evalue 0.01 -reward 5 -penalty -4 -gapopen 25 -gapextend 10
- изменили весовую матрицу
blastn -task blastn -query trna_bacsu.fasta -db gt_genome.fasta -out trna3.out -outfmt 6 -evalue 0.01 -reward 5 -penalty -4 -gapopen 25 -gapextend 10 -word_size 4
- уменьшили размер затравки до минимальной
blastn -task blastn -query trna_bacsu.fasta -db gt_genome.fasta -out trna4.out -outfmt 6 -evalue 0.01 -word_size 4
В результате в файле Excel появилось еще 3 столбца.
Анализ результатов
Наблюдается тенденция к увеличению количества хитов при изменении весовой матрицы и - еще большее - при дополнительном уменьшении длины затравки. Если жа длину затравки уменьшать без изменения весовой матрицы, то количество хитов увеличивается не так сильно. Другими словами, весовая матрица играет большую роль при отборе потенциальных кандидатур на выравнивание, нежели длина затравки, хотя и она, безусловно, очень важна.Выберем тРНК, которая нашлась при одних параметрах и не нешласьпри других; одна из находок на глутаминовую тРНК (BSn5_t20894 tRNA-Gln) присутствовала при попытке 2 и 3,но отсутствовала при попытке 1 и 4. Связано это с уже вышеизложеными причинами.
Командой
seqret gt_genome.fasta -sask
вырезаем участок, выровнившийся с глутаминовой тРНК, посмотрев при этом на границы выравнивания и сделав вывод о цепи.Затем вырезаем интересующую тРНК из trna_bacsu.fasta, выравниванием командой needle, установив паарметры по умолчанию. Получаем:
# Aligned_sequences: 2 # 1: CP000557 # 2: BSn5_t20894 # Matrix: EDNAFULL # Gap_penalty: 10.0 # Extend_penalty: 0.5 # # Length: 79 # Identity: 41/79 (51.9%) # Similarity: 41/79 (51.9%) # Gaps: 29/79 (36.7%) # Score: 142.0 # # #======================================= CP000557 1 ------atagccaagtggtaaggcagaggtctgcaaaacctttacccc-- 42 |||||||||.|||||||||..|| ||.||.||.|| BSn5_t20894 1 tgggctatagccaagcggtaaggcaatgg--------actttgactccgt 42 CP000557 43 ------cggttcgaatccgggt------- 58 .|||||||||||.|.| BSn5_t20894 43 gatcgttggttcgaatccagctagcccag 71 #--------------------------------------- #---------------------------------------
Можем заметить, что качество выравнивания отличается не очень высокими показателями. Именно поэтому данная находка присутствует только при ослабленных параметрах отбора хитов - измененная весовая матрица или она же с уменьшением длины слова.
Между тем вполе FT записи полного генома бактерии Geobacillus thermodenitrificans в районе нашей находки находится много тРНК различных типов, однако нет глутаминовой тРНК, точно совпадающей с ней по координатам. Максимально перекрывающаяся с ней (complement 2966372..2966429) рамка относится к цистеиновой тРНК (complement 2966362..2966435).
Время работы программы BLAST
Проследим, как отличается время работы команды при изменении параметров:- Значения по умолчанию:
real 0m0.543s user 0m0.260s sys 0m0.252s
- Измененная весовая матрица:
real 0m0.718s user 0m0.448s sys 0m0.244s
- Изменена весовая матрица и сокращена длина слова:
real 0m41.083s user 0m40.652s sys 0m0.248s
- Сокращена длина слова:
real 0m30.623s user 0m30.332s sys 0m0.244s