Bedtools

Для определения генов и глубины их покрытия было использовано выравнивание ридов с референсной последовательностью из прошлого практикума. В таблице ниже приведены использованные команды и комментарии к ним.

$ bedtools bamtobed -i sorted.bam > sorted.bed Переводит файл с ридами из формата .bam в формат .bed.
$ bedtools coverage -counts -a gencode.genes.bed -b sorted.bed | grep -wv 0 > cov.bed Определяет глубину покрытия генов ридами. Параметр -counts нужен для получения числа перекрываний без подсчета доли перекрытых нуклеотидов.
$ bedtools intersect -c -a gencode.genes.bed -b sorted.bed | grep -wv 0 > chr17_int.bed Находит пересечения ридов с последовательностью генома. Параметр -с позволяет получить число пересечений. С помощью grep избавляемся от пересечений с покрытием 0.

Больше всего ридов попало в ген KPNB1, ниже представлено его карткое описание.

Также ридами оказались покрыты два антисмысловых транскрипта - RP11-580I16.2 (33 вхождения) и RP11-138C9.1 (2718 вхождений)

ген KPNB1
Функция Транспортный фактор, участвующий в переносе белков в ядро
Размер, bp 35,656
Координаты и цепь chr17:47,649,838-47,685,493
Цепь Сomplement
Покрытие 165074
Экзоны/интроны 22/21

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

$ bedtools bamtofastq -i sorted.bam -fq sorted.fastq Переводит .bam файл с выравниваниями в формат .fastq.

-i имя входного .bed файла

-fq имя .fastq файла для записи

$ grep -w "KPNB1" gencode.genes.bed > kpnb1.bed

$ bedtools getfasta -fi chr17.fasta -bed kpnb1.bed -fo kpnb1.fasta

Генерирует файл с нуклеотидной последовательностью для гена KPNB1 - одного из покрытых нашими выравниваниями.

-fi название файла с нуклеотидной поседовательностью хромосомы

-fo позволяет указать название .fasta файла для записи (без этого параметра результат выводится в stdout)

$ infoseq -auto -nocolumns -delimiter '\t' -only -noheading -name -length chr17.fasta > chr17_len.txt

$ bedtools makewindows -g len.bed -w 1000000 > chr17_windows.bed

С помощью infoseq узнаем длину хромосомы - 81195210 bp, затем разбиваем нашу хроосому на фрагменты по 1 млн нуклеотидов. Всего было получено 82 таких интервала.

-g файл с названием хромосомы и ее длиной в bp

-w желаемая длина интервала