Исследование сайтов связывания белок-ДНК из парных данных ChIP-Seq
Этот пример показывает, как выполнить общий для генома анализ транскрипционного фактора в организме модели Arabidopsis Thaliana (Thale Cress).
Для повышения эффективности рекомендуется запустить этот пример на 64-разрядной платформе, поскольку объем памяти близок к 2 Гб. На 32-разрядной платформе, если вы получаете "Out of memory" ошибки при запуске этого примера, попробуйте увеличить виртуальную память (или поменять пространство) вашей операционной системы или попробуйте задать 3GB switch (только 32-разрядная версия Windows ® XP). Эти методы описаны в этом документе.
Введение
ChIP-Seq является технологией, которая используется для идентификации факторов транскрипции, которые взаимодействуют с определенными сайтами ДНК. Первый иммунопреципитация хроматина обогащает комплексы ДНК-белка с помощью антитела, которое связывается с конкретным интересующим белком. Затем все полученные фрагменты обрабатывают с использованием высокопроизводительного секвенирования. Фрагменты секвенирования отображаются обратно в эталонный геном. При осмотре избыточно представленных областей можно отметить геномное расположение взаимодействий ДНК-белка.
В этом примере короткие чтения производятся парной платформой Illumina ®. Каждый фрагмент восстанавливается из двух коротких чтений, успешно отображенных, при этом точная длина фрагмента может быть вычислена. Использование информации о парном конце из показаний последовательности максимизирует точность предсказания сайтов связывания ДНК-белка.
Набор данных
В этом примере исследуются парные данные ChIP-Seq, сгенерированные Wang et.al. [1] с использованием платформы Illumina ®. Набор данных был учтиво передан в репозиторий Gene Expression Omnibus с номером GSM424618. Несопоставленные считывания парного конца могут быть получены с сайта NCBI FTP.
Этот пример принимает, что вы:
(1) загрузил данные, содержащие несопоставленные короткие считанные и преобразовал их в форматированные файлы FASTQ, используя NCBI SRA Toolkit.
(2) произвел файл с форматированием SAM путем отображения кратких показаний в эталонный геном Thale Cress с помощью отображения, такого как BWA [2], Bowtie или SSAHA2 (который является отображателем, используемым авторами [1]), и,
(3) сначала упорядочил форматированный файл SAM по ссылочному наименованию, затем по геномному положению.
Для опубликованной версии этого примера 8 655 859 коротких чтений парного конца отображаются с помощью BWA mapper [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: ''
При доступе к положениям Start и Stop отображенного краткого чтения вы можете получить геномную область значений.
x1 = min(getStart(bm1)) x2 = max(getStop(bm1))
x1 = uint32 1 x2 = uint32 30427671
Изучение охвата различными разрешениями
Чтобы исследовать покрытие для всей области значений хромосом, необходим алгоритм раскладывания. The getBaseCoverage способ формирует сигнал покрытия на основе эффективных выравниваний. Это также позволяет вам задать ширину интервала для управления размером (или разрешением) выходного сигнала. Однако внутренние расчеты все еще выполняются с разрешением базовой пары (bp). Это означает, что, несмотря на установку большого размера интервала, узкий peaks в сигнале покрытия все еще могут наблюдаться. После нанесения сигнала покрытия можно запрограммировать Data Cursor рисунка, чтобы отобразить геномное положение при использовании подсказки. Можно масштабировать и панорамировать рисунок, чтобы определить положение и высоту peaks 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 (также называемом профилем свайной). Исследуйте один из большого 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, что указывает на то, что mapper, скорее всего, обнаружил несколько попаданий в геноме ссылки, следовательно, присвоение более низкого качества отображения.
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: ''
Снова визуализируйте отфильтрованный набор данных, используя оба, грубое разрешение с 1000 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] предполагается, что данные секвенирования парного конца могут повысить точность идентификации сайтов связывания хромосом белков, сопоставленных с ДНК, поскольку длина фрагмента может быть получена точно, при использовании одноконечного секвенирования необходимо прибегнуть к статистическому приближению длины фрагмента и использовать его неясно для всех предполагаемых сайтов связывания.
Используйте парные чтения для восстановления фрагментов секвенирования. Сначала получите индексы для прямого и обратного показаний в каждой паре. Эта информация получена в пятом бите 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. В этом наборе данных Wang et.al. [1] исследовал основной фактор транскрипции спирали-петли-спирали (bHLH). Белки bHLH обычно связываются с консенсусной последовательностью, называемой E-кубом (с CANNTG мотив). Использование fastaread чтобы загрузить ссылочную хромосому, найдите мотив E-box в направлениях 3 'и 5', а затем наложите мотив на сигналы покрытия. Этот пример работает по 1-200 000 области, однако пределы рисунка сужены до области 3000 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')
Заметьте, что невозможно связать каждый пик в сигналах покрытия с мотивом E-box. Это связано с тем, что длина фрагментов секвенирования сопоставима со средним расстоянием мотива, размытием peaks, которые близки. Постройте график распределения расстояний между сайтами мотива 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')
Нахождение значительного Peaks в сигнале покрытия
Используйте функцию mspeaks для выполнения пикового обнаружения с шумоподавлением вейвлеты по сигналу покрытия выравниваний фрагментов. Фильтруйте предполагаемый peaks ChIP с помощью фильтра высоты для удаления peaks, которые не обогащаются рассматриваемым процессом связывания.
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 функция для нахождения ближайшего мотива к каждому из предполагаемого peaks. Как ожидалось, большинство обогащённого peaks ChIP близки к мотиву E-куба [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 Bioinformatics 11, № 1 (2010): 81.
[2] Li, H., and R. Durbin. Быстрая и точная выравнивание при коротком считывании с преобразованием Burrows-Wheeler. Биоинформатика 25, № 14 (15 июля 2009): 1754-60.
[3] Li, H., B. Handsaker, A. Wysoker, T. Fennell, J. Ruan, N. Homer, G. Marth, G. Abecasis, R. Durbin и 1000 Genome Project Данных Processing Subroup. «Sequence Alignment Формат/Map и SAMtools». Биоинформатика 25, № 16 (15 августа 2009): 2078-79.
[4] Jothi, R., S. Cuddapah, A. Barski, K. Cui, and K. Zhao. «Общегеномная идентификация сайтов связывания белок-ДНК in vivo из данных ChIP-Seq». Исследование нуклеиновых кислот 36, № 16 (1 августа 2008): 5221-31.
[5] Hoffman, Brad G, and Steven J M Jones. «Общегеномная идентификация взаимодействий ДНК-белка с использованием иммунопреципитации хроматина в сочетании с секвенированием проточных камер». Журнал эндокринологии 201, № 1 (апрель 2009): 1-13.
[6] Ramsey, Stephen A., Theo A. Knijnenburg, Kathleen A. Kennedy, Daniel E. Zak, Mark Gilchrist, Elizabeth S. Gold, Carrie D. Johnson, en et al al. «Ацетилирование гистонов по всему геному Данных улучшение прогноза сайтов связывания факторов транскрипции млекопитающих». Биоинформатика 26, № 17 (1 сентября 2010): 2071-75.
См. также
BioMap | getBaseCoverage | getgenbank | getSummary