Депонирование контигов в Whole Genome Shotgun (WGS) базу данных

Депонирование результатов полногеномного секвенирования является важным элементом научно-исследовательской работы. На сайте NCBI имеется достаточно много документации описывающей необходимые шаги для успешной публикации ваших геномов в генбанке. Однако, как показала практика, так или иначе приходится сталкиваться с различного рода проблемами, явно не указанными в документации. Итак, оставим философию и перейдем непосредственно к описанию процесса подготовки ваших нуклеотидных последовательностей к депонированию. Хотелось бы отметить, что цель данной заметки заключается в описании процесса подготовки ваших sequences к депонированию, а не в объяснении базовых принципов структуризации информации на сайте NCBI т.е. вы должны иметь представление о таких понятиях как BioProject, BioSample и т.д. Итак поехали…

Continue reading «Депонирование контигов в Whole Genome Shotgun (WGS) базу данных»

Депонирование контигов в Whole Genome Shotgun (WGS) базу данных

Сборка генома и аннотирование снипов (pipeline скрипт)

Всех биоинформатиков можно разделить на две большие группы: пользователей и писателей. Однако, в отличие от классического программирования, где пользователи практически никогда не участвуют в процессе создания программ, в биоинформатике пользователи и писатели разделены достаточно условной стеной. Скорее всего это связано с тем, что не существует единого универсального инструмента для работы с геномами. Нет, конечно есть удобные «комбайны» с удобным интерфейсом и интуитивным подходом, но они в большинстве своем платные. Более того, предпочтение работы в той или иной программе, снижает потребность специалиста к поиску новых инструментов работы, делает его зависимым от того или иного программного обеспечения. Большинство биоинформатиков используют в своей работе большое количество разнообразных программ и запускать каждый раз последовательно их друг за другом, мягко говоря утомительно и наступает именно тот самый момент, когда биоинформатик-пользователь становится чуть-чуть писателем: он садится и пишет скрипт для управления последовательностью запуска используемых в работе программ. Такие скрипы называются pipeline-скрипты. Они могут быть написаны на любом знакомом пользователю языке, часто это bash, perl, python, ruby.
Мы например, для своего pipeline-скрипта, используем язык программирования bash. Итак, что же включает в себя наш pipeline? Для удобства работы, мы разделили скрипт на две части: исполняемую и настроек. Если исполняемая часть фактически остается неизменной, то часть настроек, вынесенная в отдельный файл, активно редактируется. Это действительно удобно, для различных микроорганизмов, хранить разные файлы настроек, тем самым избавляя себя от путаницы и хаоса.
Continue reading «Сборка генома и аннотирование снипов (pipeline скрипт)»

Сборка генома и аннотирование снипов (pipeline скрипт)

Как выбрать файлы в одной папке и скопировать их в другую папку?

Дано: В исходной папке у нас находятся файлы с расширением ext1. Нам нужно, файлы, с таким же именем, но расширением ext2, находящиеся в другой папке, скопировать в третью папку.

#Создаем список файлов из исходной папки
find . -name "*.ext1" -printf '%f\n' > file_names
# Меняем расширение ext1 на ext2
sed -i 's/.ext1/.ext2/g' file_names
# Копируем файл со списком имен файлов в нужную нам папку
cp file_names ~/path/folder2 | cd ~/path/folder2
# Копируем файлы из folder2 в folder3 в соответствии со списком файлов
cat file_names | xargs -I % cp % folder3

Как выбрать файлы в одной папке и скопировать их в другую папку?

Как получить консенсусную последовательность из коротких прочтений (reads)

Для этого воспользуемся программой Samtools и немного своего кода для awk:

1. Для работы нам понадобится т.н. pileup файл. Данный файл получаем из bam файла при помощи следующей команды:

samtools mpileup -f REFERENCE_FILE SORTED_BAM_FILE > OUR_PILEIP_FILE.pileup

где REFERENCE_FILE — файл референсного генома, SORTED_BAM_FILE — сортированный bam файл, OUR_PILEIP_FILE.pileup — итоговый pileup файл

2. Воспользуемся программой awk (входит по-умолчанию в любой дистрибутив Linux)

awk '{if ($5 !~ /\.|\d+|\,|\+|\*|\$|N/ && length($5)>2) printf ("%s",toupper(substr($5,1,1)));else printf ("%s",$3)}' OUR_PILEIP_FILE.pileup > OUR_SEQUENCE_FILE.fasta

Вуаля! Наша последовательность сохранена в файл OUR_SEQUENCE_FILE.fasta

Как получить консенсусную последовательность из коротких прочтений (reads)

Проверка наличия снипа в vcf файлах.

Дано: У нас есть группа штаммов, которая после сиквенса, snp-calling-а и аннотирования объединена в одной папке в виде vcf файлов. Также, есть файл-таблица, в котором аккумулированы данные об уникальных снипах для данной группы (3-й столбец позиция снипа).
Необходимо: Для того, чтобы быть окончательно уверенным, что все снипы в таблице действительно уникальны, следует проверить их наличие во всех vcf файлах группы. Скрипт запусакем: ./script.sh таблица_со_снипами.csv. Снипы, которые содержатся во всех vcf файлах группы, будут сохранены в файл exact_list, не прошедшие проверку снипы будут сохраняться в файл fail_list.

#!/bin/bash
workdir=`pwd`
rm fail_list
rm exact_list
files_count=$(ls -l *.eff.vcf | wc -l)
printf "%s:%s\n" $1 $files_count
awk -F "\t" '{print $3}' $1 > uniq_pos
while read -r line ; do
count_words=$(grep -rlP "${line}\t" *.eff.vcf | wc -l)
if [ $count_words == $files_count ]; then
printf "%s\t%s\t%s\tOK\n" $line $count_words $files_count
printf "%s\t%s\t%s\n" $line $count_words $files_count >> exact_list
else
printf "%s\t%s\t%s\tFAIL\n" $line $count_words $files_count
fi
if [ $count_words == 0 ]; then
printf "%s\t%s\t%s\n" $line $count_words $files_count >> fail_list
fi
done < uniq_pos rm uniq_pos

Проверка наличия снипа в vcf файлах.

Чтение csv файла, выбор и сортировка данных на основе ключевого слова

В продолжение решения задачи, рассмотренной раннее , требуется выбрать несинонимичные замены в генах, относящихся к системе токсин-антитоксин.


#!/bin/bash
# удаляем старый файл
rm tox_atox_result.txt
# отбираем записи, содержащие значения Tox и aTox и сохраняем нужные нам столбцы во временный файл
awk -F ";" '{print}' *.snp.csv|grep -E "\saTox|\sTox" |awk -F ";" '{if (length($12)!=0) print $1"\t"$2"\t"$9"\t"$10"\t"$12"\t"$13"\t"$14"\t"$15}' > tmp.csv
# сортируем столбцы
sort -t\t tmp.csv > tmp_s.csv
# удаляем повторяющиеся записи
uniq tmp_s.csv > tox_uniq.txt
# определяем кол-во файлов в папке
files_count=$(ls -l *.eff.snp.csv | wc -l)
# построчно считываем значения временного файла
while read -r line; do
# указываем тип разделения записей-табуляция
IFS=$'\t'
# присваиваем переменной значения столбцов таблицы
columns=($line)
# считаем, сколько файлов содержат уникальную нуклеотидную позицию (столбец 2)
count_words=$(grep -Rcl -e ${columns[2]} *.snp.csv |wc -l)
# если общее кол-во файлов = количеству файлов с уникальной нуклеотидной позицией, то
if [ $count_words == $files_count ]; then
printf '%s/%s:%s\n' "$files_count" "$count_words" "$line"
# сохраняем результат в файл
printf '%s\n' "$line" >> tox_atox_result.txt fi
done

Чтение csv файла, выбор и сортировка данных на основе ключевого слова

Простые решения на Ruby

Как получить комплементарную последовательность ДНК?


DNA = 'AAAACCCGGT'
cDNA = DNA.tr("AGTC","TCAG")
puts cDNA

Как получить комплементарную последовательность РНК?


DNA = 'AAAACCCGGT'
RNA= DNA_compl.tr("AGTC","UCAG")
puts RNA

Как посчитать количество каждого отдельного нуклеотида в последовательности?

DNA = 'AAAACCCGGT'
puts DNA.count "A"
puts DNA.count "C"
puts DNA.count "G"
puts DNA.count "T"

Простые решения на Ruby

Поиск и сортировка уникальных значений из csv файла

Одна из наших программ (скрипт), сохраняет данные из аннотированного vcf файла в csv формат, где разделителем служит (;). Выглядит это так:


gene; pos; ref; alt; qual; effect; effect_impact; functional_class ; codon_change; aa_change; aa_length ;tang_index; mdr_genes; unique; virulence; annotation; note
Rv0003;4013;T;C;179.0;NON_SYNONYMOUS_CODING;MODERATE;MISSENSE;aTc/aCc;I245T; 385;0.750;;NO;;DNA replication and repair protein RecF (single-strand DNA binding protein) (recF) regulated_by: (Rv0001,Rv1985c,Rv2711,Rv3286c,Rv3676) ;RecF, DNA replication and repair protein (see citations below),equivalent to other mycobacterial DNA replication and repair proteins. Also highly similar to many others. Contains PS00017 ATP/GTP-binding site motif A (P-loop),PS00617 RecF protein signature 1, and PS00618 RecF protein signature 2. Belongs to the RecF family.
Rv0006;7362;G;C;127.0;NON_SYNONYMOUS_CODING;MODERATE;MISSENSE;Gag/Cag;E21Q;838; 1.634;FLQ;NO;;DNA gyrase (subunit A) GyrA (DNA topoisomerase (ATP-hydrolysing)) (DNA topoisomerase II) (type II DNA topoisomerase) (gyrA) regulated_by: (Rv3648c,Rv3676,Rv0348) ;GyrA, DNA gyrase subunit A (see citations below). Contains PS00018 EF-hand calcium-binding domain.

Continue reading «Поиск и сортировка уникальных значений из csv файла»

Поиск и сортировка уникальных значений из csv файла

Первое испытание начинающего биоинформатика

Тема данной заметки, несмотря на свою игривость, несет в себе огромный философский смысл. Основной посыл данного сообщения нацелен на тех людей, которые работали и работают на своих компьютерах под управлением операционной системы семейства Windows компании Microsoft, а таких людей по-статистике около 90% и только начинают свои первые шаги в познании основ биоинформатики.

Continue reading «Первое испытание начинающего биоинформатика»

Первое испытание начинающего биоинформатика