Работа с целыми данными о геноме

В этом примере показано, как создать файл с отображенной памятью для данных о последовательности и работать с ними, не загружая всю геномную последовательность в память. Целые геномы доступны человеку, мыши, крысе, fugu, и нескольким другим организмам модели. Для многих из этих организмов одна хромосома может составить несколько сотен миллионов пар оснований долго. Работа с такими большими наборами данных может быть сложной, когда можно столкнуться с ограничениями аппаратного и программного обеспечения, которое вы используете. Этот пример показывает один способ работать вокруг этих ограничений в MATLAB®.

Проблемы обработки большого набора данных

При решении технических вычислительных задач, которые требуют, обработка и анализ больших объемов данных помещают высокий спрос на компьютерную систему. Большие наборы данных поднимают значительную память во время обработки и могут потребовать, чтобы много операций вычислили решение. Может также требоваться много времени, чтобы получить доступ к информации из больших файлов данных.

Компьютерные системы, однако, имеют ограниченную память и конечную скорость ЦП. Имеющиеся ресурсы варьируются процессором и операционной системой, последний которой также использует ресурсы. Например:

32-битные процессоры и операционные системы могут обратиться до 2^32 = 4,294,967,296 = 4 Гбайт памяти (также известный как виртуальное адресное пространство). Windows® XP и Windows® 2000 выделяют только 2 Гбайт этой виртуальной памяти к каждому процессу (таких как MATLAB). На UNIX® виртуальная память, выделенная процессу, конфигурируема системой и обычно - приблизительно 3 Гбайт. Приложение, выполняя вычисление, такое как MATLAB, может потребовать устройства хранения данных в дополнение к пользовательской задаче. Основная проблема, когда обработка больших объемов данных состоит в том, что требования к памяти программы могут превысить это доступное на платформе. Например, MATLAB генерирует "из памяти" ошибка, когда требования к данным превышают приблизительно 1,7 Гбайт на Windows XP.

Для получения дополнительной информации об управлении памятью и больших наборах данных, см. Производительность и память.

На типичной 32-битной машине максимальный размер одного набора данных, что можно работать с в MATLAB, составляет несколько сотен Мбайт, или о размере большой хромосомы. Размещение в ОЗУ файлов позволяет MATLAB работать вокруг этого ограничения и позволяет вам работать с очень большими наборами данных интуитивным способом.

Целые наборы данных генома

Последние целые наборы данных генома могут быть загружены с Веб-сайта Ensembl. Данные обеспечиваются в нескольких форматах. Они регулярно обновляются, когда новая информация о последовательности становится доступной. Этот пример будет использовать данные о ДНК человека, хранимые в формате FASTA. Хромосома 1 является (в Релизе GRCh37.56 сентября 2009) сжатым файлом на 65,6 Мбайт. После распаковки файла это - приблизительно 250 МБ. MATLAB использует 2 байта за символ, поэтому если вы считаете файл в MATLAB, будет требоваться приблизительно 500 МБ памяти.

Этот пример принимает, что вы уже загрузили и распаковали файл FASTA в свою локальную директорию. Поменяйте имя переменной FASTAfilename в подходящих случаях.

FASTAfilename = 'Homo_sapiens.GRCh37.56.dna.chromosome.1.fa';
fileInfo = dir(which(FASTAfilename))
fileInfo = 

  struct with fields:

       name: 'Homo_sapiens.GRCh37.56.dna.chromosome.1.fa'
     folder: '/mathworks/hub/qe/test_data/Bioinformatics_Toolbox/v000/demoData/biomemorymapdemo'
       date: '01-Feb-2013 11:54:41'
      bytes: 253404851
      isdir: 0
    datenum: 7.3527e+05

Файлы с отображенной памятью

Размещение в ОЗУ позволяет MATLAB получать доступ к данным в файле, как будто это находится в памяти. Можно использовать стандартные операции индексации MATLAB, чтобы получить доступ к данным. См. документацию для memmapfile для получения дополнительной информации.

Вы могли только сопоставить файл FASTA и получить доступ к данным непосредственно оттуда. Однако файл формата FASTA включает символы новой строки. memmapfile функционируйте обрабатывает эти символы таким же образом как все другие символы. При удалении их перед размещением в ОЗУ файл сделает операции индексации более простыми. Кроме того, размещение в ОЗУ не работает непосредственно с символьными данными, таким образом, необходимо будет обработать данные как 8-битные целые числа (uint8 класс). Функциональный nt2int в Bioinformatics Toolbox™ может использоваться, чтобы преобразовать символьную информацию в целочисленные значения. int2nt используется, чтобы преобразовать назад в символы.

Откройте сначала файл FASTA и извлеките заголовок.

fidIn = fopen(FASTAfilename,'r');
header = fgetl(fidIn)
header =

    '>1 dna:chromosome chromosome:GRCh37:1:1:249250621:1'

Откройте файл, чтобы быть с отображенной памятью.

[fullPath, filename, extension] = fileparts(FASTAfilename);
mmFilename = [filename '.mm']
fidOut = fopen(mmFilename,'w');
mmFilename =

    'Homo_sapiens.GRCh37.56.dna.chromosome.1.mm'

Считайте файл FASTA в блоках 1 МБ, удалите символы новой строки, преобразуйте в uint8 и запишите в файл MM.

newLine = sprintf('\n');
blockSize = 2^20;
while ~feof(fidIn)
    % Read in the data
    charData = fread(fidIn,blockSize,'*char')';
    % Remove new lines
    charData = strrep(charData,newLine,'');
    % Convert to integers
    intData = nt2int(charData);
    % Write to the new file
    fwrite(fidOut,intData,'uint8');
end

Закройте файлы.

fclose(fidIn);
fclose(fidOut);

Новый файл о том же размере как старый файл, но не содержит новые строки или информацию о заголовке.

mmfileInfo = dir(mmFilename)
mmfileInfo = 

  struct with fields:

       name: 'Homo_sapiens.GRCh37.56.dna.chromosome.1.mm'
     folder: '/tmp/BR2021bd_1724986_5107/publish_examples1/tpf4666575/ex57563178'
       date: '22-Jul-2021 07:33:09'
      bytes: 249250621
      isdir: 0
    datenum: 7.3836e+05

Доступ к данным в файле с отображенной памятью

Команда memmapfile создает объект memmapfile, который сопоставляет новый файл с памятью. Для того, чтобы сделать это, это должно знать формат файла. Формат этого файла прост, хотя намного более сложные форматы могут быть сопоставлены.

chr1 = memmapfile(mmFilename, 'format', 'uint8')
chr1 = 

    Filename: '/tmp/BR2021bd_1724986_5107/publish_examples1/tpf4666575/ex57563178/Homo_sapiens.GRCh37.56.dna.chromosome.1.mm'
    Writable: false
      Offset: 0
      Format: 'uint8'
      Repeat: Inf
        Data: 249250621x1 uint8 array

Объект MEMMAPFILE

Объект memmapfile имеет различные свойства. Filename хранит полный путь к файлу. Writable указывает, могут ли данные быть изменены. Обратите внимание на то, что, если вы действительно измените данные, это также изменит исходный файл. Offset позволяет вам задавать пространство, использованное любой информацией о заголовке. Format указывает на формат данных. Repeat используется, чтобы задать сколько блоков (как задано Format) сопоставлять. Это может быть полезно для ограничения, сколько памяти используется, чтобы создать карту распределения памяти. К этим свойствам можно получить доступ таким же образом как другие данные MATLAB. Для получения дополнительной информации смотрите, вводят help memmapfile или doc memmapfile.

chr1.Data(1:10)
ans =

  10x1 uint8 column vector

   15
   15
   15
   15
   15
   15
   15
   15
   15
   15

Можно получить доступ к любой области операций индексации использования данных.

chr1.Data(10000000:10000010)'
ans =

  1x11 uint8 row vector

   1   1   2   2   2   2   3   4   2   4   2

Помните, что информация о нуклеотиде была преобразована в целые числа. Можно использовать int2nt вернуть информацию о последовательности.

int2nt(chr1.Data(10000000:10000010)')
ans =

    'AACCCCGTCTC'

Или используйте seqdisp отобразить последовательность.

seqdisp(chr1.Data(10000000:10001000)')
ans =

  17x71 char array

    '   1  AACCCCGTCT CTACAATAAA TTAAAATATT AGCTGGGCAT GGTGGTGTGT GCTTGTAGTC'
    '  61  CCAGCTACTT GGCGGGCTGA GGTGGGAGAA TCATCCAAGC CTTGGAGGCA GAGGTTGCAG'
    ' 121  TGAGCTGAGA TTGTGACACT GCACTCCAGC CTGGGAGACA GAGTGAGACT CCTACTCAAA'
    ' 181  AAAAAACAAA AAACAAAAAA CAAACCACAA AACTTTCCAG GTAACTTATT AAAACATGTT'
    ' 241  TTTTGTTTGT TTTGAGACAG AGTCTTGCTC TGTCGCCCAG GCTGGAGTGC AGTGGAGCAA'
    ' 301  TCTCAGCTCA CTGCAAGCTC CGCCTCCCGG GTTCACACCA TTCTCCTGCC TCAGCCTCCC'
    ' 361  GAGTAGCTAG GACTATAGGC ACCCGCCACC ACGCCCAGCT TATTTTTTTT GTATTTTTTA'
    ' 421  GTAGAGACGG GGTTTCATCG TGTTAGCCAG GATGGTCTCG ATCTCCTGAC CTCGTGATCC'
    ' 481  GCCCACCTCA GCCTCCCAAA GTGCTGGGAT TACAGGCGTG AGCCACTGCA CCCGGCCTAG'
    ' 541  TTTTTGTATA TTTTTTTTAG TAGAGACAGG GTTTCACCAT GTTAGCCAGG ATGGTCTCAA'
    ' 601  TCTCCTGACC TCGTGATCCG CCCGCCTCGG CCTCCCAAAG TGCTGGGGTT ACAGGCGTGA'
    ' 661  GCCACCGCAC ACAGCATTAA AGCATGTTTT ATTTTCCTAC ACATAATGAA ATCATTACCA'
    ' 721  GATGATTTGA CATGTGTACT TCATTGGAGA GGATTCTTAC AGTATATTCA AAATTAAATA'
    ' 781  TAATGACAAA AAATTACTAC CTAATCTATT AAAATTGGCA TAAGTCATCT ATGATCATTA'
    ' 841  ATGATATGCA AACATAAACA AGTATTATAC CCAGAAGTGT AATTTATTGT AGCTACATCT'
    ' 901  TATGTATAAT AGTTTAGTGG ATTTTTCCTG GAAATTGTCC ATTTTAATTT TTCTCTTAAG'
    ' 961  TCTGTGGAAT TTTCCAGTAA AAGTCAAGGC AAACCCAAGA T                    '

Анализ целой хромосомы

Теперь, когда можно легко получить доступ к целой хромосоме, можно анализировать данные. Этот пример показывает один способ посмотреть на содержимое GC вдоль хромосомы.

Вы извлекаете блоки 500000bp и вычисляете содержимое GC.

Вычислите сколько блоков, чтобы использовать.

numNT = numel(chr1.Data);
blockSize = 500000;
numBlocks = floor(numNT/blockSize);

Один способ помочь эффективности MATLAB, когда работа с большими наборами данных должна "предварительно выделить" место для данных. Это позволяет MATLAB выделять достаточно места для всех данных вместо того, чтобы иметь необходимость вырастить массив в маленьких фрагментах. Это ускорит вещи и также защитит вас от проблем данных, становящихся слишком большим, чтобы сохранить. Для получения дополнительной информации о предварительном выделении массивов см.: http://www.mathworks.com/support/solutions/data/1-18150.html? solution=1-18150

Простой способ предварительно выделить массив состоит в том, чтобы использовать zeros функция.

ratio = zeros(numBlocks+1,1);

Цикл по данным, ища C или G и затем делит этот номер на общее количество A, T, C, и G. Это займет приблизительно минуту, чтобы запуститься.

A = nt2int('A'); C = nt2int('C'); G = nt2int('G'); T = nt2int('T');

for count = 1:numBlocks
    % calculate the indices for the block
    start = 1 + blockSize*(count-1);
    stop = blockSize*count;
    % extract the block
    block = chr1.Data(start:stop);
    % find the GC and AT content
    gc = (sum(block == G | block == C));
    at = (sum(block == A | block == T));
    % calculate the ratio of GC to the total known nucleotides
    ratio(count) = gc/(gc+at);
end

Итоговый блок меньше, так обработайте это как особый случай.

block = chr1.Data(stop+1:end);
gc = (sum(block == G | block == C));
at = (sum(block == A | block == T));
ratio(end) = gc/(gc+at);

График содержимого GC для хромосомы человека разумного 1

xAxis = [1:blockSize:numBlocks*blockSize, numNT];
plot(xAxis,ratio)
xlabel('Base pairs');
ylabel('Relative GC content');
title('Relative GC content of Homo Sapiens Chromosome 1')

Область в центре графика вокруг 140Mbp является большой областью нс.

seqdisp(chr1.Data(140000000:140001000))
ans =

  17x71 char array

    '   1  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    '  61  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    ' 121  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    ' 181  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    ' 241  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    ' 301  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    ' 361  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    ' 421  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    ' 481  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    ' 541  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    ' 601  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    ' 661  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    ' 721  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    ' 781  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    ' 841  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    ' 901  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN'
    ' 961  NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN N                    '

Нахождение областей высокого содержимого GC

Можно использовать find идентифицировать области высокого содержимого GC.

indices = find(ratio>0.5);
ranges = [(1 + blockSize*(indices-1)), blockSize*indices];
fprintf('Region %d:%d has GC content %f\n',[ranges ,ratio(indices)]')
Region 500001:1000000 has GC content 0.501412
Region 1000001:1500000 has GC content 0.598332
Region 1500001:2000000 has GC content 0.539498
Region 2000001:2500000 has GC content 0.594508
Region 2500001:3000000 has GC content 0.568620
Region 3000001:3500000 has GC content 0.584572
Region 3500001:4000000 has GC content 0.548137
Region 6000001:6500000 has GC content 0.545072
Region 9000001:9500000 has GC content 0.506692
Region 10500001:11000000 has GC content 0.511386
Region 11500001:12000000 has GC content 0.519874
Region 16000001:16500000 has GC content 0.513082
Region 17500001:18000000 has GC content 0.513392
Region 21500001:22000000 has GC content 0.505598
Region 156000001:156500000 has GC content 0.504446
Region 156500001:157000000 has GC content 0.504090
Region 201000001:201500000 has GC content 0.502976
Region 228000001:228500000 has GC content 0.511946

Если вы хотите удалить временный файл, необходимо сначала очистить memmapfile объект.

clear chr1
delete(mmFilename)