Введение

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

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

В этом примере мы рассматриваем небольшое подмножество (100 чтений) набора данных Sargasso Sea [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 и его связанным таксономическим идентификатором (таксидом). Поскольку в настоящее время насчитывается более 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 с таксономической информацией

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

Каждый таксид соответствует определённому таксону, которому присвоено научное имя и, возможно, различные синонимы. В наших классификационных целях нас интересуют только научные имена. Таким образом, мы извлекаем эту информацию и аннотируем каждый удар 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'            }

Сохранение аннотированного отчета 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 является изучение их таксономического распределения. Мы можем легко создать список организмов, которые представлены в отчете, их таксиды и их частота следующим образом:

% === 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 (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 не удивительно, потому что это доминирующая форма жизни в море Саргассо, Буркхолдерия и Шеванелла не ожидается, чтобы присутствовать в этой морской выборке, где питательные вещества и ресурсы низки, потому что они живут или в наземных условиях или в водных, богатых питательными веществами окружениях соответственно. Подробное обсуждение присутствия этих бактерий в море Саргассо смотрите в [1].

Размещение в ОЗУ информации о узле таксона

Часто, чтобы получить четкое видение таксономического распределения набора последовательностей, рассматриваются линнеевские категории выше, чем виды. Чтобы выполнить этот анализ, нам нужно создать карту между каждым таксидом и его присвоенным рангом, а также карту между каждым таксидом и таксидом его родительского узла, согласно схеме базы данных таксономии NCBI. Файлы, содержащие эту информацию, могут быть созданы с помощью функции helper 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 по более высокому таксономическому рангу

После создания карт для каждого попадания, которое связано с таксидом, соответствующим категории Linnean, более конкретным, чем целевой ранг, мы определяем родительского таксида и его ранг до достижения цели. Затем аннотируем попадание таксидом его более отдаленного предка. Синтетические конструкции или узлы без ранга считаются потомками корня. Эта процедура обхода таксономической иерархии выполняется вспомогательной функцией 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)

Теперь мы можем рассмотреть все хиты, классифицированные как Proteobacteria (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. В графике метки включают имена классов и количество вхождений в отчете 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., «Environmental genome shotgun секвенирования of the Sargasso sea», Science, 304 (5667): 66-74, 2004.