Практикум 12

Задача: картировать чтения, полученные в результате секвенирования транскриптома (человек, версия сборки генома hg19)

Часть I: подготовка чтений

1. Чтения, картирующиеся на участок хромосомы человека (получены путем секвенирования РНК). Файлы с одноконцевыми чтениями в формате fastq лежат на kodomo в директории /P/y14/term3/block4/SNP/rnaseq_reads. Распределение файлов по студентам см. в табл.. Берем ту же хромосому, что и была при выполнении заданий практикума 11 (!) Для каждой хромосомы есть 2 (!) файла fastq - это 2 биологические реплики. В рамках обязательного задания достаточно работать с одной (любой) репликой.

2. Разметка человеческого генома по версии Gencode19 для сборки hg19. Разметка в формате .gtf лежит на kodomo в директории /P/y14/term3/block4/SNP/rnaseq_reads (gff и gtf формат - это почти одно и то же)

1. Анализ качества чтений.

Сделайте контроль качества Ваших чтений с помощью программы FastQC.
Комментарий: программа FastQC установлена на kodomo, её можно вызвать командой "fastqc chr21.1.fastq".
Версию с графическим интерфейсом можно поставить на свой компьютер.

В результате работы программы Вы получите архив (.zip), который содержит отчет о программе в виде html файла.

Команда: fastqc chr21.1.fastq

Ссылка на отчет fastqc

Рис.1. "Per base quality"

2. Картирование чтений.

Откартируйте очищенные чтения с помощью программы Hisat2.
Этапы
Сначала необходимо проиндексировать референсную последовательность; команда hisat2-build

Команда:export PATH=${PATH}:/home/students/y06/anastaisha_w/hisat2-2.0.5
Команда:hisat2-build chr21.fasta chr21

В результате возникли файлы chr21.1.ht2,...,chr21.6.ht2

Затем построить выравнивание прочтений и референса в формате .sam. Запустите hisat2 с параметром --no-softclip. Параметр "--no-spliced-alignment" так как при сплайсинге перестановки могут быть важны.

Команда: hisat2 -x chr21 -U chr21afterafter.fastq --no-softclip chr21.sam

3. Анализ выравнивания

Переведите выравнивание чтений с референсом в бинарный формат .bam. Используйте пакет samtools, команда view: samtools view;

Команда: samtools view -b chr21.1.sam > chr21.1.bam Отсортируйте выравнивание чтений с референсом (получившийся после картирования .bam файл) по координате в референсе начала чтения; команда samtools sort;

Команда: samtools sort chr21.1.bam chr21sort.1
В итоге имеем файл chr21sort.1.bam

Проиндексируйте отсортированный .bam файл командой samtools index

Команда: samtools index chr21sort.1.bam
Появился новый файл chr21sort.1.bam.bai

Рис.2. Информация из Hisat2

4. Подсчет чтений

Файл с разметкой генов по версии Gencode для генома человека сборки hg19 лежит тут: /P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed
Программа bedtools лежит тут: /P/y14/term3/block4/SNP/bedtools2/bin
Пример запуска: /P/y14/term3/block4/SNP/bedtools2/bin/bedtools intersect -h
Подсказка: предлагается из файла с выравниваем достать координаты каждого выровненного чтения (.bed), а затем пересечь эти координаты с разметкой генов.
Подумайте, в каком порядке пересекать файлы и какую опцию стоит добавить, чтобы сразу посчитать количество чтений.

Первым делом нужно перевести файл из формата .bam в формат .bed

Команда:/P/y14/term3/block4/SNP/bedtools2/bin/bedtools bamtobed -i chr21sort.1.bam > chr21.1.bed

Далее воспользуемся командой intersect

Команда: /P/y14/term3/block4/SNP/bedtools2/bin/bedtools intersect -a chr21sort.1.bed -b /P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed

Таким образом были найдены 2 гена: Usp16 и Cct8.

Usp16:
Ubiquitin carboxyl-terminal hydrolase 16 is an enzyme that in humans is encoded by the USP16 gene.
Координаты: 30,396,950-30,426,809
Длина:29,859
Экзоны:20
Чтений на ген:14

Ссt8:
T-complex protein 1 subunit theta is a protein that in humans is encoded by the CCT8 gene.
Координаты: 30,428,126-30,446,118
Длина: 17,992
Экзоны:16
Чтений на ген:188

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