Практикум 14. Анализ данных NGS: RNA-seq
В ходе данного практикума будут пройдены основные этапы обработки данных RNA-seq: подготовка референсной последовательности, оценка качества прочтений и их картирование на референс, поиск генов по разметке.
Подготовка референса
В качестве референсной последовательности снова использовалась последовательность четвёртой хромосомы человека. Референс был проиндексирован при помощи программы hisat2:
hisat2-build chr4.fna chr4 &> hisat2_run.log
На вход алгоритма был подан файл с референсной последовательностью, на выходе было получено 8 бинарных файлов chr4.?.ht2 (? = 1, ..., 8). С логами выполнения команды можно ознакомиться здесь.
Оценка качества чтений
Качество прочтений было проанализировано с использованием программы fastqc:
fastqc SRR2015720_1.fastq.gz &> fastqc_run.log
Далее приведены ссылки на лог-файл программы и на html-страницу SRR2015720_1_fastqc, полученную на выходе.
Общее число прочтений получилось равным 30 271 388.
На рисунке 1 представлен график зависимости качества прочтения нуклеотида от его позиции в чтении.
Можно видеть, что качество прочтений в целом достаточно хорошее (минимальное значение лишь для двух позиций опускается ниже 28), однако на концевых участках наблюдается достаточно значительное снижение качества по сравнению с центральными участками.
На рисунке 2 приведён график распределения длины прочтений. Очевидно, что 100% прочтений имеет длину 101 нуклеотид.
Картирование чтений на референс
Для этого также использовалась программа hisat2:
hisat2 -p 10 -x ../../reference_hisat2/chr4 -k 3 -U SRR2015720_1.fastq.gz > mapped_reads.sam 2> hisat2_map_run.log
Здесь hisat2 - вызываемая программа, -p 10 - опция, указывающая использовать параллельно 10 потоков,
Затем полученный sam-файл был переведён в сортированный bam-файл, который следом был проиндексирован:
samtools sort -o mapped_reads.bam mapped_reads.sam
samtools index mapped_reads.bam
Далее были отобраны прочтения, которые при картировании легли на четвёртую хромосому:
samtools view -h -bS mapped_reads.bam NC_000004.12 > chr4_mapped_reads.bam
Поиск экспрессирующихся генов
Для этой задачи потребуется gtf-файл с разметкой генов. Далее приведены первые несколько строк файла gencode.chr4.gtf:
##description: evidence-based annotation of the human genome (GRCh38), version 35 (Ensembl 101) ##provider: GENCODE ##contact: gencode-help@ebi.ac.uk ##format: gtf ##date: 2020-06-03 NC_000004.12 HAVANA gene 49096 50124 . - . gene_id "ENSG00000248302.3"; gene_type "transcribed_processed_pseudogene"; gene_name "BNIP3P41"; level 1; hgnc_id "HGNC:49721"; tag "pseudo_consens"; havana_gene "OTTHUMG00000159940.3"; NC_000004.12 HAVANA transcript 49096 49956 . - . gene_id "ENSG00000248302.3"; transcript_id "ENST00000616656.1"; gene_type "transcribed_processed_pseudogene"; gene_name "BNIP3P41"; transcript_type "processed_transcript"; transcript_name "BNIP3P41-202"; level 2; transcript_support_level "NA"; hgnc_id "HGNC:49721"; tag "basic"; havana_gene "OTTHUMG00000159940.3"; havana_transcript "OTTHUMT00000476945.1"; NC_000004.12 HAVANA exon 49096 49956 . - . gene_id "ENSG00000248302.3"; transcript_id "ENST00000616656.1"; gene_type "transcribed_processed_pseudogene"; gene_name "BNIP3P41"; transcript_type "processed_transcript"; transcript_name "BNIP3P41-202"; exon_number 1; exon_id "ENSE00003749735.1"; level 2; transcript_support_level "NA"; hgnc_id "HGNC:49721"; tag "basic"; havana_gene "OTTHUMG00000159940.3"; havana_transcript "OTTHUMT00000476945.1"; NC_000004.12 HAVANA transcript 49554 50124 . - . gene_id "ENSG00000248302.3"; transcript_id "ENST00000503774.1"; gene_type "transcribed_processed_pseudogene"; gene_name "BNIP3P41"; transcript_type "transcribed_processed_pseudogene"; transcript_name "BNIP3P41-201"; level 1; transcript_support_level "NA"; hgnc_id "HGNC:49721"; ont "PGO:0000004"; ont "PGO:0000019"; tag "pseudo_consens"; tag "basic"; havana_gene "OTTHUMG00000159940.3"; havana_transcript "OTTHUMT00000358399.2"; NC_000004.12 HAVANA exon 49554 50124 . - . gene_id "ENSG00000248302.3"; transcript_id "ENST00000503774.1"; gene_type "transcribed_processed_pseudogene"; gene_name "BNIP3P41"; transcript_type "transcribed_processed_pseudogene"; transcript_name "BNIP3P41-201"; exon_number 1; exon_id "ENSE00002063424.1"; level 1; transcript_support_level "NA"; hgnc_id "HGNC:49721"; ont "PGO:0000004"; ont "PGO:0000019"; tag "pseudo_consens"; tag "basic"; havana_gene "OTTHUMG00000159940.3"; havana_transcript "OTTHUMT00000358399.2"; NC_000004.12 HAVANA gene 53286 88208 . + . gene_id "ENSG00000272602.6"; gene_type "protein_coding"; gene_name "ZNF595"; level 2; hgnc_id "HGNC:27196"; havana_gene "OTTHUMG00000159865.6"; NC_000004.12 HAVANA transcript 53286 87523 . + . gene_id "ENSG00000272602.6"; transcript_id "ENST00000608255.2"; gene_type "protein_coding"; gene_name "ZNF595"; transcript_type "protein_coding"; transcript_name "ZNF595-203"; level 2; protein_id "ENSP00000476367.2"; transcript_support_level "2"; hgnc_id "HGNC:27196"; tag "basic"; tag "CCDS"; ccdsid "CCDS75077.1"; havana_gene "OTTHUMG00000159865.6"; havana_transcript "OTTHUMT00000357816.4"; NC_000004.12 HAVANA exon 53286 53491 . + . gene_id "ENSG00000272602.6"; transcript_id "ENST00000608255.2"; gene_type "protein_coding"; gene_name "ZNF595"; transcript_type "protein_coding"; transcript_name "ZNF595-203"; exon_number 1; exon_id "ENSE00003752731.1"; level 2; protein_id "ENSP00000476367.2"; transcript_support_level "2"; hgnc_id "HGNC:27196"; tag "basic"; tag "CCDS"; ccdsid "CCDS75077.1"; havana_gene "OTTHUMG00000159865.6"; havana_transcript "OTTHUMT00000357816.4"; NC_000004.12 HAVANA exon 85731 87523 . + . gene_id "ENSG00000272602.6"; transcript_id "ENST00000608255.2"; gene_type "protein_coding"; gene_name "ZNF595"; transcript_type "protein_coding"; transcript_name "ZNF595-203"; exon_number 2; exon_id "ENSE00003704644.2"; level 2; protein_id "ENSP00000476367.2"; transcript_support_level "2"; hgnc_id "HGNC:27196"; tag "basic"; tag "CCDS"; ccdsid "CCDS75077.1"; havana_gene "OTTHUMG00000159865.6"; havana_transcript "OTTHUMT00000357816.4"; NC_000004.12 HAVANA CDS 86054 87448 . + 0 gene_id "ENSG00000272602.6"; transcript_id "ENST00000608255.2"; gene_type "protein_coding"; gene_name "ZNF595"; transcript_type "protein_coding"; transcript_name "ZNF595-203"; exon_number 2; exon_id "ENSE00003704644.2"; level 2; protein_id "ENSP00000476367.2"; transcript_support_level "2"; hgnc_id "HGNC:27196"; tag "basic"; tag "CCDS"; ccdsid "CCDS75077.1"; havana_gene "OTTHUMG00000159865.6"; havana_transcript "OTTHUMT00000357816.4";
Файл gtf состоит из строк заголовка, начинающихся с «##», и расположенных ниже записей в 9 полях:
• sequence - название последовательности
• source - источник аннотации локальных особенностей (features)
• feature - тип особенности (ген, экзон и пр.)
• start - координата начала особенности
• end - координата конца особенности
• score - значение, характеризующее «уверенность» источника в аннотации
• strand - направление цепи, на которой расположена особенность
• phase - число, указывающее сдвиг рамки считывания относительно координаты начала (для кодирующих последовательностей)
• attributes - дополнительная информация
Далее при помощи программы htseq-count был проведён подсчёт числа прочтений, картированных на каждый из генов в разметке:
htseq-count -f bam -s no -m intersection-nonempty -t exon chr4_mapped_reads.bam gencode.chr4.gtf > exon_reads.txt 2> htseq_run.log
Здесь htseq-count - используемая пргорамма, -f bam - формат файла с прочтениями, поданного на вход, -s no - опция, указывающая, что прочтение считается перекрывающимся с размеченной особенностью вне зависимости от направления цепи, -m intersection-nonempty - опция, задающая режим обработки прочтений, перекрывающихся с более чем одной особенностью, -t exon - опция, указывающая тип особенностей, пересечения с которыми надо искать. Логи выполнения команды можно просмотреть здесь.
В результате получилось, что 230 866 прочтений не попали в границы генов (информация в конце файла). Число прочтений в границах генов было подсчитано при помощи приведённой ниже команды и составило 717 938.
awk 'match($1,"_")!=1{s+=$2}END{print s}' exon_reads.txt