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

Введение в транскриптомный анализ

Подготовка референса

Для начала была создана директория hisat2, куда была скопирована референсная последовательность десятой хромосомы человека. Затем референс был проиндексирован следующим образом: $ hisat2-build -p 10 chr10.fna chr10 На вход мы подали последовательность десятой хромосомы человека, а на выходе получилось 8 бинарных файлов, которые называются chr10.N.ht2, где N - цифра от 1 до 8.

Проверка качества чтений

Использовать будем снова fastqc. Применим такую команду: fastqc SRR2015719_1.fastq.gz Был получен файл с данными, по которым можно оценить качество чтений. Всего в файле содержится 34642046 чтения. Понуклеотидное качество чтений можно увидеть на Рисунке 1. Видно, что абсолютное большинство нуклеотидов в чтениях находятся в "зелёной" зоне, т.е. они очень качественные и им можно доверять. Распределение по длинам чтений можно видеть на Рисунке 2. Все чтения имеют длину 101 нуклеотид.

Картирование чтений на геном

Чтобы картировать геном, будем использовать следующую команду: $ hisat2 -p 10 -x ../hisat2/chr10 -k 3 -U SRR2015719_1.fastq.gz > rna_reads.sam 2> map_log.txt Посмотрим на файл с логами:

$ cat map_log.txt 
34642046 reads; of these:
  34642046 (100.00%) were unpaired; of these:
    33242050 (95.96%) aligned 0 times
    1359657 (3.92%) aligned exactly 1 time
    40339 (0.12%) aligned >1 times
4.04% overall alignment rate

Видим, что всего картировалось 4.04% чтений, а больше одного раза картировалось всего 0.12% чтений от общего числа чтений.
Теперь конвертируем sam-файл в bam-файл и последний проиндексируем. Сразу отберём такие чтения, которые картировались на нашу хромосому. $ samtools sort -o rna_reads.bam rna_reads.sam
$ samtools index rna_reads.bam
$ samtools view -h rna_reads.bam NC_000010.11 | samtools sort -o rna_chr10.bam

Поиск экспрессирующихся генов

Файл gencode.chr10.gtf состоит из заголовка, где написано, что это вообще за файл, а также строк с разметкой. В каждой строке есть информация с координатами и о том, что по этим координатам находится.
Чтобы посчитать количество чтений, попавших на каждый из генов, применим следующую команду: $ htseq-count -f bam -s no -t exon -m union -o count.sam rna_chr10.bam gencode.chr10.gtf 1> htseq_log1.txt 2> htseq_log2.txt htseq-count – программа, которая считает.
-f bam – эта опция говорит, исходный формат файла с картированными чтениями в формате bam.
-s no – эта опция говорит, что чтения могут попадать как на прямую, так и на обратную цепь.
-t exon – эта опция говорит считать только те чтения, которые попали именно в границы экзонов.
-m union – эта опция говорит объединять чтения, если они перекрылись.
-o count.sam – это выходной файл, куда опять запишутся картированные чтения
rna_chr10.bam – это исходный файл.
gencode.chr10.gtf – это файл с разметкой генома.
1> htseq_log1.txt – это первый файл с логами.
2> htseq_log2.txt – это второй файл с логами.
Информация по генам содержится в первом файле (htseq_log1.txt). Итак, мимо границ чтений (__no_feature) попало 315456 чтения. Чтобы посчитать количество чтений, попавших в границы генов, я написал простенький python-скрипт:

with open('htseq_log1.txt', 'r') as file:
    count = 0
    for line in file:
        if line.startswith('__'):
            break
        count += int(line.strip().split('\t')[1])
    print(count)

Всего их получилось 999631.