Введение

Человеческий дистальный кишечник является самой высокой плотностью, природной бактериальной экосистемой, известной на сегодняшний день. Его размер - до 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 и Metanobacteria.

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

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

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

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

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')

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

Ссылки

[1] Gill, S. et al., «Metagenomic Analysis of the Human Distal Gut Microbiome», Science, 312 (5778): 1355-9, 2006.

[2] Kurokawa, K., et al., «Сравнительная метагеномика, выявленная обычно обогащенные наборы генов в микробиомах кишечника человека», DNA Research, 14 (4): 169-81, 2007.