This bash script was created as a side project during my first steps in the NGS Illumina data processing and shell scripting. It consists of simple free programs and code. The framework includes QC check, filtering, alignment, processing, calling; all the steps mostly use soft cut-offs.
Sorry, but I wrote it for Russian-speaking people that are new in the field.
NB: since at that time I wanted to automatize the process, the code is overcrowded with not smart conditionings.
The commands itself are simple; if somebody interest in the script for any reason,
please write to danatyermakovich@gmail.com
####################RU######################
Релаксированный алгоритм обработки данных NGS
от fastq до vcf при помощи простых команд
- Для работы вам потребуются следующие программы определенных версий
установленные на операционной системе GNU/Linux
(наверное, где возможно, лучше скачивать binary и т.п., ибо source нужно компилировать
и у вас может что-то пойти не так)
Trimmomatic-0.36
http://www.usadellab.org/cms/?page=trimmomatic
bowtie2-2.3.3 #linux-x84_64
https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.3/
Samtools 1.5 & Bcftools 1.5
https://sourceforge.net/projects/samtools/files/samtools/1.5/
vcffilter #vcflib
https://github.com/vcflib/vcflib
ANNOVAR
http://annovar.openbioinformatics.org/en/latest/user-guide/download/
#annovar досутпен только последней версии
#следует проверить команду в скрипте, вдруг они добавили изменения
#у меня была устновлена версия 2016-02-01
на компьютере должна быть установлена java
В случае использование этих программ других версий
в некоторых случаях опции и даже команды могут не совпадать
(например, bcftools 0.1.19 и 1.5 вызвавают варианты разными командами - view и call соответсвенно)
если у вас установлены другие версии
проверьте правильность вызова команд, измените пути к программам
- Программы следует скачать и сложить в одну папку
я не буду заморачиваться и просить прописать путь для каждой из программ
на случай того, что они находятся в разных местах
это усложнит этот маленький скрипт и потребует много взаимодействия между
пользователем и программой
если вам тяжело сложить все программы в одну папку...
Если вы понимаете, что здесь нет ничего сложного,
и у вас программы разбросаны по папкам и часть из них добавлена в PATH:
измените скрипт - удалите пути для "PATH"-программ
допишите конкретные пути к другим
(.jar не может быть вызван при добавлении пути его расположения в PATH)
удалите уже ненужную перменную
будет проще пользоваться, да
- #for_newbie
Для запуска скрипта, стоит зайти в содержащую его папку через терминал
при помощи команды cd:
например, cd /media/illumina
#пусть он у вас лежит в этойпапке#
и затем вызвать его напрямую bash ./NGS_pipeline.sh
тут же его можно сделать исполняемым при помощи следующей команды:
chmod +x ./NGS_pipeline.sh
и запускать без слова bash -> ./NGS_pipeline.sh
также можно вызвать из любого места, указав путь сразу
bash /media/illumina/NGS_pipeline.sh
если скрипт уже исполняемый -> /media/illumina/NGS_pipeline.sh
При необходимости изменения опций программ и иного редактирования
можно открыть файл в любом текстовом редакторе
(желательно с подсветкой bash синтаксиса; тот же gedit на GNU/Linux).
- Я не обрабатывал исключения - это слишком долго и бессмысленно
весь stderr относительно запуска скрипта
будет постпупать в терминал
и вы можете сами догадаться, почитав
где вы забыли пробел, где путь не существует, где папка существует и т.д.
почти все используемые программы записывают файлы в stdout
а лог - в stderr; чтобы на засорять терминал и для его (лога) сохранности
stderr перенаправлен в лог-файлы
и вы можете всегда почитать, как был выполнен тот или иной этап
правильно/не_правильно и с статисткиой рассчетов
-
В начале алогритма есть маленькая команда
превращающая согласно регулярному выражению стандартные названия файлов
поступающие с прибора MiSeq Illumina в имена необходимой формы
если имена имеют другую стуктуру - в надлежащий вид их можно превратить
просто переименовав в ручную или изменив регулярное выражение внутри скрипта -
Коротко по программам:
А) Trimmomatic обрезает адаптеры на основании информации из файла TruSeq3-PE.fa
Остальные опции программы усредненнны.
Если использовалась другая пробоподготовка и другие адаптеры
можно изменить файл - найти более подходящий в папке Trimmomatic-0.36/adapters
но мне кажется, это не слишком существенно,
чтобы прописывать отдельно ввод файла с адаптерами, учитывая что не всегда очевидно, какой использовать
и огромной разницы при варьировании файлов получено не было
В) Новый референс - тогда с ним идет следующуая работа: samtools делает почти мгновенную индексацию, в то время как bowtie2-build долго строит индексы для выравнивания. Индексация референса для выравнивания занимает много времени и в этот момент онлайн-отгрузка информации в файлы может подвисать - может не отмечаться увеличения файлов-индексов в Krusader. Но, скорее всего, программа не зависла. У меня индексация рефернса hg19 размером в 3 Gb заняла 2 часа на одном ядре в 2 Gb. Подвисала раза 3. B процессе можно анализировать Log.txt.
C) Bowti2 выравнивает очищенные прочтения. Алгоритмы выравнивания парных и непарных прочтений разнятся между собой. Так как после trimmomatic получается и те, и другие, их следует выравнивать двумя разными вызовами bowtie2.
Весь вопрос состоит в том, насколько это целесообразно делать. Количество непарных прочтений обычно невелико. Парные и непарные прочтения выравниваются различными алгоритмами. Соответсвенно, при сливании получаемых из этих алгоритмов файлов в один может быть привнесена некоторая погрешность. Объединение файлов требуется для дальнейшего анализа.
Опять таки, вопрос в том, насколько вы хотите выжать максимум из ваших данных и насколько это соотносится с вносимой ошибкой. В случае чего, вы можете поробовать оба варианта
D) samtools конвертирует и вызывает статистику. Здесь происходит этап удаления дуплексов. Он необходим, но из-за необоснованных опасений потерять данные я не производил его раньше на своих данных. Предлагаю обработать данные с/без удаления, сравнить результаты по кол-ву и качеству выхода и решить для себя. Думаю в дальнейшем я буду удалять дуплексы. Но вам следует определиться самостоятельно.
E) Samtools собирает позиции, Bcftools вызывает варианты. В данном месте происходит фильтрация по глубине покрытия > 15, качества картирования > 30 и уверененности прибора > 40:
вот команда (строка 245)
vcffilter -f "DP > 15 & MQ > 30 & QUAL > 40 " ${i}.raw.vcf >
если хотите изменить параметр - откройте скрипт и измените цифры
также можно вызвать справку по vcffilter
(путем указания команды)
и узнать, по каким параметрам еще можно фильтровать и как удалять лишние
F) annovar аннотирует :)
Перед запуском скрипта НЕОБХОДИМО глянуть мануал ANNOVARa, открыть скрипт, сравнить команды. В моей версии ANNOVAR аннотирует при помощи скаченных баз данных - их нужно предварительно все скачать!! Нужно внимательно проверить опции в мануале и скрипте, чтобы не было никаких ляпов.
В команде после -protocol идут названия баз данных
после чего в -operation идут опции
кол-во баз данных должно равняться опциям!
скаченные базы данных следует распаковать в папку annovar/humandb/
команда для скачивания бд выглядит примерно так
(стоит уточнить на сайте; .pl - скрипт, лежащий в папке annovara)
annotate_variation.pl -buildver hg19 -downdb -webfrom annovar <база_данных> humandb/
выбрать базы данных можно на сайте
https://doc-openbio.readthedocs.io/projects/annovar/en/latest/user-guide/download/#-for-filter-based-annotation
если вы не знаете, какие скачать - скачивайте и распаковывайте те, которые находятся у нас в скрипте
refGene,esp6500siv2_all,avsnp147,clinvar_20170501,revel,intervar_20170202,1000g2015aug_all,1000g2015aug_afr,1000g2015aug_eas,1000g2015aug_eur,gnomad_genome