Этот пример показывает, как выполнить общий для генома анализ метилирования ДНК у человека с помощью секвенирования генома.
Примечание.Для повышения эффективности MathWorks рекомендует запускать этот пример на 64-разрядной платформе, поскольку объем памяти близок к 2 ГБ. На 32-разрядной платформе, если вы получаете "Out of memory" ошибки при запуске этого примера, попробуйте увеличить виртуальную память (или поменять пространство) вашей операционной системы или попробуйте задать 3GB switch (только 32-разрядная версия Windows ® XP). Эти методы описаны в этом документе.
Метилирование ДНК является эпигенетической модификацией, которая модулирует экспрессию генов и поддержание геномной организации в нормальных и болезненных процессах. Метилирование ДНК может определять различные состояния камеры, и это наследуемо во время репликации камеры. Аберрантные шаблоны метилирования ДНК были связаны с генами раковых и опухолевых супрессоров.
В этом примере вы исследуете профили метилирования ДНК двух раковых камер человека: родительских клеток HCT116 камер рака толстой кишки и DICERex5 камер. DICERex5 камер получают из HCT116 камер после усечения DICER1 аллелей. Serre et al. в [1] предложен для изучения профилей метилирования ДНК с использованием MBD2 белка в качестве метил CpG-связывающего области и затем использовано высокопроизводительное секвенирование (HTseq). Этот метод обычно известен как MBD-Seq. Краткие чтения для двух повторений двух выборок были отправлены в архив SRA NCBI авторами [1]. Существуют и другие технологии для изучения статуса метилирования ДНК сайтов CpG в сочетании с HTseq, например MeDIP-seq или использования рестрикционных ферментов. Можно также анализировать наборы данных этого типа, следуя подходу, представленному в этом примере.
Можно получить несопоставленные одноконцевые чтения для четырех экспериментов по секвенированию из NCBI. Короткие показания были получены с использованием анализатора II генома Illumina ®. Средний размер вставки составляет 120 п.о., а длина коротких показаний составляет 36 п.о.
Этот пример принимает, что вы:
(1) загрузил файлы SRR030222.sra, SRR030223.sra, SRR030224.sra и SRR030225.sra содержащий несопоставленные короткие чтения для двух повторений из DICERex5 выборки и двух повторений из HCT116 выборки соответственно, из NCBI SRA Run Selector и преобразовал их в FASTQ-форматированные файлы с помощью NCBI SRA Toolkit.
(2) произвели файлы с форматированием SAM путем сопоставления коротких показаний с эталонным геномом человека (NCBI Build 37.5) с помощью алгоритма Bowtie [2]. Сообщаются только уникально сопоставленные чтения.
(3) сжал форматированные файлы SAM в BAM и упорядочил их сначала по ссылочному наименованию, затем по геномному положению с помощью SAMtools [3].
Этот пример также предполагает, что вы загрузили эталонный геном человека (GRCh37.p5). Можно использовать bowtie-inspect команда, чтобы восстановить человеческую ссылку непосредственно из индексов боути. Или вы можете загрузить ссылку из репозитория NCBI, исключив следующую линию:
% getgenbank('NC_000009','FileFormat','fasta','tofile','hsch9.fasta');
Чтобы исследовать охват сигнала HCT116 выборок необходимо создать BioMap. BioMap имеет интерфейс, который обеспечивает прямой доступ к отображенным коротким чтениям, хранящимся в файле в формате BAM, таким образом минимизируя объем данных, который фактически загружается в память. Используйте функцию baminfo для получения списка существующих ссылок и фактического количества коротких чтений, сопоставленных с каждым.
info = baminfo('SRR030224.bam','ScanDictionary',true); fprintf('%-35s%s\n','Reference','Number of Reads'); for i = 1:numel(info.ScannedDictionary) fprintf('%-35s%d\n',info.ScannedDictionary{i},... info.ScannedDictionaryCount(i)); end
Reference Number of Reads gi|224589800|ref|NC_000001.10| 205065 gi|224589811|ref|NC_000002.11| 187019 gi|224589815|ref|NC_000003.11| 73986 gi|224589816|ref|NC_000004.11| 84033 gi|224589817|ref|NC_000005.9| 96898 gi|224589818|ref|NC_000006.11| 87990 gi|224589819|ref|NC_000007.13| 120816 gi|224589820|ref|NC_000008.10| 111229 gi|224589821|ref|NC_000009.11| 106189 gi|224589801|ref|NC_000010.10| 112279 gi|224589802|ref|NC_000011.9| 104466 gi|224589803|ref|NC_000012.11| 87091 gi|224589804|ref|NC_000013.10| 53638 gi|224589805|ref|NC_000014.8| 64049 gi|224589806|ref|NC_000015.9| 60183 gi|224589807|ref|NC_000016.9| 146868 gi|224589808|ref|NC_000017.10| 195893 gi|224589809|ref|NC_000018.9| 60344 gi|224589810|ref|NC_000019.9| 166420 gi|224589812|ref|NC_000020.10| 148950 gi|224589813|ref|NC_000021.8| 310048 gi|224589814|ref|NC_000022.10| 76037 gi|224589822|ref|NC_000023.10| 32421 gi|224589823|ref|NC_000024.9| 18870 gi|17981852|ref|NC_001807.4| 1015 Unmapped 6805842
В этом примере вы сосредоточитесь на анализе хромосомы 9. Создайте BioMap для этих двух HCT116 выборка повторяется.
bm_hct116_1 = BioMap('SRR030224.bam','SelectRef','gi|224589821|ref|NC_000009.11|') bm_hct116_2 = BioMap('SRR030225.bam','SelectRef','gi|224589821|ref|NC_000009.11|')
bm_hct116_1 =
BioMap with properties:
SequenceDictionary: 'gi|224589821|ref|NC_000009.11|'
Reference: [106189x1 File indexed property]
Signature: [106189x1 File indexed property]
Start: [106189x1 File indexed property]
MappingQuality: [106189x1 File indexed property]
Flag: [106189x1 File indexed property]
MatePosition: [106189x1 File indexed property]
Quality: [106189x1 File indexed property]
Sequence: [106189x1 File indexed property]
Header: [106189x1 File indexed property]
NSeqs: 106189
Name: ''
bm_hct116_2 =
BioMap with properties:
SequenceDictionary: 'gi|224589821|ref|NC_000009.11|'
Reference: [107586x1 File indexed property]
Signature: [107586x1 File indexed property]
Start: [107586x1 File indexed property]
MappingQuality: [107586x1 File indexed property]
Flag: [107586x1 File indexed property]
MatePosition: [107586x1 File indexed property]
Quality: [107586x1 File indexed property]
Sequence: [107586x1 File indexed property]
Header: [107586x1 File indexed property]
NSeqs: 107586
Name: ''
Использование алгоритма раскладывания, предоставленного getBaseCoverage метод, можно построить график покрытия обоих повторений для начальной проверки. Для ссылки можно также добавить идеограмму для хромосомы 9 человека к графику с помощью chromosomeplot функция.
figure ha = gca; hold on n = 141213431; % length of chromosome 9 [cov,bin] = getBaseCoverage(bm_hct116_1,1,n,'binWidth',100); h1 = plot(bin,cov,'b'); % plots the binned coverage of bm_hct116_1 [cov,bin] = getBaseCoverage(bm_hct116_2,1,n,'binWidth',100); h2 = plot(bin,cov,'g'); % plots the binned coverage of bm_hct116_2 chromosomeplot('hs_cytoBand.txt', 9, 'AddToPlot', ha) % plots an ideogram along the x-axis axis(ha,[1 n 0 100]) % zooms-in the y-axis fixGenomicPositionLabels(ha) % formats tick labels and adds datacursors legend([h1 h2],'HCT116-1','HCT116-2','Location','NorthEast') ylabel('Coverage') title('Coverage for two replicates of the HCT116 sample') fig = gcf; fig.Position = max(fig.Position,[0 0 900 0]); % resize window
Поскольку короткие чтения представляют метилированные области ДНК, существует корреляция между выровненным покрытием и метилированием ДНК. Наблюдают повышенное метилирование ДНК, близкое к хромосомным теломерам; известно, что существует связь между метилированием ДНК и ролью теломер для поддержания целостности хромосом. На графике покрытия вы также можете увидеть длинную щель над центромером хромосомы. Это связано с повторяющимися последовательностями, присутствующими в центромере, которые мешают нам выравнивать короткие чтения в уникальное положение в этой области. В наборах данных, используемых в этом примере, только около 30% коротких показаний были однозначно сопоставлены с геномом ссылки.
Метилирование ДНК обычно происходит в динуклеотидах CpG. Изменение шаблонов метилирования ДНК может привести к транскрипционному молчанию, особенно на островах промотора генов CpG. Но также известно, что метилирование ДНК может блокировать связывание CTCF и может замолчать транскрипцию miRNA среди других соответствующих функций. В целом ожидается, что отображенные показания должны предпочтительно совпадать с богатыми CpG областями.
Загрузите хромосому 9 человека из файла ссылки hs37.fasta. В данном примере принято, что вы восстановили ссылку из индексов Боути, используя bowtie-inspect команда; поэтому hs37.fasta содержит все хромосомы человека. Чтобы загрузить только хромосому 9, можно использовать пару опция nave-value BLOCKREAD с fastaread функция.
chr9 = fastaread('hs37.fasta','blockread',9); chr9.Header
ans =
'gi|224589821|ref|NC_000009.11| Homo sapiens chromosome 9, GRCh37 primary reference assembly'
Используйте cpgisland функция для поиска кластеров CpG. Используя стандартное определение для островов CpG [4], 200 или более островов bp с 60% или более соотношением CpGobserved/CpGexpected, приводит к 1682 островам GpG, обнаруженным в хромосоме 9.
cpgi = cpgisland(chr9.Sequence)
cpgi =
struct with fields:
Starts: [1x1682 double]
Stops: [1x1682 double]
Используйте getCounts метод для вычисления отношения выровненных основ, находящихся внутри островов CpG. Для первого повторения HCT116 выборки отношение близко к 45%.
aligned_bases_in_CpG_islands = getCounts(bm_hct116_1,cpgi.Starts,cpgi.Stops,'method','sum') aligned_bases_total = getCounts(bm_hct116_1,1,n,'method','sum') ratio = aligned_bases_in_CpG_islands ./ aligned_bases_total
aligned_bases_in_CpG_islands =
1724363
aligned_bases_total =
3822804
ratio =
0.4511
Можно исследовать графики покрытия высокого разрешения двух выборок, и наблюдать, как сигнал коррелирует с островами CpG. Для примера исследуйте область между 23 820 000 и 23 830 000 bp. Это область 5 'генов человека, ELAVL2.
r1 = 23820001; % set the region limits r2 = 23830000; fhELAVL2 = figure; % keep the figure handle to use it later hold on % plot high-resolution coverage of bm_hct116_1 h1 = plot(r1:r2,getBaseCoverage(bm_hct116_1,r1,r2,'binWidth',1),'b'); % plot high-resolution coverage of bm_hct116_2 h2 = plot(r1:r2,getBaseCoverage(bm_hct116_2,r1,r2,'binWidth',1),'g'); % mark the CpG islands within the [r1 r2] region for i = 1:numel(cpgi.Starts) if cpgi.Starts(i)>r1 && cpgi.Stops(i)<r2 % is CpG island inside [r1 r2]? px = [cpgi.Starts([i i]) cpgi.Stops([i i])]; % x-coordinates for patch py = [0 max(ylim) max(ylim) 0]; % y-coordinates for patch hp = patch(px,py,'r','FaceAlpha',.1,'EdgeColor','r','Tag','cpgi'); end end axis([r1 r2 0 20]) % zooms-in the y-axis fixGenomicPositionLabels(gca) % formats tick labels and adds datacursors legend([h1 h2 hp],'HCT116-1','HCT116-2','CpG Islands') ylabel('Coverage') xlabel('Chromosome 9 position') title('Coverage for two replicates of the HCT116 sample')
Чтобы найти области, которые содержат больше отображенных показаний, чем ожидалось бы случайно, можно следовать аналогичному подходу, описанному Serre et al. [1]. Количество отсчётов для неперекрывающихся смежных окон 100 bp статистически смоделировано.
Во-первых, используйте getCounts метод для подсчета количества сопоставленных чтений, которые начинаются в каждом окне. В этом примере вы используете подход раскладывания, который рассматривает только начальное положение каждого отображенного чтения, следуя подходу Serre et al. Однако вы также можете использовать OVERLAP и METHOD Пары "имя-значение" в getCounts для вычисления более точной статистики. Например, чтобы получить максимальное покрытие для каждого окна с учетом разрешения базовой пары, установите OVERLAP по 1 и METHOD на MAX.
n = numel(chr9.Sequence); % length of chromosome w = 1:100:n; % windows of 100 bp counts_1 = getCounts(bm_hct116_1,w,w+99,'independent',true,'overlap','start'); counts_2 = getCounts(bm_hct116_2,w,w+99,'independent',true,'overlap','start');
Во-первых, попробуйте смоделировать счетчики, принимая, что все окна с счетчиками являются биологически значимыми и, следовательно, из одного и того же распределения. Используйте отрицательное биономиальное распределение, чтобы соответствовать модели, данным отсчета.
nbp = nbinfit(counts_1);
Постройте график подобранной модели по гистограмме эмпирических данных.
figure hold on emphist = histc(counts_1,0:100); % calculate the empirical distribution bar(0:100,emphist./sum(emphist),'c','grouped') % plot histogram plot(0:100,nbinpdf(0:100,nbp(1),nbp(2)),'b','linewidth',2); % plot fitted model axis([0 50 0 .001]) legend('Empirical Distribution','Negative Binomial Fit') ylabel('Frequency') xlabel('Counts') title('Frequency of counts, 100bp windows (HCT116-1)')
Плохой подбор кривой указывает, что наблюдаемое распределение может быть связано со смесью двух моделей, одной, которая представляет фон, и одной, которая представляет данные подсчета в метилированных окнах ДНК.
Более реалистичным сценарием было бы предположение, что окна с небольшим количеством сопоставленных чтений являются главным образом фоном (или нулевой моделью). Serre et al. предположил, что 100-bp окна, содержащие четыре или более чтений, вряд ли будут сгенерированы случайно. Чтобы оценить хорошее приближение к нулевой модели, можно подогнать левое тело эмпирического распределения к усеченному отрицательному биномиальному распределению. Для соответствия усеченному распределению используйте mle функция. Сначала вам нужно задать анонимную функцию, которая задает правую усеченную версию nbinpdf.
rtnbinpdf = @(x,p1,p2,t) nbinpdf(x,p1,p2) ./ nbincdf(t-1,p1,p2);
Определите функцию аппроксимации с помощью другой анонимной функции.
rtnbinfit = @(x2,t) mle(x2,'pdf',@(x3,p1,p2) rtnbinpdf(x3,p1,p2,t),'start',nbinfit(x2),'lower',[0 0]);
Прежде чем подгонять реальные данные, давайте оценим процедуру подгонки с некоторыми выборочными данными из известного распределения.
nbp = [0.5 0.2]; % Known coefficients x = nbinrnd(nbp(1),nbp(2),10000,1); % Random sample trun = 6; % Set a truncation threshold nbphat1 = nbinfit(x); % Fit non-truncated model to all data nbphat2 = nbinfit(x(x<trun)); % Fit non-truncated model to truncated data (wrong) nbphat3 = rtnbinfit(x(x<trun),trun); % Fit truncated model to truncated data figure hold on emphist = histc(x,0:100); % Calculate the empirical distribution bar(0:100,emphist./sum(emphist),'c','grouped') % plot histogram h1 = plot(0:100,nbinpdf(0:100,nbphat1(1),nbphat1(2)),'b-o','linewidth',2); h2 = plot(0:100,nbinpdf(0:100,nbphat2(1),nbphat2(2)),'r','linewidth',2); h3 = plot(0:100,nbinpdf(0:100,nbphat3(1),nbphat3(2)),'g','linewidth',2); axis([0 25 0 .2]) legend([h1 h2 h3],'Neg-binomial fitted to all data',... 'Neg-binomial fitted to truncated data',... 'Truncated neg-binomial fitted to truncated data') ylabel('Frequency') xlabel('Counts')
Для двух повторений выборки HCT116 подбирайте правое усеченное отрицательное биномиальное распределение к наблюдаемой нулевой модели с помощью rtnbinfit анонимная функция, ранее заданная.
trun = 4; % Set a truncation threshold (as in [1]) pn1 = rtnbinfit(counts_1(counts_1<trun),trun); % Fit to HCT116-1 counts pn2 = rtnbinfit(counts_2(counts_2<trun),trun); % Fit to HCT116-2 counts
Вычислите значение p для каждого окна в нулевом распределении.
pval1 = 1 - nbincdf(counts_1,pn1(1),pn1(2)); pval2 = 1 - nbincdf(counts_2,pn2(1),pn2(2));
Вычислите скорость ложного обнаружения с помощью mafdr функция. Используйте пару "имя-значение" BHFDR использование процедуры линейного шага вверх (LSU) ([6]) для вычисления скорректированных p-значений FDR. Установка FDR < 0,01 позволяет вам идентифицировать окна с 100 б.п., которые значительно метилированы.
fdr1 = mafdr(pval1,'bhfdr',true); fdr2 = mafdr(pval2,'bhfdr',true); w1 = fdr1<.01; % logical vector indicating significant windows in HCT116-1 w2 = fdr2<.01; % logical vector indicating significant windows in HCT116-2 w12 = w1 & w2; % logical vector indicating significant windows in both replicates Number_of_sig_windows_HCT116_1 = sum(w1) Number_of_sig_windows_HCT116_2 = sum(w2) Number_of_sig_windows_HCT116 = sum(w12)
Number_of_sig_windows_HCT116_1 =
1662
Number_of_sig_windows_HCT116_2 =
1674
Number_of_sig_windows_HCT116 =
1346
В целом, вы идентифицировали 1662 и 1674 неперекрывающихся 100-bp окон в двух репликатах выборок HCT116, что указывает на наличие значительных доказательств метилирования ДНК. Существует 1346 окон, которые значительны в обоих репликатах.
Для примера, снова глядя в промоторную область ELAVL2 гена человека, можно заметить, что в обеих выборки репликатах было отмечено значительное количество окон с 100 п.о.
figure(fhELAVL2) % bring back to focus the previously plotted figure plot(w(w1)+50,counts_1(w1),'bs', 'HandleVisibility','off') % plot significant windows in HCT116-1 plot(w(w2)+50,counts_2(w2),'gs', 'HandleVisibility','off') % plot significant windows in HCT116-2 axis([r1 r2 0 100]) title('Significant 100-bp windows in both replicates of the HCT116 sample')
Метилирование ДНК участвует в модуляции экспрессии генов. Например, хорошо известно, что гиперметилирование связано с инактивацией нескольких генов-супрессоров опухолей. Вы можете изучить в этом наборе данных метилирование областей промотора генов.
Во-первых, загрузите из Ensembl вкладки таблицу с разделенными значениями (TSV) со всеми генами, кодирующими белок, в текстовый файл, ensemblmart_genes_hum37.txt. В данном примере мы используем Ensembl release 64. Используя сервис Ensembl BioMart, можно выбрать таблицу со следующими атрибутами: имя хромосомы, биотип гена, имя гена, начало/конец гена и направление цепи.
Используйте предоставленную функцию helper ensemblmart2gff для преобразования загруженного TSV-файла в GFF-форматированный файл. Затем используйте GFFAnnotation загрузить файл в MATLAB и создать подмножество с генами, присутствующими только в хромосоме 9. Это приводит к 800 аннотированным генам, кодирующим белок, в базе данных Ensembl.
GFFfilename = ensemblmart2gff('ensemblmart_genes_hum37.txt'); a = GFFAnnotation(GFFfilename) a9 = getSubset(a,'reference','9') numGenes = a9.NumEntries
a =
GFFAnnotation with properties:
FieldNames: {1x9 cell}
NumEntries: 21184
a9 =
GFFAnnotation with properties:
FieldNames: {1x9 cell}
NumEntries: 800
numGenes =
800
Найдите области промотора для каждого гена. В этом примере мы рассматриваем проксимальный промотор как область -500/100 в восходящем направлении.
downstream = 500; upstream = 100; geneDir = strcmp(a9.Strand,'+'); % logical vector indicating strands in the forward direction % calculate promoter's start position for genes in the forward direction promoterStart(geneDir) = a9.Start(geneDir) - downstream; % calculate promoter's end position for genes in the forward direction promoterStop(geneDir) = a9.Start(geneDir) + upstream; % calculate promoter's start position for genes in the reverse direction promoterStart(~geneDir) = a9.Stop(~geneDir) - upstream; % calculate promoter's end position for genes in the reverse direction promoterStop(~geneDir) = a9.Stop(~geneDir) + downstream;
Использование dataset как контейнер для информации о промоторе, так как мы можем позже добавить новые столбцы для хранения количества генов и p-значений.
promoters = dataset({a9.Feature,'Gene'});
promoters.Strand = char(a9.Strand);
promoters.Start = promoterStart';
promoters.Stop = promoterStop';
Найдите гены со значительным метилированием ДНК в промоторной области, посмотрев на количество отображенных коротких показаний, которые перекрывают, по крайней мере, одну пару оснований в определенной промоторной области.
promoters.Counts_1 = getCounts(bm_hct116_1,promoters.Start,promoters.Stop,... 'overlap',1,'independent',true); promoters.Counts_2 = getCounts(bm_hct116_2,promoters.Start,promoters.Stop,... 'overlap',1,'independent',true);
Подгонка нулевого распределения для каждой выборки репликации и вычисление p-значений:
trun = 5; % Set a truncation threshold pn1 = rtnbinfit(promoters.Counts_1(promoters.Counts_1<trun),trun); % Fit to HCT116-1 promoter counts pn2 = rtnbinfit(promoters.Counts_2(promoters.Counts_2<trun),trun); % Fit to HCT116-2 promoter counts promoters.pval_1 = 1 - nbincdf(promoters.Counts_1,pn1(1),pn1(2)); % p-value for every promoter in HCT116-1 promoters.pval_2 = 1 - nbincdf(promoters.Counts_2,pn2(1),pn2(2)); % p-value for every promoter in HCT116-2 Number_of_sig_promoters = sum(promoters.pval_1<.01 & promoters.pval_2<.01) Ratio_of_sig_methylated_promoters = Number_of_sig_promoters./numGenes
Number_of_sig_promoters =
74
Ratio_of_sig_methylated_promoters =
0.0925
Заметим, что только 74 (из 800) гена в хромосоме 9 имеют значительно метилированные области ДНК (pval < 0,01 в обоих репликатах). Отобразите отчет о 30 генах с наиболее значимыми метилированными промоторными областями.
[~,order] = sort(promoters.pval_1.*promoters.pval_2); promoters(order(1:30),[1 2 3 4 5 7 6 8])
ans =
Gene Strand Start Stop Counts_1
{'DMRT3' } + 976464 977064 223
{'CNTFR' } - 34590021 34590621 219
{'GABBR2' } - 101471379 101471979 404
{'CACNA1B' } + 140771741 140772341 454
{'BARX1' } - 96717554 96718154 264
{'FAM78A' } - 134151834 134152434 497
{'FOXB2' } + 79634071 79634671 163
{'TLE4' } + 82186188 82186788 157
{'ASTN2' } - 120177248 120177848 141
{'FOXE1' } + 100615036 100615636 149
{'MPDZ' } - 13279489 13280089 129
{'PTPRD' } - 10612623 10613223 145
{'PALM2-AKAP2'} + 112542089 112542689 134
{'FAM69B' } + 139606522 139607122 112
{'WNK2' } + 95946698 95947298 108
{'IGFBPL1' } - 38424344 38424944 110
{'AKAP2' } + 112542269 112542869 107
{'C9orf4' } - 111929471 111930071 102
{'COL5A1' } + 137533120 137533720 84
{'LHX3' } - 139096855 139097455 74
{'OLFM1' } + 137966768 137967368 75
{'NPR2' } + 35791651 35792251 68
{'DBC1' } - 122131645 122132245 61
{'SOHLH1' } - 138591274 138591874 56
{'PIP5K1B' } + 71320075 71320675 59
{'PRDM12' } + 133539481 133540081 53
{'ELAVL2' } - 23826235 23826835 50
{'ZFP37' } - 115818939 115819539 59
{'RP11-35N6.1'} + 103790491 103791091 60
{'DMRT2' } + 1049854 1050454 54
pval_1 Counts_2 pval_2
6.6613e-16 253 5.5511e-16
6.6613e-16 226 5.5511e-16
6.6613e-16 400 5.5511e-16
6.6613e-16 408 5.5511e-16
6.6613e-16 286 5.5511e-16
6.6613e-16 499 5.5511e-16
1.4e-13 165 6.0352e-13
3.5649e-13 151 4.7347e-12
4.3566e-12 163 8.0969e-13
1.2447e-12 133 6.7598e-11
2.8679e-11 148 7.3682e-12
2.3279e-12 127 1.6448e-10
1.3068e-11 135 5.0276e-11
4.1911e-10 144 1.3295e-11
7.897e-10 125 2.2131e-10
5.7523e-10 114 1.1364e-09
9.2538e-10 106 3.7513e-09
2.0467e-09 96 1.6795e-08
3.6266e-08 97 1.4452e-08
1.8171e-07 91 3.5644e-08
1.5457e-07 69 1.0074e-06
4.8093e-07 73 5.4629e-07
1.5082e-06 62 2.9575e-06
3.4322e-06 67 1.3692e-06
2.0943e-06 63 2.5345e-06
5.6364e-06 61 3.4518e-06
9.2778e-06 62 2.9575e-06
2.0943e-06 47 3.0746e-05
1.7771e-06 42 6.8037e-05
4.7762e-06 46 3.6016e-05
Serre et al. [1] сообщили, что в этих наборах данных приблизительно 90% уникально отображенных показаний выпадают за пределы 5 'областей промотора генов. Используя подобный подход, как и прежде, можно найти гены, которые имеют межгенные метилированные области. Чтобы компенсировать различные длины генов, можно использовать максимальное покрытие, вычисленное по основаниям, вместо необработанного количества отображенных коротких показаний. Другой альтернативный подход к нормализации счетчиков по длине генов - установить METHOD Пара "имя-значение" с rpkm в getCounts функция.
intergenic = dataset({a9.Feature,'Gene'});
intergenic.Strand = char(a9.Strand);
intergenic.Start = a9.Start;
intergenic.Stop = a9.Stop;
intergenic.Counts_1 = getCounts(bm_hct116_1,intergenic.Start,intergenic.Stop,...
'overlap','full','method','max','independent',true);
intergenic.Counts_2 = getCounts(bm_hct116_2,intergenic.Start,intergenic.Stop,...
'overlap','full','method','max','independent',true);
trun = 10; % Set a truncation threshold
pn1 = rtnbinfit(intergenic.Counts_1(intergenic.Counts_1<trun),trun); % Fit to HCT116-1 intergenic counts
pn2 = rtnbinfit(intergenic.Counts_2(intergenic.Counts_2<trun),trun); % Fit to HCT116-2 intergenic counts
intergenic.pval_1 = 1 - nbincdf(intergenic.Counts_1,pn1(1),pn1(2)); % p-value for every intergenic region in HCT116-1
intergenic.pval_2 = 1 - nbincdf(intergenic.Counts_2,pn2(1),pn2(2)); % p-value for every intergenic region in HCT116-2
Number_of_sig_genes = sum(intergenic.pval_1<.01 & intergenic.pval_2<.01)
Ratio_of_sig_methylated_genes = Number_of_sig_genes./numGenes
[~,order] = sort(intergenic.pval_1.*intergenic.pval_2);
intergenic(order(1:30),[1 2 3 4 5 7 6 8])
Number_of_sig_genes =
62
Ratio_of_sig_methylated_genes =
0.0775
ans =
Gene Strand Start Stop Counts_1
{'AL772363.1'} - 140762377 140787022 106
{'CACNA1B' } + 140772241 141019076 106
{'SUSD1' } - 114803065 114937688 88
{'C9orf172' } + 139738867 139741797 99
{'NR5A1' } - 127243516 127269709 86
{'BARX1' } - 96713628 96717654 77
{'KCNT1' } + 138594031 138684992 58
{'GABBR2' } - 101050391 101471479 65
{'FOXB2' } + 79634571 79635869 51
{'NDOR1' } + 140100119 140113813 54
{'KIAA1045' } + 34957484 34984679 50
{'ADAMTSL2' } + 136397286 136440641 55
{'PAX5' } - 36833272 37034476 48
{'OLFM1' } + 137967268 138013025 55
{'PBX3' } + 128508551 128729656 45
{'FOXE1' } + 100615536 100618986 49
{'MPDZ' } - 13105703 13279589 51
{'ASTN2' } - 119187504 120177348 43
{'ARRDC1' } + 140500106 140509812 49
{'IGFBPL1' } - 38408991 38424444 45
{'LHX3' } - 139088096 139096955 44
{'PAPPA' } + 118916083 119164601 44
{'CNTFR' } - 34551430 34590121 41
{'DMRT3' } + 976964 991731 40
{'TUSC1' } - 25676396 25678856 46
{'ELAVL2' } - 23690102 23826335 35
{'SMARCA2' } + 2015342 2193624 36
{'GAS1' } - 89559279 89562104 34
{'GRIN1' } + 140032842 140063207 36
{'TLE4' } + 82186688 82341658 36
pval_1 Counts_2 pval_2
8.6597e-15 98 1.8763e-14
8.6597e-15 98 1.8763e-14
2.2904e-12 112 7.7716e-16
7.4718e-14 96 3.5749e-14
4.268e-12 90 2.5457e-13
7.0112e-11 62 2.569e-09
2.5424e-08 73 6.9019e-11
2.9078e-09 58 9.5469e-09
2.2131e-07 58 9.5469e-09
8.7601e-08 55 2.5525e-08
3.0134e-07 55 2.5525e-08
6.4307e-08 45 6.7163e-07
5.585e-07 49 1.8188e-07
6.4307e-08 42 1.7861e-06
1.4079e-06 51 9.4566e-08
4.1027e-07 46 4.8461e-07
2.2131e-07 42 1.7861e-06
2.6058e-06 43 1.2894e-06
4.1027e-07 36 1.2564e-05
1.4079e-06 39 4.7417e-06
1.9155e-06 36 1.2564e-05
1.9155e-06 35 1.7377e-05
4.8199e-06 37 9.0815e-06
6.5537e-06 37 9.0815e-06
1.0346e-06 31 6.3417e-05
3.0371e-05 41 2.4736e-06
2.2358e-05 40 3.4251e-06
4.1245e-05 41 2.4736e-06
2.2358e-05 38 6.5629e-06
2.2358e-05 37 9.0815e-06
Например, исследуйте профиль метилирования гена BARX1, шестого значимого гена с межгенным метилированием в предыдущем списке. Форматированный файл GTF ensemblmart_barx1.gtf содержит структурную информацию для этого гена, полученную от Ensembl с помощью сервиса BioMart.
Использование GTFAnnotation для загрузки структурной информации в MATLAB. Существует два аннотированных транскрипта для этого гена.
barx1 = GTFAnnotation('ensemblmart_barx1.gtf')
transcripts = getTranscriptNames(barx1)
barx1 =
GTFAnnotation with properties:
FieldNames: {1x11 cell}
NumEntries: 18
transcripts =
2x1 cell array
{'ENST00000253968'}
{'ENST00000401724'}
Постройте график метилирования ДНК для обеих выборок HCT116 с разрешением пары оснований. Наложите острова CpG и постройте графики экзонов для каждого из двух транскриптов вдоль дна графика.
range = barx1.getRange; r1 = range(1)-1000; % set the region limits r2 = range(2)+1000; figure hold on % plot high-resolution coverage of bm_hct116_1 h1 = plot(r1:r2,getBaseCoverage(bm_hct116_1,r1,r2,'binWidth',1),'b'); % plot high-resolution coverage of bm_hct116_2 h2 = plot(r1:r2,getBaseCoverage(bm_hct116_2,r1,r2,'binWidth',1),'g'); % mark the CpG islands within the [r1 r2] region for i = 1:numel(cpgi.Starts) if cpgi.Starts(i)>r1 && cpgi.Stops(i)<r2 % is CpG island inside [r1 r2]? px = [cpgi.Starts([i i]) cpgi.Stops([i i])]; % x-coordinates for patch py = [0 max(ylim) max(ylim) 0]; % y-coordinates for patch hp = patch(px,py,'r','FaceAlpha',.1,'EdgeColor','r','Tag','cpgi'); end end % mark the exons at the bottom of the axes for i = 1:numel(transcripts) exons = getSubset(barx1,'Transcript',transcripts{i},'Feature','exon'); for j = 1:exons.NumEntries px = [exons.Start([j j]);exons.Stop([j j])]'; % x-coordinates for patch py = [0 1 1 0]-i*2-1; % y-coordinates for patch hq = patch(px,py,'b','FaceAlpha',.1,'EdgeColor','b','Tag','exon'); end end axis([r1 r2 -numel(transcripts)*2-2 80]) % zooms-in the y-axis fixGenomicPositionLabels(gca) % formats tick labels and adds datacursors ylabel('Coverage') xlabel('Chromosome 9 position') title('High resolution coverage in the BARX1 gene') legend([h1 h2 hp hq],'HCT116-1','HCT116-2','CpG Islands','Exons','Location','NorthWest')
Наблюдайте высоко метилированную область в области промотора 5 '(самый правый остров CpG). Напомним, что для этого гена транскрипция происходит в обратной цепи. Более интересно наблюдать высоко метилированные области, которые перекрывают инициацию каждого из двух аннотированных транскриптов (два средних острова CpG).
В исследовании Serre et al. анализируют также другую клеточную линию. Новые камеры (DICERex5) получают из тех же HCT116 камер рака толстой кишки после усечения DICER1 аллелей. В литературе [5] сообщалось, что происходит локализованное изменение метилирования ДНК у небольшого числа генных промоторов. В этом примере вы найдете значительные 100-bp окна в двух выборочных репликатах DICERex5 клеток, следуя тому же подходу, что и родительские HCT116 клетки, и затем вы будете искать статистически значимые различия между двумя клеточными линиями.
Функция помощника getWindowCounts захватывает аналогичные шаги для поиска окон со значительным покрытием, как и прежде. getWindowCounts возвращает векторы с счетчиками, p-значениями и частотой ложных открытий для каждой новой репликации.
bm_dicer_1 = BioMap('SRR030222.bam','SelectRef','gi|224589821|ref|NC_000009.11|'); bm_dicer_2 = BioMap('SRR030223.bam','SelectRef','gi|224589821|ref|NC_000009.11|'); [counts_3,pval3,fdr3] = getWindowCounts(bm_dicer_1,4,w,100); [counts_4,pval4,fdr4] = getWindowCounts(bm_dicer_2,4,w,100); w3 = fdr3<.01; % logical vector indicating significant windows in DICERex5_1 w4 = fdr4<.01; % logical vector indicating significant windows in DICERex5-2 w34 = w3 & w4; % logical vector indicating significant windows in both replicates Number_of_sig_windows_DICERex5_1 = sum(w3) Number_of_sig_windows_DICERex5_2 = sum(w4) Number_of_sig_windows_DICERex5 = sum(w34)
Number_of_sig_windows_DICERex5_1 =
908
Number_of_sig_windows_DICERex5_2 =
1041
Number_of_sig_windows_DICERex5 =
759
Для выполнения дифференциального анализа используются окна 100-bp, которые значительны по крайней мере в одной из выборок (HCT116 или DICERex5).
wd = w34 | w12; % logical vector indicating windows included in the diff. analysis counts = [counts_1(wd) counts_2(wd) counts_3(wd) counts_4(wd)]; ws = w(wd); % window start for each row in counts
Используйте функцию manorm для нормализации данных. The PERCENTILE Пара "имя-значение" позволяет вам отфильтровать окна с очень большим количеством отсчётов при нормализации, поскольку эти окна в основном вызваны программными продуктами, такими как повторяющиеся области в ссылку хромосоме.
counts_norm = round(manorm(counts,'percentile',90).*100);
Используйте функцию mattest для выполнения t-теста с двумя образцами для идентификации дифференциально покрытых окон из двух различных клеточных линий.
pval = mattest(counts_norm(:,[1 2]),counts_norm(:,[3 4]),'bootstrap',true,... 'showhist',true,'showplot',true);
Составьте отчет с 25 наиболее значимыми дифференциально покрытыми окнами. При создании отчета используйте функцию helper findClosestGene определить, является ли окно межгенным, внутригенным или находится ли оно в проксимальной промоторной области.
[~,ord] = sort(pval); fprintf('Window Pos Type p-value HCT116 DICERex5\n\n'); for i = 1:25 j = ord(i); [~,msg] = findClosestGene(a9,[ws(j) ws(j)+99]); fprintf('%10d %-25s %7.6f%5d%5d %5d%5d\n', ... ws(j),msg,pval(j),counts_norm(j,:)); end
Window Pos Type p-value HCT116 DICERex5
140311701 Intergenic (EXD3) 0.000026 13 13 104 105
139546501 Intragenic 0.001827 21 21 91 93
10901 Intragenic 0.002671 258 257 434 428
120176801 Intergenic (ASTN2) 0.002733 266 270 155 155
139914801 Intergenic (ABCA2) 0.002980 64 63 26 25
126128501 Intergenic (CRB2) 0.003193 94 93 129 130
71939501 Prox. Promoter (FAM189A2) 0.005549 107 101 0 0
124461001 Intergenic (DAB2IP) 0.005618 77 76 39 37
140086501 Intergenic (TPRN) 0.006520 47 42 123 124
79637201 Intragenic 0.007512 52 51 32 31
136470801 Intragenic 0.007512 52 51 32 31
140918001 Intergenic (CACNA1B) 0.008115 176 169 71 68
100615901 Intergenic (FOXE1) 0.008346 262 253 123 118
98221901 Intergenic (PTCH1) 0.009934 26 30 104 99
138730601 Intergenic (CAMSAP1) 0.010273 26 21 97 93
89561701 Intergenic (GAS1) 0.010336 77 76 6 12
977401 Intergenic (DMRT3) 0.010369 236 245 129 124
37002601 Intergenic (PAX5) 0.010559 133 127 207 211
139744401 Intergenic (PHPT1) 0.010869 47 46 32 31
126771301 Intragenic 0.011458 43 46 97 93
93922501 Intragenic 0.011486 34 34 149 161
94187101 Intragenic 0.011507 73 80 6 6
136044401 Intragenic 0.011567 39 34 110 105
139611201 Intergenic (FAM69B) 0.011567 39 34 110 105
139716201 Intergenic (C9orf86) 0.011832 73 72 136 130
Постройте график профиля метилирования ДНК для промоторной области генной FAM189A2, наиболее значительной дифференциально покрытой области промотора по сравнению с предыдущим списком. Наложите острова CpG и ген FAM189A2.
range = getRange(getSubset(a9,'Feature','FAM189A2')); r1 = range(1)-1000; r2 = range(2)+1000; figure hold on % plot high-resolution coverage of all replicates h1 = plot(r1:r2,getBaseCoverage(bm_hct116_1,r1,r2,'binWidth',1),'b'); h2 = plot(r1:r2,getBaseCoverage(bm_hct116_2,r1,r2,'binWidth',1),'g'); h3 = plot(r1:r2,getBaseCoverage(bm_dicer_1,r1,r2,'binWidth',1),'r'); h4 = plot(r1:r2,getBaseCoverage(bm_dicer_2,r1,r2,'binWidth',1),'m'); % mark the CpG islands within the [r1 r2] region for i = 1:numel(cpgi.Starts) if cpgi.Starts(i)>r1 && cpgi.Stops(i)<r2 % is CpG island inside [r1 r2]? px = [cpgi.Starts([i i]) cpgi.Stops([i i])]; % x-coordinates for patch py = [0 max(ylim) max(ylim) 0]; % y-coordinates for patch hp = patch(px,py,'r','FaceAlpha',.1,'EdgeColor','r','Tag','cpgi'); end end % mark the gene at the bottom of the axes px = range([1 1 2 2]); py = [0 1 1 0]-2; hq = patch(px,py,'b','FaceAlpha',.1,'EdgeColor','b','Tag','gene'); axis([r1 r1+4000 -4 30]) % zooms-in fixGenomicPositionLabels(gca) % formats tick labels and adds datacursors ylabel('Coverage') xlabel('Chromosome 9 position') title('DNA Methylation profiles along the promoter region of the FAM189A2 gene') legend([h1 h2 h3 h4 hp hq],... 'HCT116-1','HCT116-2','DICERex5-1','DICERex5-2','CpG Islands','FAM189A2 Gene',... 'Location','NorthEast')
Заметьте, что острова CpG явно неметилированы для обоих DICERex5 повторяются.
[1] Serre, D., Lee, B.H., and Ting A.H. «, MBD-выделенное секвенирование генома обеспечивает высокопроизводительное и комплексное исследование метилирования ДНК в геноме человека» Nucleic Acids Research 38 (2): 391-9, 2010.
[2] Langmead, B., Trapnell, C., Pop, M. and Salzberg, S.L., «Ultrafast and Memory-effective Выравнивания of Short ДНК Sequences to the Human Genome», Genome Biology, 10 (3): R25, 2009.
[3] Li, H., et al., «The Sequence Alignment/map (SAM) Формата и SAMtools», Bioinformatics, 25 (16): 2078-9, 2009.
[4] Gardiner-Garden, M. and Frommer, M «., острова CpG в геномах позвоночных», Journal of Molecular Biology, 196 (2): 261-82, 1987.
[5] Ting, A.H., et al., «A Requirement for DICER to Maintain Full Promoter CpG Island Hypermethilation in Human Cancer Камер», Cancer Research, 68 (8): 2570-5, 2008.
[6] Benjamini, Y. and Hochberg, Y., «Controlling the false discovery rate: a practical and english approach to multiple проверки», Journal of the Royal Statistical Society, 57 (1): 289-300, 1995.