Bedtools

Задание 0. Гены из pr13. Покрытие и особенности

Input Output Программа Что делает
sortedbam.bam aln.bed 
bedtools bamtobed -i sortedbam.bam > aln.bed
 Переводим файл с выравниванием в формате .bam в формат, удобный
для работы bedtools - .bed
aln.bed sorted.int.chr8.bed 
bedtools intersect -a /P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed -b aln.bed -c >
 intersect.bed ; grep -r "^chr8"  intersect.bed  | grep -r -v "0$" > sorted.int.chr8.bed
 Получаем пересечение наших чтений с последовательностью генома (по координатам), для каждой хромосомы;
далее отбираем строки с интересующими нас данными - chr8, с покрытием более 0
в программу закралась ошибка и, скорее всего, позже она будет исправлена
aln.bed mb.genes.bed 
bedtools intersect -a /P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed -b aln.bed -u  > mb.genes.bed
Практически то же самое, но число находок больше в сравнении с предыдущим способом; не выводится глубина покрытия
aln.bed wawb.bed 
bedtools intersect -a /P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed -b aln.bed -wa -wb > wawb.bed
 Каждое перекрывание выписывается в output одной строкой;
gencode.genes.bed, aln.bed cover.bed 
bedtools coverage -a /P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed -b aln.bed |
 grep -r "^chr8" > cover.bed
 Определяем покрытие разметки нашими ридами

Комментарии к таблице

$ bedtools bamtobed -i sortedbam.bam > aln.bed
-i Указываем имя файла в формате .bam, который хотим перевести в формат .bed 

$ bedtools intersect -a /P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed -b aln.bed -c -v ...
-a имя файла с разметкой генов для генома человека сборки hg19 в формате .bed -b имя файла с выравниванием; из него будут браться координаты для
пересечения с координатами разметки генов
-c флажок, который дает для каждого элемента (в нашем случае для каждой хромосомы) число
пересечений с чтениями из файла с выравниваниями 

$ bedtools intersect -a /P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed -b aln.bed -u  > mb.genes.bed
-u флажок позволяет выписывать данные о пересечении в том случае, если оно является ненулевым 

bedtools intersect -a /P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed -b aln.bed -wa -wb > wawb.bed
-wa, -wb формат вывода, когда строка output-файла соответствует одному покрытию ридом

Ниже приведена таблица с основной информацией по найденным генам
С помощью команды bedtools coverage было определено покрытие последовательности хромосомы
нашими ридами (данные по ней были "вырезаны" с помощью grep из общего файла с результатами
по всему геному
Можно заметить, что файл содержит много дублирующихся строк, перекрывающихся интервалов, что значительно снижает
число покрытия. Для нашей работы просто удалим строки-дубликаты (результаты представлены в таблице ниже)

Поиск информации о количестве экзонов и полном количестве гена производился с помощью NCBI, Ensembl
Для AC103686.1 и Y_RNA в базах данных не было найдено информации, которую можно было бы привести здесь даннные по
интересующим нас вопросам были найдены только для белок-кодирующих генов)

Ген Цепь Положение в геноме Продукт Покрытие (без удаления дубликатов) Покрытие (с удалением дубликатов Экзоны Полное название Функция
AC103686.1 + miRNA 29 29 Под этим идентификатором в NCBI мы получаем chr8 полностью
(complete sequence); позможно, этот участок не является геном
COL22A1 - chr8:8q24.23-q24.3 protein_coding 34 18 71 collagen type XXII alpha 1 chain Белок способствует стабилизации миотендиновых соединений, прикреплению
скелетных мышц во время сократительной активности
MCM4 - 8q11.21 protein_coding 4031 1311 16 minichromosome maintenance complex component 4 Белок, кодируемый этим геном, является одним из высококонсервативных небольших хромосомных
поддерживающих белков, которые необходимы для инициирования репликации эукариот
PRKDC + 8q11.21 protein_coding 217954 59912 86 protein kinase, DNA-activated, catalytic subunit Этот ген кодирует каталитическую субъединицу ДНК-зависимой протеинкиназы
Y_RNA + misc_RNA 114 114 На NCBI и Ensembl записей не найдено; на Ensembl нечто, наиболее подходящее по
координатам это compmerge.1114.K562.chr8, но он не соответствует нашему варианту по
направлению цепи

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

Задание 1 bam.bam aln.fq 
bedtools bamtofastq -i bam.bam -fq aln.fq
 Переводим выравнивание из формата .bam в формат .fastq
Задание 2 chr8.fasta, coord.bed COL22A1.fasta 
bedtools getfasta -fi ../chr8.fasta -bed coord.bed > COL22A1.fasta
 Получаем файл в формате .fasta для гена COL22A1; в файле chr8.fasta содержится последовательность
chr8, в файле coord.bed - координаты гена, которые были проверены в базе данных Ensembl
Задание 4 aln.bed clusters.bed 
bedtools cluster -i aln.bed -s -d 30 > clusters.bed
 Ксластеризуем наши риды; объединяем только по одноименным цепочкам и с максимальным расстоянием
в 30 bp; всего получилось 15 кластеров

Комментарии к заданиям


Задание 1 
$ bedtools bamtofastq -i bam.bam -fq aln.fq
-i указываем имя файла в формате .bam, над которым мы хотим работать
-fq указываем имя файла, который будет получен на выходе; разумеется, в формате .fastq

Задание 2 
$ bedtools getfasta -fi ../chr8.fasta -bed coord.bed > COL22A1.fasta
-fi указываем имя файла в формате .fasta, над которым мы хотим работать
-bed указываем имя файла, в котором лежат координаты интересующего нас участка в формате .bed

Задание 4 
$ bedtools cluster -i aln.bed -s -d 30 > clusters.bed
-i указываем имя файла в формате .bed, над которым мы хотим работать
(выравнивание из обязательного задания (см. Задание 0 ранее)
-s флажок, который помогает проводить кластеризацию по одноименным цепям
-d после флажка указываем целое число - максимальное расстояние между ридами,
которые должны будут попасть в один кластер. В нашем случае оно было принято за 30

© Беляева Юлия, 2018