Этот пример показывает, как выполнить анализ всего генома транскрипционного фактора в модельном организме Arabidopsis Thaliana (Thale Cress).
Для повышения производительности рекомендуется запустить этот пример на 64-разрядной платформе, поскольку объем памяти близок к 2 Гбайт. На 32-разрядной платформе, если вы получаете "Out of memory" ошибки при выполнении этого примера, попробуйте увеличить виртуальную память (или пространство подкачки) операционной системы или попробуйте установить 3GB коммутатор (только 32-разрядная ОС Windows ® XP). Эти методы описаны в этом документе.
ChIP-Seq - это технология, которая используется для идентификации факторов транскрипции, которые взаимодействуют со специфическими сайтами ДНК. Первая иммунопреципитация хроматина обогащает комплексы ДНК-белок с использованием антитела, которое связывается с конкретным белком, представляющим интерес. Затем все полученные фрагменты обрабатывают с использованием высокопроизводительного секвенирования. Секвенирующие фрагменты отображают обратно на эталонный геном. Проверяя чрезмерно представленные области, можно отметить геномное расположение ДНК-белковых взаимодействий.
В этом примере короткие чтения производятся парной платформой Illumina ®. Каждый фрагмент восстанавливается из двух коротких считываний, успешно отображенных, при этом точная длина фрагмента может быть вычислена. Использование информации парного конца от считывания последовательности максимизирует точность прогнозирования сайтов связывания ДНК-белка.
В этом примере рассматриваются парные данные ChIP-Seq, генерируемые Wang et.al. [1] с использованием платформы Illumina ®. Набор данных был любезно передан в хранилище Gene Expression Omnibus с регистрационным номером GSM424618. Несопоставленные парные чтения можно получить с FTP-сайта NCBI.
В этом примере предполагается, что:
(1) загрузил данные, содержащие несопоставленное короткое чтение, и преобразовал их в файлы в формате FASTQ, используя набор инструментов NCBI SRA.
(2) создал файл в формате SAM путем сопоставления коротких чтений с эталонным геномом Thale Cress, используя маппер, такой как BWA [2], Bowtie или SSAHA2 (который является маппером, используемым авторами [1]), и,
(3) упорядочил отформатированный SAM-файл сначала по имени ссылки, затем по геномной позиции.
Для опубликованной версии этого примера 8655859 спаренных коротких считываний отображаются с помощью BWA mapper [2]. BWA создала файл в формате SAM (aratha.sam) с 17 311 718 записями (8 655 859 x 2). Повторяющиеся попадания были выбраны случайным образом, и сообщается только об одном попадании, но с более низким качеством картирования. Перед загрузкой в MATLAB файл SAM был упорядочен и преобразован в файл в формате BAM с помощью SAMtools [3].
В последней части примера также предполагается, что вы загрузили эталонный геном для модельного организма Thale Cress (который включает пять хромосом). Раскомментируйте следующие строки кода для загрузки ссылки из репозитория NCBI:
% getgenbank('NC_003070','FileFormat','fasta','tofile','ach1.fasta'); % getgenbank('NC_003071','FileFormat','fasta','tofile','ach2.fasta'); % getgenbank('NC_003074','FileFormat','fasta','tofile','ach3.fasta'); % getgenbank('NC_003075','FileFormat','fasta','tofile','ach4.fasta'); % getgenbank('NC_003076','FileFormat','fasta','tofile','ach5.fasta');
Чтобы создать локальные трассы и посмотреть на покрытие, нам нужно построить BioMap. BioMap имеет интерфейс, который обеспечивает прямой доступ к отображенным коротким чтениям, хранящимся в форматированном файле BAM, таким образом минимизируя объем данных, которые фактически загружаются в рабочую область. Создать BioMap для доступа ко всем коротким чтениям, сопоставленным в форматированном файле BAM.
bm = BioMap('aratha.bam')
bm =
BioMap with properties:
SequenceDictionary: {5x1 cell}
Reference: [14637324x1 File indexed property]
Signature: [14637324x1 File indexed property]
Start: [14637324x1 File indexed property]
MappingQuality: [14637324x1 File indexed property]
Flag: [14637324x1 File indexed property]
MatePosition: [14637324x1 File indexed property]
Quality: [14637324x1 File indexed property]
Sequence: [14637324x1 File indexed property]
Header: [14637324x1 File indexed property]
NSeqs: 14637324
Name: ''
Используйте getSummary способ получения списка существующих ссылок и фактического количества коротких чтений, сопоставленных каждой из них.
getSummary(bm)
BioMap summary:
Name: ''
Container_Type: 'Data is file indexed.'
Total_Number_of_Sequences: 14637324
Number_of_References_in_Dictionary: 5
Number_of_Sequences Genomic_Range
Chr1 3151847 1 30427671
Chr2 3080417 1000 19698292
Chr3 3062917 94 23459782
Chr4 2218868 1029 18585050
Chr5 3123275 11 26975502
Остальная часть этого примера сосредоточена на анализе одной из пяти хромосом, Chr1. Создание нового BioMap для доступа к коротким чтениям, отображенным на первую хромосому, путем подстановки первой.
bm1 = getSubset(bm,'SelectReference','Chr1')
bm1 =
BioMap with properties:
SequenceDictionary: 'Chr1'
Reference: [3151847x1 File indexed property]
Signature: [3151847x1 File indexed property]
Start: [3151847x1 File indexed property]
MappingQuality: [3151847x1 File indexed property]
Flag: [3151847x1 File indexed property]
MatePosition: [3151847x1 File indexed property]
Quality: [3151847x1 File indexed property]
Sequence: [3151847x1 File indexed property]
Header: [3151847x1 File indexed property]
NSeqs: 3151847
Name: ''
Доступ к позициям Пуск (Start) и Стоп (Stop) отображенного короткого чтения позволяет получить геномный диапазон.
x1 = min(getStart(bm1)) x2 = max(getStop(bm1))
x1 = uint32 1 x2 = uint32 30427671
Чтобы исследовать охват для всего диапазона хромосомы, требуется алгоритм связывания. getBaseCoverage способ формирует сигнал покрытия на основе эффективных выравниваний. Он также позволяет задать ширину ячейки для управления размером (или разрешением) выходного сигнала. Однако внутренние вычисления все еще выполняются при разрешении базовой пары (bp). Это означает, что, несмотря на установку большого размера ячейки, узкие пики в сигнале покрытия все еще могут наблюдаться. После построения графика сигнала покрытия можно запрограммировать курсор данных рисунка для отображения положения генома при использовании всплывающей подсказки. Можно увеличить и панорамировать рисунок, чтобы определить положение и высоту пиков ChIP-Seq.
[cov,bin] = getBaseCoverage(bm1,x1,x2,'binWidth',1000,'binType','max'); figure plot(bin,cov) axis([x1,x2,0,100]) % sets the axis limits fixGenomicPositionLabels % formats tick labels and adds datacursors xlabel('Base position') ylabel('Depth') title('Coverage in Chromosome 1')
Также можно исследовать сигнал покрытия с разрешением bp (также называемый профилем свай). Изучите один из больших пиков, наблюдаемых в данных в положении 4598837.
p1 = 4598837-1000; p2 = 4598837+1000; figure plot(p1:p2,getBaseCoverage(bm1,p1,p2)) xlim([p1,p2]) % sets the x-axis limits fixGenomicPositionLabels % formats tick labels and adds datacursors xlabel('Base position') ylabel('Depth') title('Coverage in Chromosome 1')
Наблюдайте большой пик с глубиной покрытия 800 + между позициями 4599029 и 4599145. Исследуйте, как эти чтения выравниваются с эталонной хромосомой. Можно извлечь подмножество этих считываний, достаточное для удовлетворения глубины охвата 25, поскольку этого достаточно для понимания того, что происходит в этом регионе. ИспользоватьgetIndex для получения индексов к этому подмножеству. Затем использовать getCompactAlignment для отображения соответствующего множественного выравнивания коротких считываний.
i = getIndex(bm1,4599029,4599145,'depth',25); bmx = getSubset(bm1,i,'inmemory',false) getCompactAlignment(bmx,4599029,4599145)
bmx =
BioMap with properties:
SequenceDictionary: 'Chr1'
Reference: [62x1 File indexed property]
Signature: [62x1 File indexed property]
Start: [62x1 File indexed property]
MappingQuality: [62x1 File indexed property]
Flag: [62x1 File indexed property]
MatePosition: [62x1 File indexed property]
Quality: [62x1 File indexed property]
Sequence: [62x1 File indexed property]
Header: [62x1 File indexed property]
NSeqs: 62
Name: ''
ans =
35x117 char array
'AGTT AATCAAATAGAAAGCCCCGAGGGCGCCATATCCTAGGCGC AAACTATGTGATTGAATAAATCCTCCTCTATCTGTTGCGG GAGGACTCCTTCTCCTTCCCCTTTTGG'
'AGTGC TCAAATAGAAAGCCCCGAGGGCGCCATATTCTAGGAGCCC GAATAAATCCTCCTCTATCTGTTGCGGGTCGAGGACTCCT CTCCTGCCCCTTTTGG'
'AGTTCAA CCCGAGGGCGCCATATTCTAGGAGCCCAAACTATGTGATT TATCTGTTGCGGGTCGAGGACTCCTTCTCCTTCCCCTTCT '
'AGTTCAATCAAATAGAAAGC TTCTAGGAGCCCAAACTATGTGATTGAATAAATCCTCCTC AGGACTCCTTCTCCTTCCCCTTTTGG'
'AGTT AAGGAGCCCAAAATATGTGATTGAATAAATCCACCTCTAT GGACTCCTTCTCCTTCCCCTTTTGG'
'AGTACAATCAAATAGAAAGCCCCGAGGGCGCCATA TAGGAGCCCAAACTATGTGATTGAATAAATCCTCCTCTAT CCTTCACCTTCCCCTTTTGG'
'CGTACAATCAAATAGAAAGCCCCGAGGGCGCCATATTC GGAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCT TTCCCCTTTTGG'
'CGTACAATCAAATAGAAAGCCCCGAGGGCGCCATATTC GGAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCT '
'CGTACAATCAAATAGAAAGCCCCGAGGGCGCCATATTC GGAGCCCAAGCTATGTGATTGAATAAATCCTCCTCTATCT '
'CGTACAATCAAATAGAAAGCCCCGAGGGCGCCATATTC GGAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCT '
'AGTTCAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTA GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG '
'GATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTA GAGCCCAAACTATGTGATTGAATAAATCTTCCTCTATCTG '
'GATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTA GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG '
'GATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTA GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG '
'GATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTA GAGCCCAAATTATGTGATTGAATAAATCCTCCTCTATCTG '
' ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAG CCCAAACTATGTGATTGAATAAATCCTCCTCTATCTGTTG '
' ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAG CACAAACTATGTGATTGAATAAATCCTCCTCTATCTGTTG '
' ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAG CCAAACTATGTGATTGAATAAATCCTCCTCTATCTGTTGC '
' ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAG '
' ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTCG '
' ATACAATCAAATAGAAAGCCCCGGGGGCGCCATATTCTAG '
' ATTGAGTCAAATAGAAAGCCCCGAGGGCGCCATATTCTAG '
' ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAG '
' ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAG '
' ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAG '
' CAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAGGAG '
' CAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAGGAG '
' TAGGAGCCCAAACTATGTGATTGAATAAATCCTCCTCTAT '
' TAGGAGCCCAAACTATGCCATTGAATAAATCCTCCGCTAT '
' GGAGCCCAAGCTATGTGATTGAATAAATCCTCCTCTATCT '
' GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG '
' GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG '
' GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG '
' GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG '
' GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG '
В дополнение к визуальному подтверждению выравнивания, вы также можете изучить качество отображения для всех коротких чтений в этой области, так как это может указывать на потенциальную проблему. В этом случае менее одного процента коротких считываний имеют качество Phred 60, что указывает на то, что сопоставитель, скорее всего, обнаружил несколько попаданий в эталонный геном, следовательно, назначая более низкое качество отображения.
figure i = getIndex(bm1,4599029,4599145); hist(double(getMappingQuality(bm1,i))) title('Mapping Quality of the reads between 4599029 and 4599145') xlabel('Phred Quality Score') ylabel('Number of Reads')
Большинство больших пиков в этом наборе данных происходит из-за повторяющихся областей спутника или из-за его близости к центромере [4], и показывают характеристики, аналогичные только что исследованному примеру. С помощью этой же процедуры можно исследовать другие области с большими пиками.
Для предотвращения этих проблемных областей используются две методики. Во-первых, учитывая, что предоставленный набор данных использует последовательность парного конца, удаление считываний, которые не выровнены в правильной паре, уменьшает количество потенциальных ошибок выравнивания или неоднозначностей. Вы можете достичь этого, изучив flag поле форматированного SAM-файла, в котором второй менее значащий бит используется для указания, отображено ли короткое чтение в правильной паре.
i = find(bitget(getFlag(bm1),2)); bm1_filtered = getSubset(bm1,i)
bm1_filtered =
BioMap with properties:
SequenceDictionary: 'Chr1'
Reference: [3040724x1 File indexed property]
Signature: [3040724x1 File indexed property]
Start: [3040724x1 File indexed property]
MappingQuality: [3040724x1 File indexed property]
Flag: [3040724x1 File indexed property]
MatePosition: [3040724x1 File indexed property]
Quality: [3040724x1 File indexed property]
Sequence: [3040724x1 File indexed property]
Header: [3040724x1 File indexed property]
NSeqs: 3040724
Name: ''
Во-вторых, рассмотрим только уникально сопоставленные чтения. Вы можете обнаружить чтения, которые одинаково сопоставлены с различными областями ссылочной последовательности, глядя на качество отображения, потому что BWA назначает более низкое качество отображения (менее 60) для этого типа короткого чтения.
i = find(getMappingQuality(bm1_filtered)==60); bm1_filtered = getSubset(bm1_filtered,i)
bm1_filtered =
BioMap with properties:
SequenceDictionary: 'Chr1'
Reference: [2313252x1 File indexed property]
Signature: [2313252x1 File indexed property]
Start: [2313252x1 File indexed property]
MappingQuality: [2313252x1 File indexed property]
Flag: [2313252x1 File indexed property]
MatePosition: [2313252x1 File indexed property]
Quality: [2313252x1 File indexed property]
Sequence: [2313252x1 File indexed property]
Header: [2313252x1 File indexed property]
NSeqs: 2313252
Name: ''
Снова визуализируйте отфильтрованный набор данных, используя оба, грубое разрешение с 1000 bp для всей хромосомы и тонкое разрешение для небольшой области 20000 bp. Большая часть больших пиков из-за артефактов была удалена.
[cov,bin] = getBaseCoverage(bm1_filtered,x1,x2,'binWidth',1000,'binType','max'); figure plot(bin,cov) axis([x1,x2,0,100]) % sets the axis limits fixGenomicPositionLabels % formats tick labels and adds datacursors xlabel('Base Position') ylabel('Depth') title('Coverage in Chromosome 1 after Filtering') p1 = 24275801-10000; p2 = 24275801+10000; figure plot(p1:p2,getBaseCoverage(bm1_filtered,p1,p2)) xlim([p1,p2]) % sets the x-axis limits fixGenomicPositionLabels % formats tick labels and adds datacursors xlabel('Base Position') ylabel('Depth') title('Coverage in Chromosome 1 after Filtering')
В статье Вана [1] предполагается, что данные секвенирования парных концов могут повысить точность идентификации сайтов связывания хромосом ДНК-ассоциированных белков, поскольку длина фрагмента может быть получена точно. при использовании одноконцевого секвенирования необходимо прибегнуть к статистической аппроксимации длины фрагмента и использовать ее неотличимо для всех предполагаемых сайтов связывания.
Для восстановления секвенирующих фрагментов используйте парные операции чтения. Сначала получите индексы для прямого и обратного считываний в каждой паре. Эта информация захватывается в пятом бите flag в соответствии с форматом файла SAM.
fow_idx = find(~bitget(getFlag(bm1_filtered),5)); rev_idx = find(bitget(getFlag(bm1_filtered),5));
Файлы в формате SAM используют одни и те же строки заголовка для идентификации пар. Объединив строки заголовка, можно определить способ чтения краткого текста в BioMap парные. Чтобы объединить строки заголовка, просто упорядочьте их по возрастанию и используйте индексы сортировки (hf и hr) для связывания несортированных строк заголовка.
[~,hf] = sort(getHeader(bm1_filtered,fow_idx)); [~,hr] = sort(getHeader(bm1_filtered,rev_idx)); mate_idx = zeros(numel(fow_idx),1); mate_idx(hf) = rev_idx(hr);
Использовать результирующий fow_idx и mate_idx переменные для извлечения пар пар. Например, извлеките спаренные чтения для первых 10 фрагментов.
for j = 1:10 disp(getInfo(bm1_filtered, fow_idx(j))) disp(getInfo(bm1_filtered, mate_idx(j))) end
SRR054715.sra.6849385 163 20 60 40M AACCCTAAACCTCTGAATCCTTAATCCCTAAATCCCTAAA BBBBBBBBBBCBCB?2?BBBBB@7;BBC?7=7?BCC4*)3 SRR054715.sra.6849385 83 229 60 40M CCTATTTCTTGTGGTTTTCTTTCCTTCACTTAGCTATGGA 06BBBB=BBBBBBBBBBBBBBA6@@@9<*9BBA@>BBBBB SRR054715.sra.6992346 99 20 60 40M AACCCTAAACCTCTGAATCCTTAATCCCTAAATCCCTAAA =B?BCB=2;BBBBB=B8BBCCBBBBBBBC66BBB=BC8BB SRR054715.sra.6992346 147 239 60 40M GTGGTTTTCTTTCCTTCACTTAGCTATGGATGGTTTATCT BBCBB6B?B0B8B<'.BBBBBBBB=BBBBB6BBBBB;*6@ SRR054715.sra.8438570 163 47 60 40M CTAAATCCCTAAATCTTTAAATCCTACATCCATGAATCCC BC=BBBBCBB?==BBB;BB;?BBB8BCB??B-BB<*<B;B SRR054715.sra.8438570 83 274 60 40M TATCTTCATTTGTTATATTGGATACAAGCTTTGCTACGAT BBBBB=;BBBBBBBBB;6?=BBBBBBBB<*9BBB;8BBB? SRR054715.sra.1676744 163 67 60 40M ATCCTACATCCATGAATCCCTAAATACCTAATCCCCTAAA BBCB>4?+<BB6BB66BBC?77BBCBC@4ABB-BBBCCBB SRR054715.sra.1676744 83 283 60 40M TTGTTATATTGGATACAAGCTTTGCTACGATCTACATTTG CCB6BBB93<BBBB>>@B?<<?BBBBBBBBBBBBBBBBBB SRR054715.sra.6820328 163 73 60 40M CATCCATGAATCCCTAAATACCTAATTCCCTAAACCCGAA BB=08?BB?BCBBB=8BBB8?CCB-B;BBB?;;?BB8B;8 SRR054715.sra.6820328 83 267 60 40M GTTGGTGTATCTTCATTTGTTATATTGGATACGAGCTTTG BBBBB646;BB8@44BB=BBBB?C8BBBB=B6.9B8CCCB SRR054715.sra.1559757 163 103 60 40M TAAACCCGAAACCGGTTTCTCTGGTTGAAACTCATTGTGT BBBBBCBBBBBBBBBBBCBBBB?BBBB<;?*?BBBBB7,* SRR054715.sra.1559757 83 311 60 40M GATCTACATTTGGGAATGTGAGTCTCTTATTGTAACCTTA <?BBBBB?7=BBBBBBBBBBBBBB@;@>@BBBBBBBBBBB SRR054715.sra.5658991 163 103 60 40M CAAACCCGAAACCGGTTTCTCTGGTTGAAACTCATTGTGT 7?BBBBBB;=BBBB?8B;B-;BCB-B<49<6B8-BB?+?B SRR054715.sra.5658991 83 311 60 40M GATCTACATTTGGGAATGTGAGTCTCTTATTGTAACCTTA 3,<-BBCBBBBBB?=BBBBA<ABBBBBBBBB?79BBB?BB SRR054715.sra.4625439 163 143 60 40M ATATAATGATAATTTTAGCGTTTTTATGCAATTGCTTATT BBBBB@,*<8BBB++2B6B;+6B8B;8+9BB0,'9B=.=B SRR054715.sra.4625439 83 347 60 40M CTTAGTGTTGGTTTATCTCAAGAATCTTATTAATTGTTTG +BB8B0BBB?BBBB-BBBB22?BBB-BB6BB-BBBBBB?B SRR054715.sra.1007474 163 210 60 40M ATTTGAGGTCAATACAAATCCTATTTCTTGTGGTTTGCTT BBBBBBBB;.>BB6B6',BBBCBB-08BBBBB;CB9630< SRR054715.sra.1007474 83 408 60 40M TATTGTCATTCTTACTCCTTTGTGGAAATGTTTGTTCTAT BBB@AABBBCCCBBBBBBB=BBBCB8BBBBB=B6BCBB77 SRR054715.sra.7345693 99 213 60 40M TGAGGTCAATACAAATCCTATTTCTTGTGGTTTTCTTTCT B>;>BBB9,<6?@@BBBBBBBBBBBBBB7<9BBBBBB6*' SRR054715.sra.7345693 147 393 60 40M TTATTTTTGGACATTTATTGTCATTCTTACTCCTTTGGGG BB-?+?C@>9BBBBBB6.<BBB-BBB94;A4442+49';B
Использовать индексы парного конца для создания нового BioMap с минимальной информацией, необходимой для представления фрагментов секвенирования. Сначала вычислите размеры вставки.
J = getStop(bm1_filtered, fow_idx); K = getStart(bm1_filtered, mate_idx); L = K - J - 1;
Получите новую сигнатуру (или строку CIGAR) для каждого фрагмента, используя короткие прочитанные исходные сигнатуры, разделенные соответствующим количеством пропускаемых символов CIGAR (N).
n = numel(L); cigars = cell(n,1); for i = 1:n cigars{i} = sprintf('%dN' ,L(i)); end cigars = strcat( getSignature(bm1_filtered, fow_idx),... cigars,... getSignature(bm1_filtered, mate_idx));
Реконструировать последовательности для фрагментов путем конкатенации соответствующих последовательностей спаренных коротких считываний.
seqs = strcat( getSequence(bm1_filtered, fow_idx),...
getSequence(bm1_filtered, mate_idx));
Рассчитайте и постройте график распределения размера фрагмента.
J = getStart(bm1_filtered,fow_idx); K = getStop(bm1_filtered,mate_idx); L = K - J + 1; figure hist(double(L),100) title(sprintf('Fragment Size Distribution\n %d Paired-end Fragments Mapped to Chromosome 1',n)) xlabel('Fragment Size') ylabel('Count')
Создание нового BioMap для представления фрагментов секвенирования. Благодаря этому вы сможете исследовать сигналы покрытия, а также локальные выравнивания фрагментов.
bm1_fragments = BioMap('Sequence',seqs,'Signature',cigars,'Start',J)
bm1_fragments =
BioMap with properties:
SequenceDictionary: {0x1 cell}
Reference: {0x1 cell}
Signature: {1156626x1 cell}
Start: [1156626x1 uint32]
MappingQuality: [0x1 uint8]
Flag: [0x1 uint16]
MatePosition: [0x1 uint32]
Quality: {0x1 cell}
Sequence: {1156626x1 cell}
Header: {0x1 cell}
NSeqs: 1156626
Name: ''
Сравните сигнал покрытия, полученный при использовании восстановленных фрагментов, с сигналом покрытия, полученным при использовании индивидуальных парных считываний. Обратите внимание, что обогащенные сайты связывания, представленные пиками, могут быть лучше различимы от фонового сигнала.
cov_reads = getBaseCoverage(bm1_filtered,x1,x2,'binWidth',1000,'binType','max'); [cov_fragments,bin] = getBaseCoverage(bm1_fragments,x1,x2,'binWidth',1000,'binType','max'); figure plot(bin,cov_reads,bin,cov_fragments) xlim([x1,x2]) % sets the x-axis limits fixGenomicPositionLabels % formats tick labels and adds datacursors xlabel('Base position') ylabel('Depth') title('Coverage Comparison') legend('Short Reads','Fragments')
Выполните то же сравнение при разрешении bp. В этом наборе данных Ван et.al. [1] исследовали основной фактор транскрипции спираль-петля-спираль (bHLH). Белки bHLH обычно связываются с консенсусной последовательностью, называемой E-box (с CANNTG мотив). Использовать fastaread чтобы загрузить эталонную хромосому, найдите мотив E-box в направлениях 3 'и 5', а затем наложите положения мотива на сигналы покрытия. Этот пример работает над областью 1-200000, однако пределы рисунка сужены до области 3000 б.п., чтобы лучше изобразить детали.
p1 = 1; p2 = 200000; cov_reads = getBaseCoverage(bm1_filtered,p1,p2); [cov_fragments,bin] = getBaseCoverage(bm1_fragments,p1,p2); chr1 = fastaread('ach1.fasta'); mp1 = regexp(chr1.Sequence(p1:p2),'CA..TG')+3+p1; mp2 = regexp(chr1.Sequence(p1:p2),'GT..AC')+3+p1; motifs = [mp1 mp2]; figure plot(bin,cov_reads,bin,cov_fragments) hold on plot([1;1;1]*motifs,[0;max(ylim);NaN],'r') xlim([111000 114000]) % sets the x-axis limits fixGenomicPositionLabels % formats tick labels and adds datacursors xlabel('Base position') ylabel('Depth') title('Coverage Comparison') legend('Short Reads','Fragments','E-box motif')
Обратите внимание, что невозможно связать каждый пик в сигналах покрытия с мотивом E-box. Это происходит потому, что длина секвенирующих фрагментов сравнима со средним расстоянием мотива, размывая пики, которые близки. Постройте график распределения расстояний между площадками движения электронного ящика.
motif_sep = diff(sort(motifs)); figure hist(motif_sep(motif_sep<500),50) title('Distance (bp) between adjacent E-box motifs') xlabel('Distance (bp)') ylabel('Counts')
Используйте функцию mspeaks для выполнения обнаружения пика с помощью вейвлетов, деноизирующих сигнал покрытия выравниваний фрагментов. Фильтрация предполагаемых пиков ChIP с использованием фильтра высоты для удаления пиков, которые не обогащены рассматриваемым процессом связывания.
putative_peaks = mspeaks(bin,cov_fragments,'noiseestimator',20,... 'heightfilter',10,'showplot',true); hold on legend('off') plot([1;1;1]*motifs(motifs>p1 & motifs<p2),[0;max(ylim);NaN],'r') xlim([111000 114000]) % sets the x-axis limits fixGenomicPositionLabels % formats tick labels and adds datacursors legend('Coverage from Fragments','Wavelet Denoised Coverage','Putative ChIP peaks','E-box Motifs') xlabel('Base position') ylabel('Depth') title('ChIP-Seq Peak Detection')
Используйте knnsearch функция, чтобы найти ближайший мотив к каждому из предполагаемых пиков. Как и ожидалось, большинство обогащенных пиков ChIP близки к мотиву E-box [1]. Это усиливает важность выполнения обнаружения пика при максимально точном разрешении (разрешение bp), когда ожидаемая плотность сайтов связывания высока, как это происходит в случае мотива E-box. Этот пример также иллюстрирует, что для этого типа анализа парное секвенирование должно рассматриваться как секвенирование на одном конце [1].
h = knnsearch(motifs',putative_peaks(:,1)); distance = putative_peaks(:,1)-motifs(h(:))'; figure hist(distance(abs(distance)<200),50) title('Distance to Closest E-box Motif for Each Detected Peak') xlabel('Distance (bp)') ylabel('Counts')
[1] Ван, Конгмао, Цзе Сюй, Дашэн Чжан, Зои А Уилсон и Дабин Чжан. «Эффективный подход к идентификации сайтов связывания белка-ДНК in vivo из данных ChIP-Seq с парным концом». BMC Биоинформатика 11, № 1 (2010): 81.
[2] Литий, H. и Р. Дербин. «Быстрое и точное короткое выравнивание чтения с преобразованием Burrows-Wheeler». Биоинформатика 25, № 14 (15 июля 2009): 1754-60.
Li, H., B. Handsaker, A. Wysoker, T. Fennell, J. Ruan, N. Homer, G. Marth, G. Abecasis, R. Durbin и 1000 Подгруппа обработки данных проекта генома. «Формат выравнивания/сопоставления последовательностей и средства SAMtools». Биоинформатика 25, № 16 (15 августа 2009): 2078-79.
[4] Джоти, Р., С. Куддапах, А. Барски, К. Цуй и К. Чжао. «Общегеномная идентификация сайтов связывания белка-ДНК in vivo из данных ChIP-Seq». Исследование нуклеиновых кислот 36, № 16 (1 августа 2008 г.): 5221-31.
[5] Хоффман, Брэд Джи и Стивен Дж. М. Джонс. «Общегеномная идентификация взаимодействий ДНК-белка с использованием иммунопреципитации хроматина в сочетании с секвенированием проточных клеток». Журнал эндокринологии 201, № 1 (апрель 2009): 1-13.
[6] Рэмзи, Стивен А., Тео А. Книйненбург, Кэтлин А. Кеннеди, Дэниел Э. Зак, Марк Гилкрист, Элизабет С. Голд, Кэрри Д. Джонсон и др. «Данные по ацетилированию гистонов по всему геному улучшают прогнозирование сайтов связывания факторов транскрипции млекопитающих». Биоинформатика 26, № 17 (1 сентября 2010): 2071-75.
BioMap | getBaseCoverage | getgenbank | getSummary