Идентификация дифференцированно описанных генов из данных RNA-Seq

В этом примере показано, как протестировать данные RNA-Seq на дифференцированно описанные гены с помощью отрицательной биномиальной модели.

Введение

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

Набор данных для этого примера включает данные RNA-Seq, полученные в эксперименте, описанном Бруксом и др. [1]. Авторы исследовали эффект сокрушительного удара siRNA pasilla, ген, который, как известно, играл важную роль в регулировании соединения у Дрозофилы melanogaster. Набор данных состоит из 2 биологических, реплицирует управления, которое (необработанные) выборки и 2 биологических реплицируют (обработанных) выборок сокрушительного удара.

Осмотр чтения считает таблицы для геномных функций

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

Включенный файл pasilla_count_noMM.mat содержит две таблицы с матрицами количества на генном уровне и на уровне экзона для каждой из продуманных выборок. Можно получить подобные матрицы с помощью функционального featurecount.

load pasilla_count_noMM.mat
% preview the table of read counts for genes
geneCountTable(1:10,:)
ans =

  10x6 table

          ID           Reference    untreated3    untreated4    treated2    treated3
    _______________    _________    __________    __________    ________    ________

    {'FBgn0000003'}     {'3R'}            0             1            1           2  
    {'FBgn0000008'}     {'2R'}          142           117          138         132  
    {'FBgn0000014'}     {'3R'}           20            12           10          19  
    {'FBgn0000015'}     {'3R'}            2             4            0           1  
    {'FBgn0000017'}     {'3L'}         6591          5127         4809        6027  
    {'FBgn0000018'}     {'2L'}          469           530          492         574  
    {'FBgn0000024'}     {'3R'}            5             6           10           8  
    {'FBgn0000028'}     {'X' }            0             0            2           1  
    {'FBgn0000032'}     {'3R'}         1160          1143         1138        1415  
    {'FBgn0000036'}     {'3R'}            0             0            0           1  

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

% preview the table of read counts for exons
exonCountTable(1:10,:)
ans =

  10x6 table

                   ID                    Reference    untreated3    untreated4    treated2    treated3
    _________________________________    _________    __________    __________    ________    ________

    {'FBgn0000003_2648220_2648518'  }     {'3R'}            0             0            0           1  
    {'FBgn0000008_18024938_18025756'}     {'2R'}            0             1            0           2  
    {'FBgn0000008_18050410_18051199'}     {'2R'}           13             9           14          12  
    {'FBgn0000008_18052282_18052494'}     {'2R'}            4             2            5           0  
    {'FBgn0000008_18056749_18058222'}     {'2R'}           32            27           26          23  
    {'FBgn0000008_18058283_18059490'}     {'2R'}           14            18           29          22  
    {'FBgn0000008_18059587_18059757'}     {'2R'}            1             4            3           0  
    {'FBgn0000008_18059821_18059938'}     {'2R'}            0             0            2           0  
    {'FBgn0000015_12758093_12760298'}     {'3R'}            1             2            0           0  
    {'FBgn0000017_16615461_16618374'}     {'3L'}         1807          1572         1557        1702  

Можно аннотировать и сгруппировать выборки путем создания логического вектора можно следующим образом:

samples = geneCountTable(:,3:end).Properties.VariableNames;
untreated = strncmp(samples,'untreated',length('untreated'))
treated = strncmp(samples,'treated',length('treated'))
untreated =

  1x4 logical array

   1   1   0   0


treated =

  1x4 logical array

   0   0   1   1

Графический вывод присвоений функции

Включенный файл также содержит таблицу geneSummaryTable со сводными данными присвоенных и неприсвоенных записей SAM. Можно построить основное распределение результатов подсчета путем рассмотрения количества чтений, которые присвоены данным геномным функциям (экзоны или гены для этого примера), а также количество чтений, которые являются неприсвоенными (i.e. не перекрывая функции) или неоднозначный (i.e. наложение нескольких функций).

st = geneSummaryTable({'Assigned','Unassigned_ambiguous','Unassigned_noFeature'},:)
bar(table2array(st)','stacked');
legend(st.Properties.RowNames','Interpreter','none','Location','southeast');
xlabel('Sample')
ylabel('Number of reads')
st =

  3x4 table

                            untreated3    untreated4     treated2      treated3 
                            __________    __________    __________    __________

    Assigned                1.5457e+07    1.6302e+07    1.5146e+07    1.8856e+07
    Unassigned_ambiguous    1.5708e+05    1.6882e+05    1.6194e+05    1.9977e+05
    Unassigned_noFeature    7.5455e+05    5.8309e+05    5.8756e+05    6.8356e+05

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

geneSummaryTable{'TotalEntries',:} - sum(geneSummaryTable{2:end,:})
ans =

       89516       95885       98207      104629

Графический вывод покрытия чтения через данную хромосому

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

% consider chromosome arm 2L
chr2L = strcmp(exonCountTable.Reference,'2L');
exonCount = exonCountTable{:,3:end};

% retrieve exon start positions
exonStart = regexp(exonCountTable{chr2L,1},'_(\d+)_','tokens');
exonStart = [exonStart{:}];
exonStart = cellfun(@str2num, [exonStart{:}]');

% sort exon by start positions
[~,idx] = sort(exonStart);

% plot read coverage along the genomic coordinates
figure;
plot(exonStart(idx),exonCount(idx,treated),'.-r',...
exonStart(idx),exonCount(idx,untreated),'.-b');
xlabel('Genomic position');
ylabel('Read count (exon level)');
title('Read count on Chromosome arm 2L');

% plot read coverage for each group separately
figure;
subplot(2,1,1);
plot(exonStart(idx),exonCount(idx,untreated),'.-r');
ylabel('Read count (exon level)');
title('Chromosome arm 2L (treated samples)');
subplot(2,1,2);
plot(exonStart(idx),exonCount(idx,treated),'.-b');
ylabel('Read count (exon level)');
xlabel('Genomic position');
title('Chromosome arm 2L (untreated samples)');

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

% load gene start and stop position information
load pasilla_geneLength
geneLength(1:10,:)
ans =

  10x5 table

          ID              Name        Reference    Start    Stop 
    _______________    ___________    _________    _____    _____

    {'FBgn0037213'}    {'CG12581'}       3R          380    10200
    {'FBgn0000500'}    {'Dsk'    }       3R        15388    16170
    {'FBgn0053294'}    {'CR33294'}       3R        17136    21871
    {'FBgn0037215'}    {'CG12582'}       3R        23029    30295
    {'FBgn0037217'}    {'CG14636'}       3R        30207    41033
    {'FBgn0037218'}    {'aux'    }       3R        37505    53244
    {'FBgn0051516'}    {'CG31516'}       3R        44179    45852
    {'FBgn0261436'}    {'DhpD'   }       3R        53106    54971
    {'FBgn0037220'}    {'CG14641'}       3R        56475    58077
    {'FBgn0015331'}    {'abs'    }       3R        58765    60763

% consider chromosome 3 ('Reference' is a categorical variable)
chr3 = (geneLength.Reference == '3L') | (geneLength.Reference == '3R');
sum(chr3)

% consider the counts for genes in chromosome 3
counts = geneCountTable{:,3:end};
[~,j,k] = intersect(geneCountTable{:,'ID'},geneLength{chr3,'ID'});
gstart = geneLength{k,'Start'};
gcounts = counts(j,:);

% sort according to ascending start position
[~,idx] = sort(gstart);

% plot read coverage by genomic position
figure;
plot(gstart(idx), gcounts(idx,treated),'.-r',...
	gstart(idx), gcounts(idx,untreated),'.-b');
xlabel('Genomic position')
ylabel('Read count');
title('Read count on Chromosome 3');
ans =

        6360

Нормализация количеств чтения

Количество чтения в данных RNA-Seq, как находили, было линейно связано с распространенностью расшифровок стенограммы [2]. Однако счет чтения для данного гена зависит не только от уровня экспрессии гена, но также и на общем количестве упорядоченных чтений и длина генной расшифровки стенограммы. Поэтому для того, чтобы вывести уровень экспрессии гена от количества чтения, нам нужно с учетом глубины секвенирования и генной длины расшифровки стенограммы. Один общий метод, чтобы нормировать количество чтения должен использовать RPKM (Чтение На Сопоставленный Kilobase) значения, где количество чтения нормировано на общее количество чтений, к которым дают (в миллионах) и длина каждой расшифровки стенограммы (в kilobases). Этот метод нормализации, однако, является не всегда эффективным, поскольку немногие, очень высоко описанные гены могут доминировать над общим количеством маршрута и скосить анализ выражения.

Лучший метод нормализации состоит из вычисления эффективного размера библиотеки путем рассмотрения фактора размера для каждой выборки. Путем деления количеств каждой выборки на соответствующие факторы размера мы приносим все значения количества к общей шкале, делая их сопоставимыми. Интуитивно, если выборка A будет упорядочена времена N глубже, чем выборка B, то количества чтения недифференцированно описанных генов, как ожидают, будут в среднем N временами выше в выборке, чем в выборке B, даже если не будет никакого различия в выражении.

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

% estimate pseudo-reference with geometric mean row by row
pseudoRefSample = geomean(counts,2);
nz = pseudoRefSample > 0;
ratios = bsxfun(@rdivide,counts(nz,:),pseudoRefSample(nz));
sizeFactors = median(ratios,1)
sizeFactors =

    0.9374    0.9725    0.9388    1.1789

% transform to common scale
normCounts = bsxfun(@rdivide,counts,sizeFactors);
normCounts(1:10,:)
ans =

   1.0e+03 *

         0    0.0010    0.0011    0.0017
    0.1515    0.1203    0.1470    0.1120
    0.0213    0.0123    0.0107    0.0161
    0.0021    0.0041         0    0.0008
    7.0315    5.2721    5.1225    5.1124
    0.5003    0.5450    0.5241    0.4869
    0.0053    0.0062    0.0107    0.0068
         0         0    0.0021    0.0008
    1.2375    1.1753    1.2122    1.2003
         0         0         0    0.0008

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

figure;

subplot(2,1,1)
maboxplot(log2(counts),'title','Raw Read Count','orientation','horizontal')
ylabel('sample')
xlabel('log2(counts)')

subplot(2,1,2)
maboxplot(log2(normCounts),'title','Normalized Read Count','orientation','horizontal')
ylabel('sample')
xlabel('log2(counts)')

Вычисляя среднее значение, дисперсию и изменение сгиба

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

% consider the mean
meanTreated = mean(normCounts(:,treated),2);
meanUntreated = mean(normCounts(:,untreated),2);

% consider the dispersion
dispTreated = std(normCounts(:,treated),0,2) ./ meanTreated;
dispUntreated = std(normCounts(:,untreated),0,2) ./ meanUntreated;

% plot on a log-log scale
figure;
loglog(meanTreated,dispTreated,'r.');
hold on;
loglog(meanUntreated,dispUntreated,'b.');
xlabel('log2(Mean)');
ylabel('log2(Dispersion)');
legend('Treated','Untreated','Location','southwest');

Учитывая небольшое количество реплицирует, не удивительно ожидать, что дисперсионные значения рассеиваются с некоторым отклонением вокруг истинного значения. Часть этого отклонения отражает отклонение выборки и некоторые отражения истинная изменчивость среди экспрессий гена выборок.

Можно посмотреть на различие экспрессии гена среди двух условий, путем вычисления изменения сгиба (FC) для каждого гена, i.e. отношение между количествами в обработанной группе по количествам в необработанной группе. Обычно эти отношения рассматриваются в шкале log2, так, чтобы любое изменение было симметрично относительно нуля (e.g. отношение 1/2 или 2/1 соответствует-1 или +1 в логарифмической шкале).

% compute the mean and the log2FC
meanBase = (meanTreated + meanUntreated) / 2;
foldChange = meanTreated ./ meanUntreated;
log2FC = log2(foldChange);

% plot mean vs. fold change (MA plot)
mairplot(meanTreated, meanUntreated,'Type','MA','Plotonly',true);
set(get(gca,'Xlabel'),'String','mean of normalized counts')
set(get(gca,'Ylabel'),'String','log2(fold change)')
Warning: Zero values are ignored 

Возможно аннотировать значения в графике с соответствующими названиями генов, в интерактивном режиме избранными генами и списками генов экспорта к рабочей области путем вызова mairplot функционируйте, как проиллюстрировано ниже:

mairplot(meanTreated,meanUntreated,'Labels',geneCountTable.ID,'Type','MA');
Warning: Zero values are ignored 

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

% create table with statistics about each gene
geneTable = table(meanBase,meanTreated,meanUntreated,foldChange,log2FC);
geneTable.Properties.RowNames = geneCountTable.ID;
% summary
summary(geneTable)
Variables:

    meanBase: 11609x1 double

        Values:

            Min          0.21206
            Median        201.24
            Max       2.6789e+05

    meanTreated: 11609x1 double

        Values:

            Min                0
            Median        201.54
            Max       2.5676e+05

    meanUntreated: 11609x1 double

        Values:

            Min                0
            Median        199.44
            Max       2.7903e+05

    foldChange: 11609x1 double

        Values:

            Min             0 
            Median    0.99903 
            Max           Inf 

    log2FC: 11609x1 double

        Values:

            Min            -Inf
            Median    -0.001406
            Max             Inf

% preview
geneTable(1:10,:)
ans =

  10x5 table

                   meanBase    meanTreated    meanUntreated    foldChange      log2FC   
                   ________    ___________    _____________    __________    ___________

    FBgn0000003     0.9475        1.3808         0.51415         2.6857           1.4253
    FBgn0000008     132.69        129.48           135.9        0.95277        -0.069799
    FBgn0000014     15.111        13.384          16.838        0.79488         -0.33119
    FBgn0000015     1.7738       0.42413          3.1234        0.13579          -2.8806
    FBgn0000017     5634.6        5117.4          6151.8        0.83186         -0.26559
    FBgn0000018     514.08        505.48          522.67        0.96711        -0.048243
    FBgn0000024     7.2354        8.7189           5.752         1.5158          0.60009
    FBgn0000028    0.74465        1.4893               0            Inf              Inf
    FBgn0000032     1206.3        1206.2          1206.4        0.99983      -0.00025093
    FBgn0000036    0.21206       0.42413               0            Inf              Inf

% access information about a specific gene
myGene = 'FBgn0261570';
geneTable(myGene,:)
geneTable(myGene,{'meanBase','log2FC'})

% access information about a given gene list
myGeneSet = {'FBgn0261570','FBgn0261573','FBgn0261575','FBgn0261560'};
geneTable(myGeneSet,:)
ans =

  1x5 table

                   meanBase    meanTreated    meanUntreated    foldChange    log2FC 
                   ________    ___________    _____________    __________    _______

    FBgn0261570     4435.5       4939.1          3931.8          1.2562      0.32907


ans =

  1x2 table

                   meanBase    log2FC 
                   ________    _______

    FBgn0261570     4435.5     0.32907


ans =

  4x5 table

                   meanBase    meanTreated    meanUntreated    foldChange    log2FC 
                   ________    ___________    _____________    __________    _______

    FBgn0261570     4435.5       4939.1          3931.8          1.2562      0.32907
    FBgn0261573     2936.9       2954.8          2919.1          1.0122      0.01753
    FBgn0261575     4.3776       5.6318          3.1234          1.8031      0.85047
    FBgn0261560     2041.1       1494.3            2588         0.57738      -0.7924

Выведение дифференциального выражения с отрицательной биномиальной моделью

Определение, отличаются ли экспрессии гена в двух условиях статистически, состоит из отклонения нулевой гипотезы, что две выборки данных прибывают из распределений с равными средними значениями. Этот анализ принимает, что количества чтения моделируются согласно отрицательному биномиальному распределению (как предложено в [3]). Функциональный nbintest выполняет этот тип тестирования гипотезы с тремя возможными вариантами задать тип рычажного устройства между отклонением и средним значением.

Путем определения ссылки между отклонением и средним значением как идентичность, мы принимаем, что отклонение равно среднему значению, и количества моделируются распределением Пуассона [4].

tIdentity = nbintest(counts(:,treated),counts(:,untreated),'VarianceLink','Identity');
h = plotVarianceLink(tIdentity);

% set custom title
h(1).Title.String = 'Variance Link on Treated Samples';
h(2).Title.String = 'Variance Link on Untreated Samples';

В качестве альтернативы путем определения отклонения как суммы термина шума выстрела (i.e. среднее значение) и константа, умноженная на среднее значение в квадрате, количества моделируются согласно распределению, описанному в [5]. Постоянный термин оценивается с помощью всех строк в данных.

tConstant = nbintest(counts(:,treated),counts(:,untreated),'VarianceLink','Constant');
h = plotVarianceLink(tConstant);

% set custom title
h(1).Title.String = 'Variance Link on Treated Samples';
h(2).Title.String = 'Variance Link on Untreated Samples';

Наконец, путем рассмотрения отклонения как сумму термина шума выстрела (i.e. среднее значение), и локально регрессировал непараметрическая сглаженная функция среднего значения, количества моделируются согласно распределению, предложенному в [3].

tLocal = nbintest(counts(:,treated),counts(:,untreated),'VarianceLink','LocalRegression');
h = plotVarianceLink(tLocal);

% set custom title
h(1).Title.String = 'Variance Link on Treated Samples';
h(2).Title.String = 'Variance Link on Untreated Samples';

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

h = plotVarianceLink(tLocal,'compare',true);

% set custom title
h(1).Title.String = 'Variance Link on Treated Samples';
h(2).Title.String = 'Variance Link on Untreated Samples';

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

figure;
histogram(tLocal.pValue,100)
xlabel('P-value')
ylabel('Frequency')
title('P-value enrichment')

Отфильтруйте те гены с относительно низким количеством, чтобы наблюдать более универсальное распространение незначащих P-значений через область значений (0,1]. Обратите внимание на то, что это не влияет на распределение значительных P-значений.

lowCountThreshold = 10;
lowCountGenes = all(counts < lowCountThreshold, 2);
histogram(tLocal.pValue(~lowCountGenes),100)
xlabel('P-value')
ylabel('Frequency')
title('P-value enrichment without low count genes')

Несколько тестирований и настроенных P-значений

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

Корректировка Benjamini-Hochberg (BH) [6] является статистическим методом, который обеспечивает настроенное P-значение, отвечающее на следующий вопрос: какова была бы часть ложных положительных сторон, если бы все гены с настроенными P-значениями ниже заданного порога были рассмотрены значительными? Установите порог 0,1 для настроенных P-значений, эквивалентных, чтобы считать 10% ложными положительными сторонами как приемлемые, и идентифицировать гены, которые значительно описываются путем рассмотрения всех генов с настроенными P-значениями ниже этого порога.

% compute the adjusted P-values (BH correction)
padj = mafdr(tLocal.pValue,'BHFDR',true);

% add to the existing table
geneTable.pvalue = tLocal.pValue;
geneTable.padj = padj;

% create a table with significant genes
sig = geneTable.padj < 0.1;
geneTableSig = geneTable(sig,:);
geneTableSig = sortrows(geneTableSig,'padj');
numberSigGenes = size(geneTableSig,1)
numberSigGenes =

        1904

Идентификация наиболее отрегулированных и вниз отрегулированных генов

Можно теперь идентифицировать наиболее отрегулированные или вниз отрегулированные гены путем рассмотрения абсолютного изменения сгиба выше выбранного сокращения. Например, сокращение 1 в шкале log2 дает к списку генов, которые отрегулированы с 2 изменениями сгиба.

% find up-regulated genes
up = geneTableSig.log2FC > 1;
upGenes = sortrows(geneTableSig(up,:),'log2FC','descend');
numberSigGenesUp = sum(up)

% display the top 10 up-regulated genes
top10GenesUp = upGenes(1:10,:)

% find down-regulated genes
down = geneTableSig.log2FC < -1;
downGenes = sortrows(geneTableSig(down,:),'log2FC','ascend');
numberSigGenesDown = sum(down)

% find top 10 down-regulated genes
top10GenesDown = downGenes(1:10,:)
numberSigGenesUp =

   129


top10GenesUp =

  10x7 table

                   meanBase    meanTreated    meanUntreated    foldChange    log2FC      pvalue         padj   
                   ________    ___________    _____________    __________    ______    __________    __________

    FBgn0030173     3.3979       6.7957               0             Inf         Inf     0.0063115      0.047764
    FBgn0036822     3.1364       6.2729               0             Inf         Inf      0.012203      0.079274
    FBgn0052548      8.158       15.269          1.0476          14.575      3.8654    0.00016945     0.0022662
    FBgn0050495     6.8315       12.635          1.0283          12.287      3.6191     0.0018949      0.017972
    FBgn0063667     20.573       38.042          3.1042          12.255      3.6153    8.5037e-08    2.3845e-06
    FBgn0033764     91.969       167.61          16.324          10.268      3.3601    1.8345e-25    2.9174e-23
    FBgn0037290     85.845       155.46          16.228          9.5801        3.26    3.5583e-23    4.6941e-21
    FBgn0033733     7.4634       13.384          1.5424          8.6773      3.1172     0.0027276      0.024283
    FBgn0037191     7.1766       12.753          1.6003          7.9694      2.9945     0.0047803      0.038193
    FBgn0033943       6.95       12.319           1.581          7.7921       2.962     0.0053635      0.041986


numberSigGenesDown =

   181


top10GenesDown =

  10x7 table

                   meanBase    meanTreated    meanUntreated    foldChange    log2FC       pvalue         padj   
                   ________    ___________    _____________    __________    _______    __________    __________

    FBgn0053498     15.469             0         30.938                0        -Inf    9.8404e-11     4.345e-09
    FBgn0259236     6.8092             0         13.618                0        -Inf    1.5526e-05    0.00027393
    FBgn0052500     4.3703             0         8.7405                0        -Inf    0.00066783     0.0075343
    FBgn0039331     3.6954             0         7.3908                0        -Inf     0.0019558      0.018474
    FBgn0040697      3.419             0         6.8381                0        -Inf     0.0027378      0.024336
    FBgn0034972     2.9145             0         5.8291                0        -Inf     0.0068564       0.05073
    FBgn0040967     2.6382             0         5.2764                0        -Inf     0.0096039      0.065972
    FBgn0031923     2.3715             0         4.7429                0        -Inf      0.016164      0.098762
    FBgn0085359     62.473        2.9786         121.97         0.024421     -5.3557    5.5813e-33    1.5068e-30
    FBgn0004854     7.4674       0.53259         14.402          0.03698     -4.7571    8.1587e-05     0.0012034

Хорошая визуализация экспрессий гена и их значения дана путем графического вывода изменения сгиба по сравнению со средним значением в логарифмической шкале и окраски точек данных согласно настроенным P-значениям.

figure
scatter(log2(geneTable.meanBase),geneTable.log2FC,3,geneTable.padj,'o')
colormap(flipud(cool(256)))
colorbar;
ylabel('log2(Fold Change)')
xlabel('log2(Mean of normalized counts)')
title('Fold change by FDR')

Вы видите здесь это для слабо описанных генов (i.e. те с низкими средними значениями), ФРГ обычно высок, потому что низко считанные количества во власти Пуассоновского шума, и следовательно любая биологическая изменчивость утоплена в неопределенности от подсчета чтения.

Ссылки

[1] Ручьи и др. Сохранение RNA регулирующая карта между Дрозофилой и млекопитающими. Исследование генома 2011. 21:193-202.

[2] Mortazavi и др. Отображение и определение количества транскриптомов млекопитающих RNA-Seq. Методы природы 2008. 5:621-628.

[3] Андерс и др. Дифференциальный анализ выражения для последовательности считает данные. Биология генома 2010. 11:R106.

[4] Мариони и др. RNA-Seq: оценка технической воспроизводимости и сравнения с массивами экспрессии гена. Исследование генома 2008. 18:1509-1517.

[5] Робинсон и др. Модерируемый статистический тест для оценки различий в распространенности тега. Биоинформатика 2007. 23 (21):2881-2887.

[6] Benjamini и др. Управление ложным уровнем открытия: практический и мощный подход к нескольким тестирование. 1995. Журнал Королевского Статистического Общества, Серии B57 (1):289-300.

Смотрите также

| | |

Похожие темы

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