Исследование различий всего генома в профилях метилирования ДНК
В этом примере показано, как выполнить анализ всего генома метилирования ДНК в человеке при помощи секвенирования генома.
Примечание: Для расширенной производительности, MathWorks рекомендует, чтобы вы запустили этот пример на 64-битной платформе, потому что объем потребляемой памяти близко к 2 Гбайт. На 32-битной платформе, если вы получаете "Out of memory" ошибки при выполнении этого примера, попытайтесь увеличить виртуальную память (или область подкачки) операционной системы или попытайтесь установить 3GB переключатель (только 32-битный Windows® XP). Эти методы описаны в этом документе.
Введение
Метилирование ДНК является эпигенетической модификацией, которая модулирует экспрессию гена и обслуживание геномной организации в процессах болезни и нормальном. Метилирование ДНК может задать различные состояния ячейки, и это является наследуемым во время репликации ячейки. Отклоняющиеся шаблоны метилирования ДНК были сопоставлены с генами-супрессорами рака и опухоли.
В этом примере вы исследуете профили метилирования ДНК двух человеческих раковых клеток: родительские ячейки рака толстой кишки HCT116 и ячейки DICERex5. Ячейки DICERex5 выведены из ячеек HCT116 после усечения аллелей DICER1. Серр и др. в [1] предложил изучить профили метилирования ДНК при помощи белка MBD2 как метил CpG обязательная область и впоследствии использовал высокопроизводительное секвенирование (HTseq). Этот метод, обычно знают как MBD-Seq. Короткие чтения для два реплицируют этих двух выборок, были представлены архиву SRA NCBI авторами [1] года. Существуют другие технологии, доступные, чтобы опросить состояние метилирования ДНК сайтов CpG в сочетании с HTseq, например, MeDIP-seq или использование ферментов ограничения. Можно также анализировать этот тип наборов данных после подхода, представленного в этом примере.
Наборы данных
Можно получить несопоставленные чтения одно конца для четырех экспериментов секвенирования от NCBI. Короткие чтения были произведены с помощью Генома Illumina® Анализатор II. Средним размером вставки является 120 BP, и продолжительностью коротких чтений является 36 BP.
Этот пример принимает что вы:
(1) загруженный файлы SRR030222.sra, SRR030223.sra, SRR030224.sra и SRR030225.sra содержание несопоставленных коротких чтений для два реплицирует от выборки DICERex5, и два реплицирует от выборки HCT116 соответственно, от Селектора Запуска SRA NCBI и преобразовал их в FASTQ-отформатированные файлы с помощью Инструментария SRA NCBI.
(2) произведенные SAM-отформатированные файлы путем отображения коротких чтений со ссылочным геномом человека (Сборка NCBI 37.5) использование Галстука-бабочки [2] алгоритм. О только исключительно сопоставленных чтениях сообщают.
(3) сжатый SAM отформатировал файлы к BAM и заказал им ссылочным именем сначала, затем геномным положением при помощи SAMtools [3].
Этот пример также принимает, что вы загрузили ссылочный геном человека (GRCh37.p5). Можно использовать bowtie-inspect команда, чтобы восстановить человеческую ссылку непосредственно от индексов галстука-бабочки. Или можно загрузить ссылку с репозитория NCBI путем некомментария следующей линии:
% getgenbank('NC_000009','FileFormat','fasta','tofile','hsch9.fasta');
Создание интерфейса MATLAB® к отформатированным BAM файлам
Чтобы исследовать покрытие сигнала выборок 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
Поскольку короткие чтения представляют methylated области ДНК, существует корреляция между выровненным покрытием и метилированием ДНК. Наблюдайте увеличенное метилирование ДНК близко к теломерам хромосомы; известно, что существует ассоциация между метилированием ДНК и ролью теломер для поддержания целостности хромосом. В графике покрытия можно также видеть длинный разрыв по центромере хромосомы. Это происходит из-за повторяющихся последовательностей, существующих в центромере, которые препятствуют тому, чтобы мы выровняли короткие чтения к уникальной позиции в этой области. В наборах данных, используемых в этом примере, только приблизительно 30% коротких чтений были исключительно сопоставлены со ссылочным геномом.
Корреляция островов CpG и метилирования ДНК
Метилирование ДНК обычно происходит в динуклеотидах CpG. Изменение шаблонов метилирования ДНК может привести к транскрипционному глушению, особенно в генном промоутере острова CpG. Но, также известно, что метилирование ДНК может блокировать привязку CTCF и может заставить miRNA запись замолчать среди других соответствующих функций. В общем случае ожидается, что сопоставленные чтения должны предпочтительно выровнять к CpG богатые области.
Загрузите человеческую хромосому 9 от ссылочного файла hs37.fasta. В данном примере это принято, что вы восстановили ссылку с индексов Галстука-бабочки с помощью bowtie-inspect команда; поэтому hs37.fasta содержит все человеческие хромосомы. Чтобы загрузить только хромосому 9, можно использовать пару значения нефа опции 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, приводит к островам на 1 682 грана на галлон, найденным в хромосоме 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')
Статистическое моделирование данных о количестве
Найти области, которые содержат более сопоставленные чтения, чем, ожидалось бы случайно, можно следовать за аналогичным подходом к тому, описанному Серром и др. [1]. Количество счетов для неналожения непрерывных 100 окон BP статистически моделируется.
Во-первых, используйте getCounts метод, чтобы считать количество сопоставленных чтений, которые запускаются в каждом окне. В этом примере вы используете подход раскладывания, который рассматривает только положение запуска каждого сопоставленного чтения, после подхода Серра и др. Однако можно также использовать 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');
Во-первых, попытайтесь смоделировать количества, принимающие, что все окна с количествами являются биологически значительными и поэтому от того же распределения. Используйте отрицательное bionomial распределение, чтобы подобрать модель данные о количестве.
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)')
Плохой подбор кривой указывает, что наблюдаемое распределение может произойти из-за смеси двух моделей, та, которая представляет фон и тот, который представляет данные о количестве в methylated окнах DNA.
Более реалистический сценарий должен был бы принять, что окна с небольшим количеством сопоставленных чтений являются в основном фоном (или пустая модель). Серр и др. принял, что окна с 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')
Идентификация значительных областей Methylated
Поскольку эти два реплицируют выборки 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 использовать процедуру ([6]) линейного шага (LSU), чтобы вычислить ФРГ настроило p-значения. Установка ФРГ <0,01 разрешения вы, чтобы идентифицировать окна с 100 BP, которые значительно являются methylated.
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
В целом, вы идентифицировали, что 1 662 и 1 674 неперекрывающихся окна с 100 BP в этих двух реплицируют выборок HCT116, который указывает, что существует значительное доказательство метилирования ДНК. Существует 1 346 окон, которые являются значительными в обоих, реплицирует.
Например, смотря снова в области промоутера человеческого гена ELAVL2 можно заметить, что в обеих выборках реплицирует, несколько окон с 100 BP были отмечены значительные.
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')
Нахождение генов со значительными областями промоутера Methylated
Метилирование ДНК вовлечено в модуляцию экспрессии гена. Например, известно, что гиперметилирование сопоставлено с деактивацией нескольких генов-супрессоров опухоли. Можно изучить в этом наборе данных метилирование генных областей промоутера.
Во-первых, загрузите с Ensembl таблицу "перейдите разделенного значения" (TSV) со всеми генами кодирования белка к текстовому файлу, ensemblmart_genes_hum37.txt. В данном примере мы используем релиз 64 Ensembl. Используя сервис BioMart Энсембла, можно выбрать таблицу со следующими атрибутами: имя хромосомы, генный биотип, название гена, ген запускает/заканчивает и скручивает направление.
Используйте обеспеченную функцию помощника 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 значительно имеют ДНК methylated, области (pval <0.01 в обоих реплицирует). Отобразите отчет этих 30 генов со старшими значащими methylated областями промоутера.
[~,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 6.6613e-16
6.6613e-16 226 6.6613e-16
6.6613e-16 400 6.6613e-16
6.6613e-16 408 6.6613e-16
6.6613e-16 286 6.6613e-16
6.6613e-16 499 6.6613e-16
1.4e-13 165 6.0363e-13
3.5649e-13 151 4.7348e-12
4.3566e-12 163 8.098e-13
1.2447e-12 133 6.7598e-11
2.8679e-11 148 7.3683e-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
Нахождение Межгенных областей, которые значительно являются Methylated
Серр и др. [1] сообщил, что в этих наборах данных приблизительно 90% исключительно сопоставленных чтений выходят за пределы 5' генных областей промоутера. Используя аналогичный подход как прежде, можно найти гены, которые имеют межгенные methylated области. Чтобы компенсировать различные длины генов, можно использовать максимальное покрытие, вычисленную основу основой, вместо необработанного количества сопоставленных коротких чтений. Другой альтернативный подход, чтобы нормировать количества генной длиной должен установить 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.3267e-15 98 1.8097e-14
8.3267e-15 98 1.8097e-14
2.2901e-12 112 1.1102e-16
7.4385e-14 96 3.5083e-14
4.2677e-12 90 2.5391e-13
7.0112e-11 62 2.5691e-09
2.5424e-08 73 6.9018e-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.0816e-06
6.5537e-06 37 9.0816e-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.0816e-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')
Наблюдайте высоко methylated область в 5' областях промоутера (самый правый остров CpG). Вспомните, что для этой транскрипции генов происходит в противоположной скрутке. Более интересный, наблюдайте высоко methylated области, которые перекрывают инициирование каждой из двух аннотируемых расшифровок стенограммы (два средних острова CpG).
Дифференциальный анализ шаблонов метилирования
В исследовании Серра и др. также анализируется другая клеточная линия. Новые ячейки (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 нормировать данные. PERCENTILE пара "имя-значение" позволяет вам отфильтровать окна с очень большим количеством количеств при нормализации, поскольку эти окна происходят в основном из-за артефактов, таких как повторяющиеся области в ссылочной хромосоме.
counts_norm = round(manorm(counts,'percentile',90).*100);
Используйте функциональный mattest выполнять 2D демонстрационный t-тест, чтобы идентифицировать дифференцированно покрытые окна от двух различных клеточных линий.
pval = mattest(counts_norm(:,[1 2]),counts_norm(:,[3 4]),'bootstrap',true,... 'showhist',true,'showplot',true);
Создайте отчет с 25 старшими значащими дифференцированно покрытыми окнами. В то время как создание отчета использует функцию помощника 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.002668 258 257 434 428
120176801 Intergenic (ASTN2) 0.002738 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.005550 107 101 0 0
124461001 Intergenic (DAB2IP) 0.005624 77 76 39 37
140086501 Intergenic (TPRN) 0.006519 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.010344 77 76 6 12
977401 Intergenic (DMRT3) 0.010381 236 245 129 124
37002601 Intergenic (PAX5) 0.010560 133 127 207 211
139744401 Intergenic (PHPT1) 0.010884 47 46 32 31
126771301 Intragenic 0.011461 43 46 97 93
93922501 Intragenic 0.011487 34 34 149 161
94187101 Intragenic 0.011508 73 80 6 6
136044401 Intragenic 0.011566 39 34 110 105
139611201 Intergenic (FAM69B) 0.011566 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 ясно unmethylated для обоих из DICERex5, реплицирует.
Ссылки
[1] Серр, D., Ли, B.H., и Звон A.H., "ИЗОЛИРОВАННОЕ ОТ МИНИМАЛЬНОЙ ДИСФУНКЦИИ МОЗГА Секвенирование Генома предоставляет высокопроизводительный и всесторонний обзор метилирования ДНК в геноме человека", Исследование Нуклеиновых кислот, 38 (2):391-9, 2010.
[2] Langmead, B., Trapnell, C., поп, M., и Salzberg, S.L., "Сверхбыстрое и эффективное памятью выравнивание коротких последовательностей ДНК к геному человека", биология генома, 10 (3): R25, 2009.
[3] Литий, H., и др., "Sequence Alignment / карта (SAM) Формат и SAMtools", Биоинформатика, 25 (16):2078-9, 2009.
[4] Сад Гардинера, M. и Фроммер, M., "острова CpG в позвоночных геномах", Журнал Молекулярной биологии, 196 (2):261-82, 1987.
[5] Звон, A.H., и др., "Требование для DICER, чтобы Обеспечить Полного Промоутера Островное Гиперметилирование CpG в Человеческих Раковых клетках", Исследования рака, 68 (8):2570-5, 2008.
[6] Benjamini, Y. и Hochberg, Y., "Управляя ложным уровнем открытия: практический и мощный подход к нескольким тестирование", Журнал Королевского Статистического Общества, 57 (1):289-300, 1995.