Отчет по практикуму 5. Ферменты и метаболические пути. База данных KEGG.

Отчет по практикуму 5. Ферменты и метаболические пути. База данных KEGG.

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

Мне была дана 4 хромосома человека chr4.fasta и файл с ридами для той же хромосомы chr4.fastq. Файл с хромосомой chr4.fasta и файл с ридами chr4.fastq был скопирован в директоию /nfs/srv/databases/ngs/tsvetcovroman. Анализ качества прочтений был проведён при помощи программы FastQC. Для улучшения качества использовался Trimmomatic. Программа FastQC была запущена через командную строку на kodomo: 
fastqc chr21.fastq
Вывод программы: 
 
Started analysis of chipseq_chunk4.fastq
Approx 10% complete for chipseq_chunk4.fastq
Approx 25% complete for chipseq_chunk4.fastq
Approx 40% complete for chipseq_chunk4.fastq
Approx 55% complete for chipseq_chunk4.fastq
Approx 65% complete for chipseq_chunk4.fastq
Approx 80% complete for chipseq_chunk4.fastq
Approx 95% complete for chipseq_chunk4.fastq
Analysis complete for chipseq_chunk4.fastq
В результате были получены 2 файла:chr21_fastqc.zip и chr21_fastqc.html.
Результаты в формате .html: chipseq_chunk4_fastqc.html.

Рис. 1. Информация о количестве ридов, откартированных на геном.


Рис. 2. Информация о количестве ридов, откартированных на геном.

Очистка чтений программой Trimmonatic не требуется.
Далее последовательности были откартированы на проиндескированный геном человека.
Использованная команда: 
bwa mem /srv/databases/ngs/hg19/GRCh37.p13.genome.fa chipseq_chunk4.fastq > chipseq_chunk4.sam
Выдача: 
[M::main_mem] read 7271 sequences (261756 bp)...
[M::mem_process_seqs] Processed 7271 reads in 1.996 CPU sec, 2.025 real sec
[main] Version: 0.7.10-r789
[main] CMD: bwa mem /srv/databases/ngs/hg19/GRCh37.p13.genome.fa chipseq_chunk4.fastq 
[main] Real time: 61.135 sec; CPU: 15.036 sec
После этого был проведен анализ полученного файла с помошью следующих команд:
Перевод в бинарный формат: samtools view -b chipseq_chunk4.sam -o chunk4.bam
Сортировка по координате в референсе начала чтения: samtools sort chunk4.bam -T chip_temp -o chipseq_chunk4.sorted.bam
Индексация: index chipseq_chunk4.sorted.bam
Получение информации о количестве откартированных ридов: samtools idxstats chipseq_chunk4.sorted.bam > chipseq_chunk4.idxstats
Просмотр количества откартированных ридов: samtools view -c chipseq_chunk4.sorted.bam
Получение количества откартированных ридов:samtools idxstats chipseq_chunk4.sorted.bam
Откартировался 7271 рид. 7235 ридов откартировалось на геном.
Подавляющее большинство ридов откартировались на chr12. Это позволяет предположить, что для анализа были предложены прочтения с 12 хромосомы. Полная информация о количестве откартированных ридов представлена в файле. Информация о количестве ридов, откартированных на геном представлена на Рисунке 1.

Рис. 1. Информация о количестве ридов, откартированных на геном.

Для поиска пиков была использована программа MACS. Попытка запуска командой: macs2 callpeak -t chipseq_chunk4.sorted.bam не дала результата, так как программа выдала 0 пиков. Далее я запустил программу со следущими параметрами: macs2 callpeak -n chunk4 -t chunk4_sorted.bam --nomodel.
Программа выдала следующие файлы:NA_peaks.narrowPeak , NA_peaks.xls, NA_summits.bed. Наиболее полная информация представлена в файле NA_peaks.xls. Программа нашла 8 пиков. Все пики находятся одном регионе 5 хромосомы.
В Таблице 1 представлена информация о пиках.
ИмяНачалоКонецШиринаabs_summitpileupfold_enrichment-log10pvalue-log10qvalue
NA_peak_193454544934547882459345464520.009.0517217.72276 10.87899
NA_peak_293481131934815364069348137728.007.5520818.99066 11.99228
NA_peak_393531850935321162679353198622.008.6466218.03990 11.17836
NA_peak_493543066935433142499354319029.0010.1351324.4788917.18385
NA_peak_593623268936236003339362342445.0013.2183940.6493331.17819
NA_peak_693855387938555872019385544721.009.0909118.28952 11.34986
NA_peak_793861258938615603039386136218.007.7235814.56367 8.03305
NA_peak_893963523939637492279396366722.006.3186813.81179 7.36060
Таблица. 1. Информация о пиках.


В Таблице 1 приведены показатели -log10pvalue и -log10qvalue. Чем выше данные параметры, тем ниже соответствующие показатели p-value и q-value, и следовательно, тем достовернее пик. Среди указанных в Таблице 1 самым достоверным является пик 5(-log10pvalue = 40.64933, -log10qvalue = 31.17819), а наименее достоверным - пик 8(-log10pvalue = 6.31868, -log10qvalue = 7.36060). Затем я поместил в начало файла NA_peaks.narrowPeak была помещена дополнительная информация так, как это описано по ссылке:http://genome.ucsc.edu/FAQ/FAQformat.html#format12.

Рис. 2. Визуализация пиков с помощью геномного браузера.


Рис. 3. Пик 5 в более крупном масштабе.

Ширина пика 5 составляет 311, достоверность пика можно оценить по p-value и q-value(-log10pvalue = 40.64933, -log10qvalue = 31.17819, следовательно p-value = 2.242*10-41, а q-value = 6.637*10-32). Вершина пика равноудалена от его начала и конца(на расстоянии 156 от начала при ширине пика 311). Начальная часть пика пересекается с повтором L1PA7.
Затем я визуализировал пик 4, который является вторым по достоверности(-log10pvalue = 24.47889, -log10qvalue = 17.18385).
Рис. 4. Пик 4 в более крупном масштабе.

Ширина пика 4 составляет 171, достоверность пика можно оценить по p-value и q-value(-log10pvalue = 24.47889, -log10qvalue = 17.18385, следовательно p-value = 3.319*10-25, а q-value = 6.546*10-18). Вершина пика смещена в сторону его конца(вершина расположена на расстоянии 124 от начала при ширине пика 171). Пик 4 не пересекается с функциональными элементами генома.

Ссылки:


[1] http://lib.stat.cmu.edu/S/bootstrap.funs
[2] Эфрон, 1979.