Этот пример иллюстрирует простой метагеномный анализ выборки данных из Саргассова моря. Требуется информация о таксономии, включенная в файлы. gi_taxid_prot.dmp, names.dmp и nodes.dmp (см. сжатый файл taxdump), которую можно загрузить с FTP-сайта таксономии NCBI.
Метагеномика - исследование таксономического состава образца организмов, полученного из общей среды обитания. Обычно он состоит из сравнения образцов последовательностей с базами данных известных последовательностей и использования информации таксономии для классификации видов образцов. Основные цели метагеномного анализа включают количественную оценку относительной численности известных видов и идентификацию неизвестных последовательностей, для которых еще не были идентифицированы родственники.
В этом примере рассматривается небольшое подмножество (100 считываний) набора данных Sargasso Sea [1], поиск которого проводился по базе данных NCBI-NR с использованием BLASTX с параметрами по умолчанию. Для удобства полученный отчет BLAST был сохранен и сжат в файл sargasso-sample1-100.rpt.gz, и он снабжен 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, мы читаем данные в блоках 1MB, сохраняем их как двоичный файл и затем используем функцию 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');
Мы можем получить доступ к таксиду первых десяти живых последовательностей 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 каждого удара и извлекаем связанный с ним таксид.
% === 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
Каждому таксиду соответствует определённый таксон, которому присвоено научное название и, возможно, различные синонимы. Для целей классификации нас интересуют только научные названия. Таким образом, мы извлекаем эту информацию и аннотируем каждый удар BLAST в отчете, используя научные названия, а не таксиды.
% === 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 является изучение их таксономического распределения. Мы можем легко создать список организмов, которые представлены в отчете, их таксиды и их частота следующим образом:
% === 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
Из простой статистики таксономического распределения мы наблюдаем, что наиболее представленным таксоном является Burkholderia sp.383 (таксид 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].
Часто для получения чёткого видения таксономического распределения набора последовательностей рассматриваются линнеевские категории выше видовых. Для выполнения этого анализа необходимо создать карту между каждым таксидом и присвоенным рангом, а также карту между каждым таксидом и таксидом родительского узла согласно схеме базы данных таксономии 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
После создания карт для каждого попадания, которое связано с таксидом, соответствующим категории Линнеева более специфической, чем целевой ранг, мы определяем родительский таксид и его ранг до достижения цели. Затем мы аннотируем хит таксидом его более далекого предка. Синтетические конструкции или узлы без ранга считаются потомками корня. Эту процедуру обхода таксономической иерархии выполняет вспомогательная функция findTaxoRank.
Предположим, мы заинтересованы в классификации наших хитов в соответствии с суперкингдом, к которому они принадлежат. После назначения таксиду сверхкингдома каждому попаданию мы группируем и подсчитываем вхождения следующим образом:
% === 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'}
Чтобы представить различные филы и класс Proteobacteria в одном графе, нам нужно создать матрицу связности, так что все филы являются потомками корня, а все классы Proteobacteria являются потомками узла Proteobacteria. На графике метки включают имена классов и количество вхождений в отчете 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] Venter, J.C., et al., «Экологическое секвенирование ружья генома Саргассова моря», Science, 304 (5667): 66-74, 2004.