Этот пример показывает, как протестировать данные РНК-Seq для дифференциально экспрессируемых генов с использованием отрицательной биномиальной модели.
Типичный дифференциальный анализ экспрессии данных RNA-Seq состоит из нормализации необработанных данных и выполнения статистических тестов для отклонения или принятия нулевой гипотезы о том, что две группы образцов не показывают существенных различий в экспрессии генов. В этом примере показано, как проверить базовую статистику необработанных данных подсчета, как определить коэффициенты размера для нормализации подсчета и как вывести наиболее дифференциально экспрессируемые гены с помощью отрицательной биномиальной модели.
Набор данных для этого примера включает данные РНК-Seq, полученные в эксперименте, описанном Brooks et al. [1]. Авторы исследовали эффект нокаута siRNA pasilla, гена, который, как известно, играет важную роль в регуляции сплайсинга у Drosophila melanogaster. Набор данных состоит из 2 биологических повторов контрольных (необработанных) образцов и 2 биологических повторов выбитых (обработанных) образцов.
Отправной точкой для этого анализа данных RNA-Seq является матрица подсчета, где строки соответствуют интересующим геномным признакам, столбцы соответствуют заданным образцам, а значения представляют количество считываний, отображенных на каждый признак в данной выборке.
Включенный файл 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. Можно построить график базового распределения результатов подсчета, учитывая количество считываний, которые назначены заданным геномным признакам (экзонам или генам для этого примера), а также количество считываний, которые не назначены (то есть не перекрывают ни один признак) или неоднозначны (то есть перекрывают несколько признаков).
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 (Read Per Kilobase Mapped), где число считываний нормализуется на общее количество считываний (в миллионах) и длину каждого транскрипта (в килобазах). Этот метод нормализации, однако, не всегда эффективен, поскольку немногие, очень высоко экспрессируемые гены могут доминировать в общем количестве дорожек и искажать анализ экспрессии.
Лучший метод нормализации состоит в вычислении эффективного размера библиотеки путем учета коэффициента размера для каждой выборки. Деля отсчеты каждой выборки на соответствующие коэффициенты размера, мы приводим все отсчеты к общей шкале, делая их сопоставимыми. Интуитивно, если образец А секвенирован в N раз глубже, чем образец В, ожидается, что считанные числа неподифференциально экспрессированных генов будут в среднем в N раз выше в образце А, чем в образце В, даже если нет разницы в экспрессии.
Для оценки размерных факторов возьмем медиану соотношений наблюдаемых отсчетов к пересчетам псевдореферентного образца, подсчеты которого можно получить, рассматривая среднее геометрическое каждого гена по всем образцам [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) для каждого гена, то есть соотношение между количеством в обработанной группе и количеством в необработанной группе. Обычно эти отношения рассматриваются в шкале log2, так что любое изменение симметрично относительно нуля (например, отношение 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';
Альтернативно, задавая дисперсию как сумму члена шумов выстрела (то есть среднее) и константу, умноженную на среднее в квадрате, счетчики моделируются в соответствии с распределением, описанным в [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';
Наконец, рассматривая дисперсию как сумму члена шумов выстрела (то есть среднего) и локально регрессионную непараметрическую гладкую функцию среднего, счетчики моделируются в соответствии с распределением, предложенным в [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-значение указывает на вероятность того, что изменение экспрессии будет столь сильным, как наблюдаемое (или даже более сильным), при нулевой гипотезе, то есть условия не будут влиять на экспрессию генов. На гистограмме 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-значений так, чтобы вероятность наблюдения, по меньшей мере, одного значимого результата из-за случайности оставалась ниже желаемого уровня значимости.
Корректировка Бенджамини-Хохберга (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')
Здесь можно видеть, что для слабо экспрессированных генов (то есть генов с низкими средствами) FDR, как правило, высока, потому что низкие счетчики считывания преобладают под действием шума Пуассона, и, следовательно, любая биологическая изменчивость утопается в неопределенностях из подсчета считывания.
[1] Брукс и др. Сохранение РНК-регуляторной карты между дрозофилой и млекопитающими. Исследование генома 2011. 21:193-202.
[2] Мортазави и др. Картирование и количественное определение транскриптомов млекопитающих с помощью РНК-Seq. Методы природы 2008. 5:621-628.
[3] Anders et al. Анализ дифференциальных выражений для данных подсчета последовательностей. Биология генома 2010. 11:R106.
[4] Marioni et al. РНК-Seq: Оценка технической воспроизводимости и сравнение с массивами экспрессии генов. Исследование генома 2008. 18:1509-1517.
[5] Робинсон и др. Модерируемый статистический тест для оценки различий в численности метки. Биоинформатика 2007. 23(21):2881-2887.
[6] Benjamini et al. Контроль частоты ложных открытий: практичный и мощный подход к множественному тестированию. 1995. Журнал Королевского статистического общества, серия B57 (1): 289-300.
featurecount | mairplot | nbintest | plotVarianceLink