Анализ транскриптомов. Bedtools.
Картирование чтения транскриптома на референсный геном. Работа с 7 хромосомой, что и в предыдущем практикуме.
| fastqc chr7.1.fastq | Проверка качества чтений |
|---|---|
| java -jar /nfs/srv/databases/ngs/suvorova/trimmomatic/trimmomatic-0.30.jar SE -phred33 chr7.1.fastq chr7_trimmed.fastq TRAILING:20 MINLEN:50 | Тримминг |
| fastqc chr7_trimmed.fastq | Проверка качества после тримминга |
| hisat2 -x chr7 -U chr7_trimmed.fastq -S chr7_align.sam --no-softclip | Выравнивание прочтений и референса |
| samtools view -b chr7_align.sam -o chr7_align.bam | Перевод выравнивания в бинарный формат |
| samtools sort chr7_align.bam chr7_align_sorted | Сортировка выравнивания по координате в референсе |
| samtools index chr7_align_sorted.bam | Индексирование отсортированного файла |
| htseq-count -f bam chr7_align.bam -i gene_id -s no gencode.v19.chr_patch_hapl_scaff.annotation.gtf > counts.txt | Подсчет чтений |
| grep -wv 0 counts.txt > counts_notzero.txt | Создание файла с числом чтений только для генов, на которые чтения в принципе были картированы. |
Оценка качества чтений и очистка мусора
Исходно было получено 15940 чтения. Качество чтений до очистки отображено на Рис. 1. С помощью Trimmomatic с конца каждого чтения были отрезаны нуклеотиды с качеством ниже 20, после чего оставлены только чтения длиной не меньше 50 нуклеотидов. После очистки осталось 15856 ридов. Качество чтений после очистки отображено на Рис. 2. В целом хочется сказать, что качество чтений по сравнению с предыдущим практикумом гораздо выше.
Картирование ридов
В результате выравнивания на последовательность хромосомы легли один раз 99.14% (15720) ридов. Больше одного раза ни один рид не лег. 99.14% картировавшихся ридов это хороший показатель.
Также при запуске программы hisat2 был убран параметр --no-spliced-alignment, так как РНК может сплайсироваться.
Подсчет чтений
Для подсчета чтений была испольована программа htseq-count, имеющая следующие параметры:
-f формат входного файла.
-s выбор цепи
-i GFF атрибут, который используется в качестве feature ID (по умолчанию gene_id, подходящий для Ensembl GTF файлов).
-m что делать с ридами, попадающими на несколько генов.
Результаты
Итоговый файл содержал следующие строчки. Это те места, куда картировались наши риды.
ENSG00000106034.13 6205 ENSG00000207090.1 4 ENSG00000212628.1 5 __no_feature 9309 __ambiguous 197 __not_aligned 136
Как можно видеть, большинство ридов(из осмысленных) картировалось на ген кадгерин-подобного белка ENSG00000106034.13. Остальные два гена это малая ядерная и рибосомальная РНК. Скорее всего это издержки пробоподготовки. Также стоит отметить что огромное количество ридов оказалось откартировано непонятно куда - feature у них нет.