В этом примере показано, как создать файл с отображенной памятью для данных о последовательности и работать с ними, не загружая всю геномную последовательность в память. Целые геномы доступны человеку, мыши, крысе, 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/BR2019bd_1170825_64229/publish_examples0/tp038a4618/ex57563178' date: '26-Jul-2019 08:39:07' bytes: 249250621 isdir: 0 datenum: 7.3763e+05
Команда memmapfile
создает объект memmapfile, который сопоставляет новый файл с памятью. Для того, чтобы сделать это, это должно знать формат файла. Формат этого файла прост, хотя намного более сложные форматы могут быть сопоставлены.
chr1 = memmapfile(mmFilename, 'format', 'uint8')
chr1 = Filename: '/tmp/BR2019bd_1170825_64229/publish_examples0/tp038a4618/ex57563178/Homo_sapiens.GRCh37.56.dna.chromosome.1.mm' Writable: false Offset: 0 Format: 'uint8' Repeat: Inf Data: 249250621x1 uint8 array
Объект 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);
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 '
Можно использовать 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)