Практикум 13. Работа с bedtools.

Задание 1. Глубина покрытия генов чтениями

Входной файл	Выходной файл	Программа	Комментарий
chr11_align.bam	chr11_align.bed	/P/y14/term3/block4/SNP/bedtools2/bin/bedtools bamtobed -i chr11_align.bam > ../pr13/chr11_align.bed	[программа была запущена из папки pr12] Меняем формат выравнивания на удобный для пакета bedtools. Заодно телепортируемся в отдельную директорию ;)
/P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed; chr11_align.bed	intersect.bed	/P/y14/term3/block4/SNP/bedtools2/bin/bedtools intersect -c -a /P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed -b chr11_align.bed > intersect.bed	Параметр -c помогает считать количество ридов, перекрывающихся с feature (с геном):"-c For each entry in A, report the number of overlaps with B." В качестве файла А берем наш референс, то есть разметку генов. В качестве файла B берем файл с координатами ридов
intersect.bed	intersect_chr11_nonull.bed	grep -r "^chr11" intersect.bed > intersect_chr11.bed \| grep -v "0$" intersect_chr11.bed > intersect_chr11_nonull.bed	Выбираем из полученного файла то, что относится к хромосоме 11, а также избавляемся от строк с генами, на которые попало 0 ридов. В принципе достаточно было бы только второго grep-a, так как пересечение с другими хромосомами не может быть ненулевым. Но пересечений с 11 хромосомой не нашлось тоже! Попробуем по-другому.
/P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed; chr11_align.bed	intersect_nonull_u.bed	/P/y14/term3/block4/SNP/bedtools2/bin/bedtools intersect -u -a /P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed -b chr11_align.bed > intersect_nonull_u.bed	Параметр -u выводит только ненулевые пересечения, хотя и не показывает количество ридов, упавших на ген. Файл ненулевой! Поищем по названию генов, которые он содержит, в файле intersect_chr11.bed... Оказывается, и там есть ненулевые наложения. Значит, проблема в grep-e. Ну во-первых, строка может кончаться на ноль, даже если кончается не на цифру ноль, а на цифру, заканчивающуюся на ноль. Во-вторых, кавычки там, кажется, ни к чему...
intersect_chr11.bed	intersect_chr11_nonull_1.bed	grep -w -v 0 intersect_chr11.bed > intersect_chr11_nonull_1.bed	Параметр -w позволяяет находить только строки, содержащие "0" как отдельное слово. Ура! Получили желаемый результат. Покрытыми оказались 4 гена: HSPA8 (heat shock protein 8) и SNORD14E, SNORD14C, SNORD14D (U14 small nucleolar RNA - гены малых ядрышковых РНК). По-видимому, наша разметка генов gencode.genes.bed содержала несколько координат для HSPA8, потому что в итоговом файле ему соответствовало несколько строк с отличающимися координатами, в связи с чем сложно посчитать, сколько всего ридов на него легло. Самое большое значение покрытия - 38733. На SNORD14E легло 561 ридов, на SNORD14D - 444 ридов, на SNORD14C - 254 ридов.

Дополнительные задачи

1. Получите из файла в выравниванием файл с чтениями в формате fastq
исходный файл: ../pr12/chr11_align.bam
Команда:

/P/y14/term3/block4/SNP/bedtools2/bin/bedtools bamtofastq -i ../pr12/chr11_align.bam -fq ctr11_align.fastq

результат: ctr11_align.fastq

7. Получите координаты одного из покрытых Вашими чтениями генов, расширенные на 1000 нуклеотидов в обе стороны.
исходный файл: intersect_chr11_nonull_1.bed
Команда:

head -1 intersect_chr11_nonull_1.bed > line1.bed | /P/y14/term3/block4/SNP/bedtools2/bin/bedtools slop -b 1000 -i line1.bed -g /P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed > line1_sloped.bed

результат: line1_sloped.bed
Точно сработано!

8. Получите координаты одного из покрытых Вашими чтениями генов, сдвинутые на 500 нуклеотидов ближе к началу хромосомы
исходный файл: intersect_chr11_nonull_1.bed
Команда:

/P/y14/term3/block4/SNP/bedtools2/bin/bedtools shift -s -500 -i line1.bed -g /P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed > line1_shifted.bed

результат: line1_shifted.bed
На самом деле начало сдвинулось на 501, а конец на 502. Скорее всего это связано с какими-то особенностями в записи генома.

9. Получите непересекающиеся фрагменты, соответствующие области, покрытой Вашими чтениями
исходный файл: intersect_chr11_nonull_1.bed
Команда:

/P/y14/term3/block4/SNP/bedtools2/bin/bedtools merge -i intersect_chr11_nonull_1.bed > chr11_intersect_merge.bed
Error: Sorted input specified, but the file intersect_chr11_nonull_1.bed has the following out of order record
chr11	122931828	122932037	protein_coding	HSPA8	9316

По-видимому, какая-то ошибка в программе: указанная в error строка соответствует формату, как и другие строки.

P.S. Все файлы, полученные в ходе выполнения практикума, можно найти в папке srv/databases/ngs/anrozina/pr13.