Этот пример иллюстрирует простой метагеномный анализ набора выборочных данных из Саргассового моря. Это запрашивает информацию таксономии, включенную в файлы gi_taxid_prot.dmp, names.dmp и nodes.dmp (см. сжатый файл taxdump), который можно загрузить с FTP-сайта таксономии NCBI.
Метагеномика является исследованием таксономического состава выборки организмов, полученных из общей среды обитания. Это обычно состоит из сравнения выборок последовательности против баз данных известных последовательностей и использования информации о таксономии, чтобы классифицировать демонстрационные разновидности. Главные цели метагеномного анализа включают квантификацию относительной распространенности известных разновидностей и идентификацию неизвестных последовательностей, для которых еще не были идентифицированы никакие родственники.
В этом примере мы рассматриваем небольшое подмножество (100 чтений) набора данных Саргассового моря [1], который искался против базы данных NCBI-NR с помощью BLASTX с параметрами по умолчанию. Для удобства получившийся отчет BLAST был сохранен и сжат в файл sargasso-sample1-100.rpt.gz, и этому предоставляют Bioinformatics Toolbox™. Мы читаем содержимое отчета и извлекаем релевантную информацию, такую как высоко выигрывающие пары, их счет, значение ожидания и идентичность процента.
% === open the blastx report reportFilename = gunzip('sargasso-sample1-100.rpt.gz',tempdir); fid = fopen(reportFilename{1}, 'rt'); % === read all strings to be able to write into xls blastInfo = textscan(fid, '%s %s %s %s %s %s %s %s %s %s %s %s'); fclose(fid); delete(reportFilename{1}); % === extract relevant information queries = blastInfo{1}; hits = blastInfo{2}; ident = str2double(blastInfo{3}); evalue = str2double(blastInfo{11}); score = str2double(blastInfo{12}); numEntries = numel(queries)
numEntries =
19817
Поскольку мы интересуемся только значительными хитами, мы фильтруем результаты на основе их счета, значения ожидания и идентичности процента с последовательностями запроса. При помощи этого процесса фильтрации мы сокращаем количество хитов приблизительно к одной четверти исходных хитов.
% === setup filter criteria scoreThreshold = 100; evalueThreshold = 10^-5; identThreshold = 50; % === consider only hits satisfying the criteria k = find(score > scoreThreshold & evalue < evalueThreshold & ident > identThreshold); queries = queries(k); hits = hits(k); evalue = evalue(k); score = score(k); numEntries = length(k) % === clear report clear blastInfo
numEntries =
5252
Таксономические классификации для всех последовательностей GenBank® хранятся в больших файлах, которые обновляются еженедельно, когда новые последовательности представлены, и таксономическая информация усовершенствована. Чтобы получить эту информацию быстрым и эффективным способом, мы создаем карту между любым возможным gi номером в базе данных GenBank и ее связанным таксономическим идентификатором (taxid). Поскольку в настоящее время существует больше чем 100 миллионов живых gi чисел, требования к памяти для загрузки такого большого набора данных могут быть очень требовательными. Таким образом использование предоставленного помощника функционирует mapTaxoFile, мы считываем данные в блоках 1 МБ, сохраняем их как двоичный файл и затем используем функциональный memmapfile, чтобы сопоставить в память содержимое самого файла, так, чтобы к данным можно было получить доступ с помощью стандартных операций индексации. Дополнительную информацию см. в справке memmapfile.
taxoFilenameIn = 'gi_taxid_prot.dmp'; taxoFilenameOut = 'gi_taxid_prot_map.dmp'; % === create map so that gi --> taxid, taxid = -1 if no live gi blockSize = 2^20; % block size (1MB) mapTaxoFile(taxoFilenameIn, taxoFilenameOut, blockSize); % === map file into memory mt = memmapfile(taxoFilenameOut, 'format', 'int32');
Мы можем получить доступ к taxid первых десяти живых последовательностей GenBank можно следующим образом:
q = find(mt.Data(1:100)>0);
mt.Data(q(1:10))
clear q
ans = 10x1 int32 column vector 9913 9913 9913 9913 9913 9913 9913 9913 9913 9913
Мы теперь интересуемся выполнением таксономической аннотации каждого хита в отчете BLAST. Мы извлекаем gi количество каждого хита и получаем его связанный taxid.
% === extract gi number for each hit gi = zeros(1, numEntries); for i = 1:numEntries g = str2double(regexpi(hits{i}, '(?<=gi\|)\d+', 'match', 'once')); if ~isempty(g) gi(i) = g; end end % === determine taxid for each hit taxid = mt.Data(gi);
Если вы выполнили поиск BLAST против базы данных, которая устарела относительно информации о таксономии, включенной в файл nodes.dmp, некоторые gi числа могут быть заменены. Поэтому необходимо исключить из анализа те последовательности, сопоставленные с замененными записями.
% === ignore dead gi numbers
livegi = (taxid > 0);
gi = gi(livegi);
taxid = taxid(livegi);
queries = queries(livegi);
hits = hits(livegi);
evalue = evalue(livegi);
score = score(livegi);
Во время поиска против Базы данных NCBI-NR первый запрос (SHAA001TR) поразил последовательности n значительным значением ожидания и счетом. Мы можем посмотреть на таксономическое присвоение этих хитов с помощью массива taxid.
SHAA001TR = strcmp('SHAA001TR', queries);
n = sum(SHAA001TR)
hits(SHAA001TR)
taxid(SHAA001TR)
n =
12
ans =
12x1 cell array
{'gi|118591585|ref|ZP_01548982.1|'}
{'gi|83951381|ref|ZP_00960113.1|' }
{'gi|86137830|ref|ZP_01056406.1|' }
{'gi|149203209|ref|ZP_01880179.1|'}
{'gi|114769111|ref|ZP_01446737.1|'}
{'gi|56709160|ref|YP_165205.1|' }
{'gi|85704868|ref|ZP_01035969.1|' }
{'gi|110681001|ref|YP_684008.1|' }
{'gi|121611410|ref|YP_999217.1|' }
{'gi|99080687|ref|YP_612841.1|' }
{'gi|84514612|ref|ZP_01001976.1|' }
{'gi|87119306|ref|ZP_01075204.1|' }
ans =
12x1 int32 column vector
384765
89187
314262
391613
367336
246200
314264
375451
391735
292414
314232
314277
Каждый taxid соответствует определенному таксону, которому дали научное название и возможно различные синонимы. В наших целях классификации мы интересуемся научными названиями только. Таким образом мы извлекаем эту информацию и аннотируем каждый BLAST, пораженный в отчете с помощью научных названий, а не taxids.
% === read taxonomy name file taxonomyFilenameIn = 'names.dmp'; fid1 = fopen(taxonomyFilenameIn,'rt'); nameInfo = textscan(fid1, '%d%s%s%s', 'delimiter', '|'); fclose(fid1); % === preallocate space for SN maxTaxid = max(double(nameInfo{1})); SN = repmat({''}, maxTaxid, 1); % === populate array so that taxid --> scientific name ind = strncmp('scientific',nameInfo{4},10); % indices of scientific names in the array SN(nameInfo{1}(ind)) = strtrim(nameInfo{2}(ind)); % === assign name to every hit sciNames = SN(taxid);
Мы можем посмотреть на научные названия организмов, последовательности которых были поражены первым запросом путем считания первых элементов n в массиве sciNames, можно следующим образом:
sciNames(1:n)
ans =
12x1 cell array
{'Labrenzia aggregata IAM 12614' }
{'Roseovarius nubinhibens ISM' }
{'Roseobacter sp. MED193' }
{'Roseovarius sp. TM1035' }
{'Rhodobacterales bacterium HTCC2255'}
{'Ruegeria pomeroyi DSS-3' }
{'Roseovarius sp. 217' }
{'Roseobacter denitrificans OCh 114' }
{'Verminephrobacter eiseniae EF01-2' }
{'Ruegeria sp. TM1040' }
{'Loktanella vestfoldensis SKA53' }
{'Marinomonas sp. MED121' }
Если мы определяем таксономическую классификацию для каждого хита, мы можем включать информацию в текстовый файл как показано ниже:
% === create annotated report for first n hits textFilename = 'sargasso-annotated-report.txt'; fid = fopen(textFilename, 'wt'); for i = 1:n fprintf(fid, '%s\t%s\t%d\t%d\t%s\n', queries{i}, hits{i}, evalue(i), taxid(i), sciNames{i}); end fclose(fid); type sargasso-annotated-report.txt
SHAA001TR gi|118591585|ref|ZP_01548982.1| 2.000000e-90 384765 Labrenzia aggregata IAM 12614 SHAA001TR gi|83951381|ref|ZP_00960113.1| 5.000000e-89 89187 Roseovarius nubinhibens ISM SHAA001TR gi|86137830|ref|ZP_01056406.1| 4.000000e-87 314262 Roseobacter sp. MED193 SHAA001TR gi|149203209|ref|ZP_01880179.1| 5.000000e-87 391613 Roseovarius sp. TM1035 SHAA001TR gi|114769111|ref|ZP_01446737.1| 8.000000e-87 367336 Rhodobacterales bacterium HTCC2255 SHAA001TR gi|56709160|ref|YP_165205.1| 1.000000e-86 246200 Ruegeria pomeroyi DSS-3 SHAA001TR gi|85704868|ref|ZP_01035969.1| 4.000000e-86 314264 Roseovarius sp. 217 SHAA001TR gi|110681001|ref|YP_684008.1| 3.000000e-84 375451 Roseobacter denitrificans OCh 114 SHAA001TR gi|121611410|ref|YP_999217.1| 4.000000e-83 391735 Verminephrobacter eiseniae EF01-2 SHAA001TR gi|99080687|ref|YP_612841.1| 4.000000e-83 292414 Ruegeria sp. TM1040 SHAA001TR gi|84514612|ref|ZP_01001976.1| 2.000000e-80 314232 Loktanella vestfoldensis SKA53 SHAA001TR gi|87119306|ref|ZP_01075204.1| 2.000000e-79 314277 Marinomonas sp. MED121
Одна причина классифицировать хиты последовательности на отчет BLAST состоит в том, чтобы изучить их таксономическое распределение. Мы можем легко создать список организмов, которые представлены в отчете, их taxids и их частоте можно следующим образом:
% === distribution by taxid taxidList = unique(taxid); % list of unique taxids T = accumarray(taxid, 1); % multiplicity of taxids taxidCount = T(unique(taxid)); % number of hits for each taxon % === simple statistics of the hit distribution numTaxa = length(taxidList) % number of distinct taxa [maxCount,maxInd] = max(taxidCount); % most represented taxon maxTaxid = taxidList(maxInd) % taxid of the most represented taxon maxSN = SN(maxTaxid) % name of the most represented taxon maxCount
numTaxa =
834
maxTaxid =
int32
269483
maxSN =
1x1 cell array
{'Burkholderia sp. 383'}
maxCount =
45
От простой статистики по таксономическому распределению мы замечаем, что наиболее представленный таксон является SP 383 Burkholderia (taxid 269483). Перепредставленность этой бактерии, которая обычно находится в наземных настройках, в Демонстрационном 1 из набора данных Саргассового моря обсуждена в [1].
Несколько таксонов в отчете, кажется, изолированные присвоения, потому что они поражены только одной последовательностью. Эти таксоны являются редко истинными членами экологического сообщества под следствием. Таким образом полезно идентифицировать их и отбросить их в случае необходимости.
t1 = taxidCount == 1; isolated = length(find(t1)) taxidList = taxidList(~t1); taxidCount = taxidCount(~t1); numTaxaFiltered = length(taxidCount)
isolated = 298 numTaxaFiltered = 536
Если мы строим таксономическое распределение хитов на столбчатой диаграмме, мы замечаем, что у большинства таксонов есть небольшое число случаев.
% === plot by sorting the counts hFig = figure(); bar(sort(taxidCount)); xlabel('Distinct taxonomic assignments'); ylabel('Number of hits'); title('Taxonomic distribution of filtered hits'); ax = gca; ax.XTickLabel = '';
Мы можем повторить вышеупомянутую процедуру путем ограничения анализа только лучше всего выигрывающим хитом для каждой последовательности запроса. Даже при том, что исследования, ограниченные лучше всего выигрывающими хитами, не могут изобразить полное и точную картину ситуации, они могут быть полезными как первое приближение и преодолеть трудность, свойственную с большими наборами данных.
% == get only best hits [queriesUnique, idx] = unique(queries, 'first'); % best hits rows bestHitTaxid = taxid(idx); bestHitSciName = sciNames(idx); % === count occurrences T = accumarray(bestHitTaxid, 1); % multiplicity of taxids bestCount = T(unique(bestHitTaxid)); % number of hits for each taxon bestCountNames = SN(unique(bestHitTaxid)); % === five most represented taxa [bestCountSorted, idx] = sort(bestCount, 'descend'); bestCountSorted(1:5) bestCountNames(idx(1:5))
ans =
16
8
7
4
3
ans =
5x1 cell array
{'Burkholderia sp. 383' }
{'Candidatus Pelagibacter ubique HTCC1062'}
{'Candidatus Pelagibacter ubique HTCC1002'}
{'Shewanella sp. ANA-3' }
{'Shewanella sp. MR-7' }
В нашем примере, когда только лучше всего выигрывающие хиты рассматриваются, Burkholderia, Candidatus pelagibacter ubique и Shewanella, кажется, наиболее представленные таксоны в отчете. При нахождении Candidatus pelagibacter ubique не удивителен, потому что это - доминирующая форма жизни в Саргассовом море, Burkholderia и Shewanella, как ожидают, не будут присутствовать в этой морской выборке, где питательные вещества и ресурсы являются низкими, потому что они живут или в наземных настройках или в водных, богатых питательным веществом средах соответственно. Для детального обсуждения относительно присутствия этих бактерий в Саргассовом море см. [1].
Часто, чтобы получить четкое представление таксономического распределения набора последовательности, категории Linnaean выше, чем разновидности рассматриваются. Чтобы выполнить этот анализ, мы должны создать карту между каждым taxid и его присвоенным рангом, а также карту между каждым taxid и taxid его родительского узла, согласно Схеме базы данных Таксономии NCBI. Файлы, содержащие эту информацию, могут быть созданы с функцией помощника mapNodeFile.
nodeFilename = 'nodes.dmp'; parentFilename = 'nodes_parent_map.dmp'; rankFilename = 'nodes_rank_map.dmp'; % === create a map mapNodeFile(nodeFilename, parentFilename, rankFilename, blockSize); % === map the files into memory mmParentObj = memmapfile(parentFilename, 'format', 'int32'); % taxid --> taxid_parent mmRankObj = memmapfile(rankFilename, 'format', 'int32'); % taxid --> rank
После того, как карты создаются для каждого хита, который сопоставлен с соответствием taxid категории Linnaean, более конкретной, чем целевой ранг, мы определяем родительский taxid и его ранг, пока цель не достигнута. Затем мы аннотируем хит taxid его большего количества удаленного предка. Синтетические построения или узлы без ранга рассматриваются потомками корня. Эта процедура хождения по таксономической иерархии выполняется функцией помощника findTaxoRank.
Предположим, что мы интересуемся классификацией наших хитов согласно суперкоролевству, которому они принадлежат. После присвоения суперкоролевства taxid к каждому хиту, мы группируем и считаем случаи можно следующим образом:
% === find superkingdom assignments skRank = findTaxoRank(taxidList, mmRankObj, mmParentObj, 1); sk = accumarray(skRank, 1); skCount = sk(unique(skRank)); skNames = SN(unique(skRank)); % === plot pie chart hFig = figure(); pie(skCount); colormap(summer) legend(skNames, 'location', 'EastOutside');
Как ожидалось большинство хитов является бактериями. Точно так же мы можем определить таксономическое распределение на уровне филюма, класса, порядка и семейства как показано ниже:
rTargetString = {'phylum', 'class', 'order', 'family'}
rTarget = [5 8 11 14];
numTarget = numel(rTarget);
rank = cell(1,numTarget);
% === annotate hits with the taxid at the target level
for i = 1:numTarget
rank{i} = findTaxoRank(taxidList, mmRankObj, mmParentObj, rTarget(i));
end
% === determine the distribution
count = cell(1,numTarget);
names = cell(1,numTarget);
for i = 1:numTarget
list = unique(rank{i});
T = accumarray(rank{i}, 1);
count{i} = T(list);
names{i} = SN(list);
end
% === plot the first two classifications
for i = 1:2
figure();
barh(count{i});
ax = gca;
ax.YTick = 1:numel(names{i});
ax.YTickLabel = names{i};
xlabel('Occurrences');
title(['Taxonomic distribution at the ' rTargetString{i} ' level'])
end
% === Draw a Pareto chart for the phyla
pnames = names{1};
pcount = count{1};
np = numel(pnames);
[ppeaks, pind] = sort(pcount, 'descend');
plabels = pnames(pind);
figure();
pareto(pcount, pnames);
ylabel('Occurrences');
text(1:numel(ppeaks), ppeaks+10, plabels, 'rotation', 90, 'clipping', 'on');
title('Pareto chart for distribution at the phylum level');
ax = gca;
ax.XTickLabel = '';
rTargetString =
1x4 cell array
{'phylum'} {'class'} {'order'} {'family'}
Таксономические дистрибутивы на разных уровнях связаны друг с другом иерархией. Предположим, что мы хотим посмотреть на распределение хитов через филюмы и визуализировать их на графике. После отфильтровывания количеств низких представленных филюмов (<5 количеств), мы создаем матрицу смежности, где все филюмы являются прямыми дочерними элементами корня.
k = count{1} > 5;
phylaNames = names{1}(k);
n1 = length(phylaNames);
CM = zeros(n1);
CM(1,2:end) = 1;
bg = biograph(CM, phylaNames);
view(bg)
Мы можем теперь считать все хиты классифицированными как Протеобактерии (taxid 1224) и выполнить тот же анализ распределения на уровне классов.
% === consider only Proteobacteria pb = taxidList(rank{1} == 1224); pbRank = findTaxoRank(pb, mmRankObj, mmParentObj, 8); pbList = unique(pbRank); pbT = accumarray(pbRank, 1); pbCount = pbT(pbList); pbNames = SN(pbList); % === filter out if less than 5 counts h = pbCount > 5; pbCount = pbCount(h); n2 = length(pbCount); pbNames = pbNames(h)
pbNames =
5x1 cell array
{'Gammaproteobacteria' }
{'Alphaproteobacteria' }
{'Betaproteobacteria' }
{'Deltaproteobacteria' }
{'Epsilonproteobacteria'}
Чтобы представлять различные филюмы и класс Протеобактерий в том же графике, мы должны создать матрицу смежности, таким образом, что все филюмы являются дочерними элементами корня, и все классы Протеобактерий являются дочерними элементами узла Протеобактерий. В графике метки включают имена классов и количество случаев в отчете BLAST.
% === find Proteobacteria node x = strcmp('Proteobacteria', phylaNames); % === combine names and counts numNodes = n1 + n2; allNames(1:n1) = phylaNames; allNames(n1+1:numNodes) = pbNames; allCount(1:n1) = count{1}(k); allCount(n1+1:numNodes) = pbCount; % === create labels for nodes (scientific name and count) labels = cell(1,numNodes); for node = 1:numNodes labels{node} = [allNames{node} ' (' num2str(allCount(node)) ')']; end % === create graph CM = zeros(numNodes); CM(1,2:n1) = 1; CM(x, n1+1:numNodes) = 1; CM(x,x) = 0; bg = biograph(CM, labels, 'showArrows', 'off'); bg.view % === clear memory mapped variables clear mmParentObj mmRankObj mt
[1] Родная мать, J.C., и др., "Экологическое упорядочивание ружья генома Саргассового моря", Наука, 304 (5667):66-74, 2004.