RNA-seq


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

    hisat2-build chr6.fna chr6 &> hisat2-build.log

Версия 2.1.0, log-файл

На вход была подана фаста-последовательност хромосомы. В итоге - восемь файлов с нумерацией от 1 до 8. Они, увы, бинарные и в документации не сказано точно, что в них именно, но, вероятно, это индексы k-меров длин 1, 2, 3, 4, ... соответственно


2. Оценка качества чтений

Программа fastqc не принимает на вход архив .fastq.gz, поэтому его необходимо предварительно переконвертировать:

    zcat ../SRR2015720_1.fastq.gz 2> zcat.log | fastqc stdin --extract &> fastqc.log

Версия zcat 1.9, log-файл zcat(пуст); версия fastqc v0.11.8, log-файл fastqc

Всего чтений 30 271 388.

Per base quality
Качество чтений не очень хорошее, но приемлимое. Практчески не выходит из зелёной зоны.
Sequence length distribution
Все риды имеют одинаковую длину 101.

3. Картирование чтений на референс

    hisat2 -x ../hisat2/chr6 -k 3 -U SRR2015720_1.fastq.gz
                -S mapped/SRR2015720_1.sam --threads 10 2> mapped/hisat2.log

Версия 2.1.0, log-файл

Картированных очень мало, но это логично: ведь в образце у нас вся мРНК, а в качестве референса только одна хромосома.


    samtools sort -o SRR2015720_1.bam SRR2015720_1.sam &> samtools_sort.log

Версия samtools здесь и далее 1.9 (htslib 1.9), log-файл


    samtools index SRR2015720_1.bam &> samtools_index.log

log-файл (пуст)


    samtools sort -o SRR2015720_1.bam SRR2015720_1.sam &> samtools_sort.log

log-файл(пуст)


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

4.1 Устройство gtf файла

Устройство файла:
В шапке файла описание и контакты.
В теле: название хромосомы,название источника, тип куска (ген, транскрипт, экзон, интрон, нетранслируемый участок и т.д.), координаты конца и начала, score(??), цепь (+ или -), свойства: id, название и тип гена и прочее.

    mawk '$3="gene"' gencode.chr6.gtf | wc -l

Всего на шестой хромосоме размечен 87 762 ген.


4.2 Анализ картирования

    htseq-count -f bam -s no -t gene -m union ../mapped/SRR2015720_chr6.bam gencode.chr6.gtf
                                                        > htseq-count.results 2> htseq-count.log

Версия 0.11.2, log-файл, результат

Параметры:
-f = тип файла с выравниванием, в нашем случае "bam"
-s = информация о цепи. Мы не знаем, поэтому "no".
-t = тип, который будет рассматривать, остальные проигнорирует. В данном случае нас интересуют все гены.
-m = что делать с теми чтениями, которые попали на фрагменты с несколькими свойствами сразу. В данном случае "union".

Результат:

Всего 1 303 216
Мимо границ генов 152 103 (12%)
На границы генов 248 929 (19%)