Этот пример показывает, как протестировать данные 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. Можно построить основное распределение результатов подсчета путем рассмотрения количества чтений, которые присвоены данным геномным функциям (экзоны или гены для этого примера), а также количество чтений, которые являются неприсвоенными (т.е. не перекрывающий функции) или неоднозначные (т.е. перекрывающий несколько функций).
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) для каждого гена, т.е. отношения между количествами в обработанной группе по количествам в необработанной группе. Обычно эти отношения рассматриваются в шкале 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
Определение, отличаются ли экспрессии гена в двух условиях статистически, состоит из отклонения нулевой гипотезы, что две выборки данных прибывают из дистрибутивов с равными средними значениями. Этот анализ принимает, что количества чтения моделируются согласно отрицательному bionomial распределению (как предложено в [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-значений так, чтобы вероятность наблюдения по крайней мере одного значительного результата случайно осталась ниже желаемого уровня значения.
Корректировка 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')
Вы видите здесь, что для слабо выраженных генов (т.е. те с низкими средними значениями), ФРГ обычно высок, потому что низко считанные количества во власти Пуассоновского шума, и следовательно любая биологическая изменчивость утоплена в неуверенности от подсчета чтения.
[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.