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