Введение

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

Чтение отчета хита BLASTX

В этом примере мы рассматриваем небольшое подмножество (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

Фильтрация хитов BLAST

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

% === 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 с таксономической информацией

Мы теперь интересуемся выполнением таксономической аннотации каждого хита в отчете 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

Классификация хитов BLAST научным названием

Каждый 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'            }

Сохранение аннотируемого отчета BLAST

Если мы определяем таксономическую классификацию для каждого хита, мы можем включать информацию в текстовый файл как показано ниже:

% === 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

Одна причина классифицировать хиты последовательности на отчет 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

Графический вывод таксономического распределения хитов BLAST

Если мы строим таксономическое распределение хитов на столбчатой диаграмме, мы замечаем, что у большинства таксонов есть небольшое число случаев.

% === 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

Классификация хитов BLAST более высоким таксономическим рангом

После того, как карты создаются для каждого хита, который сопоставлен с соответствием 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.