Исследование различий всего генома в профилях метилирования ДНК

В этом примере показано, как выполнить анализ всего генома метилирования ДНК в человеке при помощи секвенирования генома.

Примечание: Для расширенной эффективности, 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         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

Нахождение Межгенных областей, которые значительно являются 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.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')

Наблюдайте высоко 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.

Для просмотра документации необходимо авторизоваться на сайте