Исследование связывающих участков DNA белка из данных ЧИПА-Seq Парного Конца
В этом примере показано, как выполнить анализ всего генома транскрипционного фактора в организме модели Arabidopsis Thaliana (Thale Cress).
Для расширенной производительности рекомендуется, чтобы вы запустили этот пример на 64-битной платформе, потому что объем потребляемой памяти близко к 2 Гбитам. На 32-битной платформе, если вы получаете "Out of memory" ошибки при выполнении этого примера, попытайтесь увеличить виртуальную память (или область подкачки) операционной системы или попытайтесь установить 3GB переключатель (только 32-битный Windows® XP). Эти методы описаны в этом документе.
Введение
ЧИП-Seq является технологией, которая используется, чтобы идентифицировать транскрипционные факторы, которые взаимодействуют с определенными сайтами DNA. Первая иммунопреципитация хроматина обогащает комплексы белка DNA с помощью антитела, которое связывает с конкретным белком интереса. Затем все получившиеся фрагменты обрабатываются с помощью высокопроизводительного секвенирования. Упорядочивающие фрагменты сопоставлены назад со ссылочным геномом. Путем осмотра превалирующих областей возможно отметить геномное местоположение взаимодействий белка DNA.
В этом примере короткие чтения производятся к парному концу платформа Illumina®. Каждый фрагмент восстановлен от двух коротких чтений, успешно сопоставленных с этим, точная длина фрагмента может быть вычислена. Используя информацию о парном конце от последовательности чтения максимизирует точность предсказания связывающих участков белка DNA.
Набор данных
Этот пример исследует данные ЧИПА-Seq парного конца, сгенерированные Ваном et.al. [1] использование платформы Illumina®. Набор данных был вежливо представлен репозиторию Gene Expression Omnibus с инвентарным номером GSM424618. Несопоставленные чтения парного конца могут быть получены из FTP-сайта NCBI.
Этот пример принимает что вы:
(1) загруженный данные, содержащие несопоставленное короткое чтение и преобразованный, это в FASTQ отформатировало файлы с помощью Инструментария SRA NCBI.
(2) произведенный SAM отформатировал файл путем отображения коротких чтений с геномом ссылки Кресса Thale, использования картопостроителя, таких как BWA [2], Галстук-бабочка или SSAHA2 (который является картопостроителем, используемым авторами [1] года), и,
(3) упорядоченный SAM отформатировал файл ссылочным именем сначала, затем геномным положением.
Для опубликованной версии этого примера 8 655 859 парных концов короткие чтения сопоставлены с помощью картопостроителя BWA [2]. BWA произвел отформатированный файл SAM (aratha.sam) с 17 311 718 записями (8,655,859 x 2). Повторяющиеся хиты были случайным образом выбраны, и о только одном хите сообщают, но с более низким качеством отображения. Файл SAM был упорядочен и преобразован в отформатированный файл BAM с помощью SAMtools [3] прежде, чем загрузиться в MATLAB.
Последняя часть примера также принимает, что вы загрузили ссылочный геном для организма модели 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');
Создание интерфейса MATLAB® к BAM отформатированный файл
Чтобы создать локальные выравнивания и посмотреть на покрытие, мы должны создать 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: ''
Путем доступа к положениям Запуска и Остановки сопоставленного короткого чтения можно получить геномную область значений.
x1 = min(getStart(bm1)) x2 = max(getStop(bm1))
x1 = uint32 1 x2 = uint32 30427671
Исследование покрытия в различных разрешениях
Чтобы исследовать покрытие для целой области значений хромосомы, алгоритм раскладывания требуется. getBaseCoverage метод производит сигнал покрытия на основе эффективных выравниваний. Это также позволяет вам задавать ширину интервала, чтобы управлять размером (или разрешение) выходного сигнала. Однако внутренние расчеты все еще выполняются в паре оснований (BP) разрешение. Это означает, что несмотря на установку большого размера интервала, узкий peaks в сигнале покрытия может все еще наблюдаться. Если сигнал покрытия построен, можно программировать Data Cursor фигуры, чтобы отобразить геномное положение при использовании подсказки. Можно масштабировать и панорамировать фигуру, чтобы определить положение и высоту peaks ЧИПА-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 (также называемый профилем нагромождения). Исследуйте один из большого peaks, наблюдаемого в данных в положении 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')
Большая часть большого peaks в этом наборе данных происходит из-за спутниковых повторных областей или из-за его близости с центромерой [4] и показывает характеристики, похожие на пример, только исследованный. Можно исследовать другие области с большим peaks с помощью той же процедуры.
Чтобы предотвратить эти проблематичные области, два метода используются. Во-первых, учитывая, что предоставленные данные устанавливают секвенирование парного конца использования путем удаления чтений, которые не выравниваются в соответствующей паре, сокращает количество потенциальных ошибок выравнивателя или неоднозначностей. Можно достигнуть этого путем исследования 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: ''
Визуализируйте снова отфильтрованный набор данных с помощью обоих, крупного разрешения с 1 000 интервалов BP для целой хромосомы и высоким разрешением для небольшой области 20 000 BP. Большая часть большого peaks из-за артефактов была удалена.
[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] это предполагается, что данные о секвенировании парного конца имеют потенциал, чтобы увеличиться, точность идентификации связывающих участков хромосомы DNA сопоставила белки, потому что длина фрагмента может быть выведена точно, в то время как при использовании одно конца, упорядочивающего его, необходимо обратиться к статистическому приближению длины фрагмента и использовать его неотчетливо во всех предполагаемых связывающих участках.
Используйте чтения парного конца, чтобы восстановить фрагменты секвенирования. Во-первых, получите индексы для форварда и противоположных чтений в каждой паре. Эта информация получена в пятом бите 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: ''
Исследование покрытия Используя выравнивания фрагмента
Сравните сигнал покрытия, полученный при помощи восстановленных фрагментов с сигналом покрытия, полученным при помощи отдельных чтений парного конца. Заметьте, что обогащенные связывающие участки, представленные peaks, могут быть лучше различены от фонового сигнала.
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 обычно связывают с последовательностью согласия, названной электронным полем (с CANNTG мотив). Используйте fastaread чтобы загрузить ссылочную хромосому, ищите мотив электронного поля в 3' и 5' направлениях, и затем наложите позиции мотива по сигналам покрытия. Этот пример перерабатывает область 1-200 000, однако пределы фигуры сужены в 3 000 областей BP для того, чтобы лучше изобразить детали.
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')
Заметьте, что не возможно сопоставить каждый пик в сигналах покрытия с мотивом электронного поля. Это вызвано тем, что длина фрагментов секвенирования сопоставима со средним расстоянием мотива, размывая peaks, который близок. Постройте распределение расстояний между сайтами мотива электронного поля.
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')
Нахождение значительного Peaks в сигнале покрытия
Используйте функциональный mspeaks выполнять пиковое обнаружение с шумоподавлением Вейвлетов на сигнале покрытия выравниваний фрагмента. Отфильтруйте предполагаемый peaks ChIP с помощью фильтра высоты, чтобы удалить peaks, который не обогащен процессом привязки на рассмотрении.
putative_peaks = mspeaks(bin,cov_fragments,'noiseestimator',20,... 'heightfilter',10,'showplot',true); hold on 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 функционируйте, чтобы найти самый близкий мотив к каждому предполагаемому peaks. Как ожидалось большая часть обогащенного peaks ChIP близко к мотиву электронного поля [1]. Это укрепляет важность выполнения пикового обнаружения в самом прекрасном возможном разрешении (разрешение BP), когда ожидаемая плотность связывающих участков высока, как это в случае мотива электронного поля. Этот пример также иллюстрирует, что для этого типа анализа, секвенирование парного конца должно быть рассмотрено по одно концу, упорядочивающему [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] Ван К., Сюй Цз., Чжан Д., Вильсон З.Э. и Чжан Д., "Эффективный подход для идентификации в естественных условиях связывающих участков DNA белка из данных ЧИПА-Seq парного конца", Биоинформатика BMC, 11:81, 2010.
[2] Ли Х. и Дербин Р., "Быстрое и точное короткое выравнивание чтения с Норами-Wheeler преобразовывает", Биоинформатика, 25 (14):1754-60, 2009.
[3] Литий, H., и др., "Sequence Alignment / карта (SAM) Формат и SAMtools", Биоинформатика, 25 (16):2078-9, 2009.
[4] Джоти Р., и др., "Идентификация Всего генома в естественных условиях связывающих участков DNA белка из данных ЧИПА-Seq", Исследование Нуклеиновых кислот, 36 (16):5221-31, 2008.
[5] Хуфмен Б.Г. и Джонс S.J.M., "Идентификация Всего генома взаимодействий белка DNA с помощью иммунопреципитации хроматина вместе с секвенированием ячейки потока", Журнал Эндокринологии, 201 (1):1-13, 2009.
[6] Ramsey S.A., и др., "Данные об ацетилировании гистона Всего генома улучшают прогноз связывающих участков транскрипционного фактора млекопитающих", Биоинформатика, 26 (17):2071-5, 2010.