Этот пример показывает, как создать файл с отображенной памятью для данных о последовательности и работать с ними, не загружая всю геномную последовательность в память. Целые геномы доступны человеку, мыши, крысе, 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/BR2019ad_1040972_227981/publish_examples0/tp0acc86fc/ex57563178' date: '17-Jan-2019 22:44:39' bytes: 249250621 isdir: 0 datenum: 7.3744e+05
Команда memmapfile
создает memmapfile, возражает, что сопоставляет новый файл с памятью. В порядке сделать это, это должно знать формат файла. Формат этого файла прост, хотя намного более сложные форматы могут быть сопоставлены.
chr1 = memmapfile(mmFilename, 'format', 'uint8')
chr1 = Filename: '/tmp/BR2019ad_1040972_227981/publish_examples0/tp0acc86fc/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)