Этот пример показывает, как выполнить общегеномный анализ метилирования ДНК у человека с использованием секвенирования генома.
Примечание.Для повышения производительности MathWorks рекомендует запустить этот пример на 64-разрядной платформе, поскольку объем памяти близок к 2 ГБ. На 32-разрядной платформе, если вы получаете "Out of memory" ошибки при выполнении этого примера, попробуйте увеличить виртуальную память (или пространство подкачки) операционной системы или попробуйте установить 3GB коммутатор (только 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 можно получить несопоставленные односторонние чтения для четырех экспериментов секвенирования. Короткие чтения производились с использованием анализатора генома Illumina ® II. Средний размер вставки составляет 120 п.о., а длина коротких чтений - 36 п.о.
В этом примере предполагается, что:
(1) загрузил файлы SRR030222.sra, SRR030223.sra, SRR030224.sra и SRR030225.sra содержит несопоставленные короткие чтения для двух реплик из DICERex5 выборки и двух реплик из HCT116 выборки, соответственно, из селектора выполнения NCBI SRA и преобразовал их в файлы в формате FASTQ с помощью набора инструментов NCBI SRA.
(2) создание файлов в формате SAM путем отображения коротких чтений в эталонный геном человека (NCBI Build 37.5) с использованием алгоритма Bowtie [2]. Сообщается только об уникально сопоставленных операциях чтения.
(3) сжал форматированные SAM файлы в BAM и упорядочил их по имени ссылки, затем по геномной позиции с помощью SAMtools [3].
В этом примере также предполагается, что загружен эталонный геном человека (GRCh37.p5). Вы можете использовать bowtie-inspect команда для восстановления ссылки человека непосредственно из индексов bowtie. Либо можно загрузить ссылку из репозитория 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 с помощью 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 п.н. Это область 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 п.о. статистически смоделировано.
Сначала используйте getCounts метод для подсчета количества сопоставленных операций чтения, которые начинаются в каждом окне. В этом примере используется подход binning, который учитывает только начальную позицию каждого отображенного чтения, следуя подходу 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 для вычисления скорректированных значений p БУР используется процедура линейного повышения (LSU) ([6]). Установка FDR < 0,01 позволяет идентифицировать окна 100 bp, которые значительно метилированы.
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 таблицу tab-separed-value (TSV) со всеми генами, кодирующими белок, в текстовый файл. ensemblmart_genes_hum37.txt. В этом примере мы используем Ensembl версии 64. Используя сервис Ensembl's 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 имеют значительно метилированные участки ДНК (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 'промоторных областей генов. Используя подобный подход, вы можете найти гены, которые имеют межгенные метилированные области. Чтобы компенсировать изменяющиеся длины генов, можно использовать максимальный охват, вычисленный base-by-base, вместо исходного числа отображенных коротких чтений. Другим альтернативным подходом к нормализации подсчета по длине гена является установка 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 для нормализации данных. 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 наиболее значимыми различными окнами. При создании отчета используйте функцию помощника 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. и Ting A.H., «MBD-изолированное секвенирование генома обеспечивает высокопроизводительное и комплексное исследование метилирования ДНК в геноме человека», Nucleic Acids Research, 38 (2): 391-9, 2010.
[2] Langmead, B., Trapnell, C., Pop, M. и Salzberg, S.L., «Ультраскоростное и эффективное с точки зрения памяти согласование коротких последовательностей ДНК с геномом человека», Genome Biology, 10 (3): R25, 2009.
Li, H., et al., «The Sequence Alignment/Map (SAM) Format and SAMtools», Bioinformatics, 25 (16): 2078-9, 2009.
Gardiner-Garden, M. and Frommer, M., «CpG островки в геномах позвоночных», Journal of Molecular Biology, 196 (2): 261-82, 1987.
Ting, A.H., et al., «Требование к DICER для поддержания полного промоторного гиперметилирования острова CpG в раковых клетках человека», Cancer Research, 68 (8): 2570-5, 2008.
[6] Benjamini, Y. and Hochberg, Y., «Контроль частоты ложных открытий: практический и мощный подход к множественному тестированию», Journal of the Royal Statistical Society, 57 (1): 289-300, 1995.