Введение

Дистальная кишка человека является наиболее плотной природной бактериальной экосистемой, известной на сегодняшний день. Его размер - до 100 триллионов клеток - намного превышает размер всех других микробных сообществ человеческого организма. Недавние исследования показали, что кишечная микробиота помогает регулировать энергетический баланс, как путем извлечения калорий из неперевариваемых компонентов, так и путем контроля накопления энергии в адипоцитах. Кроме того, кишечная микробиота участвует в множестве биологических процессов, начиная от синтеза незаменимых витаминов и заканчивая метаболизмом углеводов, липидов и других ксенобиотиков, которые мы принимаем.

В этом примере будут использоваться два набора данных. Первый набор данных состоит из данных, полученных в результате анализа дистального кишечного микробиома двух взрослых американских субъектов [1]. Он включает филогенетическое исследование микробных сообществ и функциональный анализ метаболических функций, представленных идентифицированным генофондом. Переменные, включенные в dataset1 описаны ниже. Обратите внимание, что таксономические назначения представлены в виде номинального категориального массива.

load gutmicrobiomedata.mat

% === first data set variables
rank1 = dataset1.rank1; % superkingdom assignments of each hit
rank2 = dataset1.rank2; % phylum assignments of each hit
rank3 = dataset1.rank3; % class assignments of each hit
subjF = dataset1.subjF; % number of hits in female subject
subjM = dataset1.subjM; % number of hits in male subject

Таксономическое профилирование дистального микробиома кишечника взрослого человека

Мы выполняем таксономическое профилирование первого набора данных, рассматривая таксономическое назначение контигов в соответствии с лучшим хитом BLASTX.

Мы начинаем с вычисления количества присвоенных попаданий, которые принадлежат каждому суперкингдому для каждого субъекта.

l1 = getlabels(rank1); % superkingdom labels
n1 = numel(l1);
count1 = zeros(n1,1); % number of hits for each superkingdom and subject

for i = 1:n1
    obs = rank1 == l1{i};
    count1(i,1) = sum(subjF(obs)); count1(i,2) = sum(subjM(obs));
end

% === plot
figure()
barh(count1);
colormap(summer);
ax = gca;
ax.YTickLabel = l1;
xlabel('Number of hits');
title('Superkingdom assignment of best-BLASTX-hits');
legend({'SubjF', 'SubjM'}, 'Location','southeast');

Как вы можете видеть из графика бара, микробное сообщество, живущее в дистальном кишечном микробиоме, преимущественно бактериальной природы. Различия, наблюдаемые между субъектами мужского и женского пола, не могут быть устранены из-за ограниченного размера выборки субъекта и возможности того, что эти различия могут быть связаны с генотипом или образом жизни хозяина.

Теперь мы повторяем анализ на уровне филума.

l2 = getlabels(rank2); % phylum labels
n2 = numel(l2);
count2 = zeros(n2,1); % number of hits for each phylum and subject

for i = 1:n2
    obs = rank2 == l2{i};
    count2(i,1) = sum(subjF(obs)); count2(i,2) = sum(subjM(obs));
end

% === plot
figure()
barh(count2);
colormap(summer);
ax = gca;
ax.YTick = 1:n2;
ax.YTickLabel = l2;
xlabel('Number of hits');
title('Phylum assignment of best-BLASTX-hits');
legend('SubjF', 'SubjM');

Бактериальные филотипы отнесены в основном к двум делениям, Firmicutes и Actinobacteria. Относительная недостаточность назначений Bacteroidetes вступает в противоречие с данными других исследований, но несоответствие может быть вызвано известными отклонениями используемых методов фекального лизиса и экстракции ДНК.

Наконец, выполните те же шаги для назначений на уровне класса.

l3 = getlabels(rank3); % class labels
n3 = numel(l3);
count3 = zeros(n3,1); % number of hits for each class and subject

for i = 1:n3
    obs = rank3 == l3{i};
    count3(i,1) = sum(subjF(obs));
    count3(i,2) = sum(subjM(obs));
end

% === plot
figure();
barh(count3);
colormap(summer);
ax = gca;
ax.YTick = 1:n3;
ax.YTickLabel = l3;
xlabel('Number of hits');
title('Class assignment of best-BLASTX-hits');
legend('SubjF', 'SubjM');

Таксономическое распределение на уровне классов выявляет обилие бактериальных филотипов в группах Clostridia и Bacilli, а также Actinobacteria и Methanobacteria.

Объединение таксономического распределения и базовой классификации

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

Чтобы построить такой граф, нам нужно определить матрицу связности CM представление отношений «родитель-потомок» между узлами. Мы идентифицируем филу (детей), принадлежащую каждому суперкингдому (родителю), и, в свою очередь, классы (детей), принадлежащие каждому филуму (родителю).

L = nominal([l1 l2 l3]);
N = n1 + n2 + n3;
CM = zeros(N, N); % connectivity matrix

% === populate CM with relationships between superkingdoms and phyla
for i = 1:n1
    obs = rank1 == l1{i}; % entries classified in a given superkingdom
    from = find(L == l1{i});  % parent node
    subobs =  unique(rank2(obs)); % phyla in a given superkigdom
    for j = 1:numel(subobs)
        to = find(L == subobs(j)); % child node
        CM(from, to) = 1;
    end
end

% === populate CM with relationships between phyla and classes
for i = 1:n2
    obs = rank2 == l2{i}; % entries classified in a given phylum
    from = find(L == l2{i});
    from = from(end);
    subobs =  unique(rank3(obs)); % classes in a given phylum
    for j = 1:numel(subobs)
        to = find(L == subobs(j));
        to = to(end);
        CM(from, to) = 1;
    end
end

% === create biograph object
bg = biograph(CM-diag(diag(CM)),[],'NodeAutoSize','off','ShowTextInNodes','Label');

Результирующий граф имеет 60 узлов и 58 рёбер. Каждый уровень в графе связан с заданным таксономическим рангом, а рёбра представляют лежащую в основе таксономическую классификацию. Теперь мы можем пометить каждый узел соответствующим таксономическим назначением и повернуть весь граф против часовой стрелки на 90 градусов.

% === label each node
set(bg.Nodes,'Size',[10 100]);
for i = 1:numel(bg.Nodes)
   bg.Nodes(i).Label = char(L(i));
end
dolayout(bg);

% === rotate counterclockwise by 90 degrees
for i = 1:numel(bg.Nodes)
    bg.Nodes(i).Position = fliplr(bg.Nodes(i).Position).*[-1 1];
    bg.Nodes(i).Size = [100 15];
end

% === redraw edges without changing node positions
bg.LayoutType = 'equilibrium';
dolayout(bg,'PathsOnly',true);
view(bg)

Чтобы включить данные распределения в график, для каждого назначения мы рассматриваем среднее количество попаданий между двумя предметами и соответствующий процент. Затем мы настраиваем цвет и размер каждого узла. В частности, листовые узлы представлены коробками, в то время как внутренние - кругами, а размер каждого узла пропорционален числу попаданий (в процентах), попадающих в данное таксономическое назначение.

% === compute distribution among ranks
count = [count1; count2; count3];
count = sum(count,2)/2; % avg between subjects
pct = (count + 1)/sum(count + 1) * 100; % add pseudocounts

% === determine color schema
t = accumarray(round(pct+1),1);
t(t>0) = 1:nnz(t);
colors = flipud(summer(nnz(t)));
cindex = t(round(pct+1));

% === customize color of nodes according to distribution
for i = 1:numel(bg.Nodes)
    mynode = bg.Nodes(i);
    if (numel(getdescendants(mynode))~= 1) % leaf
        mynode.Shape = 'circle';
    end
    mynode.Color = colors(cindex(i),:);
end

view(bg)

Из этого представления сразу видно, как большинство микробных сообществ состоят из бактерий, в частности Firmicutes, включая Clostridia и Bacilli.

Сравнительный функциональный анализ с использованием категорий KEGG

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

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

genome = dataset1.genome; % reference genomes considered
keggCat = dataset1.keggCat; % KEGG category assignment
keggData = dataset1.keggData; % odds ratio for each KEGG category relative to reference genomes

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

figure()
hi = imagesc(log(keggData));
colormap(redbluecmap);
colorbar;
ha = gca;
ha.XTick = 1:numel(genome);
ha.XTickLabel = genome;
ha.YTick = 1:numel(keggCat);
ha.YTickLabel = keggCat;

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

Сравнительный функциональный анализ с использованием типов COG

Категории COG, которые используют эволюционные отношения для группирования функционально связанных генов, могут быть использованы для выполнения функционального анализа вместо категорий KEGG, которые отображают ферменты на известные метаболические пути. Объект DataMatrix dm2 состоит из данных сравнительного метагеномного анализа дистального кишечного микробиома человека нескольких японских субъектов, включая младенцев, детей и взрослых [2]. Для справки сообщаются данные американских субъектов, рассмотренных выше, а также другие наборы метагеномных данных. Строки представляют различные наблюдения COG, в то время как столбцы представляют различные группы субъектов. Числовые данные состоят из нормализованных процентов совпадений в данной категории COG для данной группы субъектов.

get(dm2)

            Name: ''
        RowNames: {3868x1 cell}
        ColNames: {1x12 cell}
           NRows: 3868
           NCols: 12
           NDims: 2
    ElementClass: 'double'

Для каждой основной категории COG мы вычисляем совокупный нормализованный процент и сохраняем результаты в новом объекте DataMatrix с именем dm2Count.

codes = {'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', ...
    'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V'}; % COG code to consider

n = numel(codes);
N = size(dm2,2);
count = zeros(n,N);

for i = 1:n
    try
        count(i,:) = sum(dm2.(codes{i}));
    catch
        sprintf('COG code %s is not found in the data set.',codes{i});
	end
end

dm2Count = bioma.data.DataMatrix(count, codes, dm2.ColNames);

Чтобы выяснить, отличаются ли модели обогащения COG среди трех возрастных групп, мы сначала рассмотрим данные, связанные со взрослыми, детьми и младенцами.

group1 = {'Adult', 'Child', 'Infant'};
figure()
plot(dm2Count.(':')(group1), '.-', 'LineWidth', 2);
haxis = gca;
haxis.XTick = 1:n;
haxis.XTickLabel = codes;
legend(group1, 'location', 'northwest')
xlabel('COG categories');
ylabel('Normalized percentage of assigned genes');

Мы наблюдаем из этого сюжета, что взрослые субъекты и дети, по-видимому, имеют аналогичную картину обогащения с точки зрения категорий COG. Младенцы, с другой стороны, проявляют некоторые сингулярности для категорий G, K и L, соответствующие переносу и метаболизму углеводов, транскрипции и репликации соответственно.

В свете этого сродства между функциональными паттернами микробиома взрослого и ребенка мы рассматриваем комбинацию двух образцов (взрослый + ребенок) при выполнении сравнения с другими микробиомами образца окружающей среды.

group2 = {'Adult+Child', 'Soil', 'Whale Fall Ave.', 'Sargasso'};
figure()
plot(dm2Count.(':')(group2), '.-', 'LineWidth', 2);
haxis = gca;
haxis.XTick = 1:n;
haxis.XTickLabel = codes;
legend(group2, 'location', 'north');
xlabel('COG categories');
ylabel('Normalized percentage of assigned genes');

Наиболее яркие различия между паттерном обогащения микробиомов человека и паттерном других микробных сообществ окружающей среды связаны с COG категории G (углеводный метаболизм). Возможно, это связано с представлением, что микробиота толстой кишки использует неперевариваемые полисахариды и пептиды в качестве основного ресурса для производства энергии и биосинтеза клеточных компонентов. Примечательно также обогащение нескольких ферментов для репарации ДНК (категория COG L).

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

r = dm2Count.(':')(group2);
colors = {'b', 'g', 'r', 'k'};
theta = (linspace(0, 2*pi, n+1))';
figure();
hold on;

for i = 1:numel(group2)
    rho = [r(:,i); r(1,i)];
    plot(rho .* cos(theta), rho .* sin(theta), '-', 'Color', colors{i}, 'LineWidth', 2);
end

% === plot outside circle and labels
m = max(max(r));
for i = 1:n
    text( (m + .5) * cos(theta(i)), (m + .5) * sin(theta(i)), codes{i}, ...
        'HorizontalAlignment', 'center');
end

theta = (linspace(0,2*pi,100))';
plot(m * cos(theta), m * sin(theta), 'k-');

axis equal
axis([-1 1 -1 1] * (m+1))
axis off

legend(group2, 'location', 'NorthEastOutside')

Кластеризация микробиомов на основе их функциональных профилей

Мы можем изучить взаимосвязь между кишечными микробиомами человека и другими микробиомами окружающей среды, используя значения обогащения для каждого COG. Мы создаем иерархическое дерево кластера, используя полный алгоритм связывания и матрицу расстояний, созданную с учетом корреляции между точками данных. Рассматриваемые образцы включают: взрослый, детский, младенческий, американский, почвенный, китовый (1, 2 и 3) и саргассо.

group3 = {'Adult', 'Child', 'Infant', 'American (SubjF)', 'American (SubjM)', ...
    'Soil', 'Whale Fall 1', 'Whale Fall 2', 'Whale Fall 3', 'Sargasso'};

z = linkage((dm2Count.(':')(group3))', 'complete', 'correlation');
dendrogram(z, 'orientation', 'left', 'labels', group3, 'colorthreshold', 'default')

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

Ссылки

Gill, S., et al., «Метагеномный анализ дистального кишечного микробиома человека», Science, 312 (5778): 1355-9, 2006.

Курокава, К., и др., «Сравнительная метагеномика выявила обычно обогащенные наборы генов в микробиомах кишечника человека», исследование ДНК, 14 (4): 169-81, 2007.