Исследование связывающих участков 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.

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