从RNA-Seq数据中识别差异表达基因

打开脚本

这个例子展示了如何使用负二项式模型测试差异表达基因的RNA-Seq数据。

简介

RNA-Seq数据的典型差异表达分析包括对原始计数进行归一化，并进行统计检验以拒绝或接受原假设，即两组样本在基因表达上没有显着差异。这个例子展示了如何检查原始计数数据的基本统计数据，如何确定计数归一化的大小因子，以及如何使用负二项式模型推断表达差异最大的基因。

本例的数据集包括Brooks等人描述的实验中获得的RNA-Seq数据。[1]。作者研究了siRNA敲除pasilla的影响，该基因已知在调控剪接中起重要作用黑腹果蝇．该数据集包括2个控制(未处理)样本的生物重复和2个敲除(处理)样本的生物重复。

检测基因组特征的读计数表

RNA-Seq数据分析的起点是一个计数矩阵，其中行对应感兴趣的基因组特征，列对应给定的样本，值表示映射到给定样本中每个特征的读取数。

包含的文件pasilla_count_noMM.mat包含两个表与计数矩阵在基因水平和外显子水平为每个考虑的样本。你可以用这个函数得到相似的矩阵featurecount．

负载pasilla_count_noMM.mat

预览基因的读计数表头(geneCountTable)

ans = 8x6表ID参考untreated3 untreated4 treated2 treated3 _____________ _________ __________ __________ ________ ________“FBgn0000003”“3R”01 12“FBgn0000008”“2R”142 117 138 132“FBgn0000014”“3R”20 12 10 19“FBgn0000015”“3R”24 01“FBgn0000017”“3L”6591 5127 4809 6027“FBgn0000018”“2L”469 530 492 574“FBgn0000024”“3R”5 6 10 8“FBgn0000028“X”002 1

请注意，当计数不进行总结时，报告单个特征(在本例中为外显子)及其元特征分配(在本例中为基因)，然后是开始和停止位置。

预览外显子的读计数表头(exonCountTable)

ans = 8 x6表ID引用untreated3 untreated4 treated2 treated3  _______________________________ _________ __________ __________ ________ ________ " FBgn0000003_2648220_2648518”“3 r”0 0 0 1”FBgn0000008_18024938_18025756”“2 r“0 1 0 2”FBgn0000008_18050410_18051199”“2 r“13 9 14 12”FBgn0000008_18052282_18052494”“2 r“4 2 5 0”FBgn0000008_18056749_18058222”“2 r“32”26日27日23 FBgn0000008_18058283_18059490”“2 r”14日18日29日22“FBgn0000008_18059587_18059757 2 r”1 4 3 0"FBgn0000008_18059821_18059938" "2R" 00 2 0

您可以通过创建一个逻辑向量来注释和分组示例，如下所示:

samples = geneCountTable(:，3:end).Properties.VariableNames;未经处理= strncmp(样本，“未经处理”长度(“未经处理”)处理= strncmp(样本，“治疗”长度(“治疗”）)

未经处理= 1x4逻辑阵列1 1 0 0处理= 1x4逻辑阵列0 0 1 1 1

绘制特征分配图

所包含的文件还包含一个表geneSummaryTable使用已分配和未分配SAM条目的摘要。您可以通过考虑分配给给定基因组特征(本例中为外显子或基因)的读取数，以及未分配(即不重叠任何特征)或不明确(即重叠多个特征)的读取数，绘制计数结果的基本分布。

st = geneSummaryTable({“分配”，“Unassigned_ambiguous”，“Unassigned_noFeature”}:)栏(table2array (st)”,“堆叠”）;传奇(st.Properties.RowNames ',“翻译”，“没有”，“位置”，“东南”）;包含(“样本”) ylabel (“读取数”）

st = 3x4表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，默认id由属性或元特征(默认为gene_id)和特征的开始和停止位置(默认为exon)组成。您可以使用外显子起始位置来绘制所考虑的任何染色体的读覆盖范围，例如染色体臂2L。

%考虑染色体臂2Lchr2L = strcmp(exonCountTable。参考,' 2 l '）;exonCount = exonCountTable{:，3:结束};检索外显子起始位置exonStart = regexp(exonCountTable{chr2L,1}，“_ (\ d +) _”，“令牌”）;exonStart = [exonStart{:}];exonStart = cellfun(@str2num， [exonStart{:}]');按起始位置排序外显子[~，idx] = sort(exonStart);沿基因组坐标绘制读取覆盖率图;情节(exonStart (idx) exonCount (idx、处理),“。r”，.．.exonStart (idx) exonCount (idx,未经处理的),“。”）;包含(“基因组的位置”）;ylabel (读取计数(外显子级别)）;标题(染色体臂2L读取计数）;%分别为每个组绘制读取覆盖率图;次要情节(2,1,1);情节(exonStart (idx) exonCount (idx,未经处理的),“。r”）;ylabel (读取计数(外显子级别)）;标题(染色体臂2L(处理过的样本)）;次要情节(2,1,2);情节(exonStart (idx) exonCount (idx、处理),“。”）;ylabel (读取计数(外显子级别)）;包含(“基因组的位置”）;标题(染色体臂2L(未经处理的样本)）;

或者，您可以考虑给定染色体中每个基因的起始位置来绘制读取覆盖率。该文件pasilla_geneLength.mat包含一个表，其中包含每个基因在相应基因注释文件中的起始和停止位置。

%负载基因启动和停止位置信息负载pasilla_geneLength头(geneLength)

ans = 8x5表ID名称参考开始停止_____________ _________ _________ __________ "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 " "DhpD" 3R 53106 54971

%考虑3号染色体(“参考”是一个分类变量)chr3 = (geneLength。参考= =' 3 l ') | (genlength .)参考= =“3 r”）;总和(chr3)%考虑3号染色体的基因计数counts = geneCountTable{:，3:结束};[~，j,k] = intersect(geneCountTable{:，“ID”}, geneLength {chr3,“ID”});gstart = geneLength{k，“开始”};Gcounts =计数(j，:);%按升序开始位置排序[~，idx] = sort(gstart);%绘图读取覆盖基因组位置图;情节(gstart (idx) gcounts (idx、处理),“。r”，.．.gstart (idx) gcounts (idx,未经处理的),“。”）;包含(“基因组的位置”) ylabel (“读计数”）;标题(“读取3号染色体计数”）;

Ans = 6360

正常化读计数

RNA-Seq数据中的读取计数已被发现与转录本[2]的丰度呈线性相关。然而，给定基因的读数不仅取决于基因的表达水平，还取决于测序的总读数和基因转录本的长度。因此，为了从读计数推断基因的表达水平，我们需要考虑测序深度和基因转录本长度。规范化读取计数的一种常用技术是使用RPKM(每千碱基映射的读取数)值，其中，读取计数由产生的读取总数(以百万为单位)和每个转录本的长度(以千碱基为单位)规范化。然而，这种归一化技术并不总是有效的，因为很少，非常高表达的基因可以支配总车道数和扭曲表达分析。

一个更好的归一化技术包括通过考虑每个样本的大小因子来计算有效库的大小。通过将每个样本的计数除以相应的大小因子，我们将所有的计数值带入一个公共尺度，使它们具有可比性。直观地，如果样品A的测序深度是样品B的N倍，即使在表达上没有差异，预计样品A的无差异表达基因的读取计数平均比样品B高N倍。

为了估计大小因素，取观察到的计数与伪参考样本的计数之比的中位数，其计数可以通过考虑所有样本中每个基因的几何平均值[3]来获得。然后，为了将观测到的计数转换为公共尺度，将每个样本中观测到的计数除以相应的大小因子。

%估计伪参考与几何平均值逐行pseudoRefSample =几何(计数，2);nz = pseudoRefSample > 0;ratio = bsxfun(@rdivide,counts(nz，:)，pseudoRefSample(nz));sizeFactors =中位数(比率，1)

sizeFactors = 0.9374 0.9725 0.9388 1.1789

%转换为公共比例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 00 0.0021 0.0008 1.2375 1.1753 1.2122 1.2003 000 0 0.0008

通过使用函数，您可以体会到这种归一化的效果箱线图表示统计措施，如中位数，四分位数，最小值和最大值。

图;次要情节(2,1,1)maboxplot (log2(计数)“标题”，“原始读计数”，“定位”，“水平”) ylabel (“样本”)包含(“log2(计数)”maboxplot(log2(normCounts)，“标题”，“规范化读计数”，“定位”，“水平”) ylabel (“样本”)包含(“log2(计数)”）

计算平均值，色散和折叠变化

为了更好地表征数据，我们考虑归一化计数的均值和离散度。读取计数的方差由两个项的和给出:跨样本的变化(原始方差)和通过计数读取测量表达式的不确定性(射击噪声或泊松)。原始方差项在高表达基因中占主导地位，而在低表达基因中占主导地位。您可以将经验分散值与规范化计数的平均值在对数刻度中绘制出来，如下所示。

%考虑平均数meanhandled = mean(normCounts(:，被处理)，2);meanunauthorized = mean(normCounts(:，未处理)，2);%考虑分散性distreating = std(normCounts(:，已处理)，0,2)./已处理;dispuntreating = std(normCounts(:，未处理)，0,2)./ mean未处理;在对数对数刻度上的%图图;重对数(meanTreated dispTreated,“r”。）;持有在；重对数(meanUntreated dispUntreated,“b”。）;包含(“log2(的意思)”）;ylabel (“log2(分散)”）;传奇(“治疗”，“未经处理”，“位置”，“西南”）;

考虑到重复的数量很少，预期离散值在真实值周围散布一些方差就不足为奇了。其中一些方差反映了抽样方差，一些反映了样本基因表达之间的真实可变性。

您可以通过计算每个基因的折叠变化(FC)来查看两种条件下基因表达的差异，即处理组的计数与未处理组的计数之间的比率。一般来说，这些比率在log2尺度下被考虑，因此任何变化都是关于零的对称(例如，1/2或2/1的比率对应于对数尺度下的-1或+1)。

%计算平均值和log2FCmeanBase = (mean processed + mean unauthorized) / 2;foldChange = meanprocessed ./ meanunauthorized;log2FC = log2(foldChange);百分比图平均与折叠变化(MA图)mairplot (meanTreated meanUntreated,“类型”，“马”，“Plotonly”,真正的);集(get (gca),“包含”)，“字符串”，“归一化计数平均值”甘氨胆酸)组(get (,“Ylabel”)，“字符串”，“log2(褶皱变化)”）

警告:忽略零值

方法可以用相应的基因名称注释图中的值，交互式地选择基因，并将基因列表导出到工作区mairplot函数如下所示:

mairplot (meanTreated meanUntreated,“标签”, geneCountTable。ID,“类型”，“马”）;

警告:忽略零值

将每个基因的均值和折叠变化信息方便地存储在一个表中。然后，您可以通过按基因名称索引表来访问有关给定基因或满足特定标准的一组基因的信息。

%创建关于每个基因的统计数据表geneTable = table(meanBase, meantreating, meantreating,foldChange,log2FC);geneTable.Properties.RowNames = genecountable . id;

%总结总结(geneTable)

变量:meanBase: 11609x1 double值:Min 0.21206中位数201.24 Max 2.6789e+05 mean治疗:11609x1 double值:Min 0中位数201.54 Max 2.5676e+05 mean治疗:11609x1 double值:Min 0中位数199.44 Max 2.7903e+05 foldChange: 11609x1 double值:Min 0中位数0.99903 Max Inf log2FC: 11609x1 double值:Min -Inf中位数-0.001406 Max Inf

%预览头(geneTable)

ans = 8x5表meanBase meanTreated mean未经处理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正正

%获取特定基因的信息myGene =“FBgn0261570”；geneTable (myGene:) geneTable (myGene, {“meanBase”，“log2FC”}）%访问给定基因列表的信息myGeneSet = [“FBgn0261570”，“FBgn0261573”，“FBgn0261575”，“FBgn0261560”];geneTable (myGeneSet:)

ans = 1x5 table meanBase mean processed mean未经处理的foldChange log2FC ________ ___________ _____________ __________ _______ FBgn0261570 4435.5 4939.1 3931.8 1.2562 0.32907 ans = 1x2 table meanBase log2FC ________ _______ FBgn0261570 4435.5 4939.1 0.32907 ans = 4x5 table meanBase mean processed mean未经处理的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.12341.8031 0.85047 FBgn0261560 2041.1 1494.3 2588 0.57738 -0.7924

用负二项式模型推导微分表达式

确定两种情况下的基因表达是否在统计上不同，包括拒绝原假设，即两个数据样本来自具有相等均值的分布。该分析假设读取计数是根据负二项分布建模的(如[3]中提出的那样)。这个函数rnaseqde使用三种可能的选项执行这种类型的假设检验，以指定方差和均值之间的联系类型。

通过指定方差和均值之间的联系作为一个恒等式，我们假设方差等于均值，并且计数由泊松分布[4]建模。“IDColumns”指定输入表中要附加到输出表中的列，以帮助保持数据的组织性。

diffTableIdentity = rnaseqde(geneCountTable，[“untreated3”，“untreated4”]、[“treated2”，“treated3”), VarianceLink =“身份”IDColumns =“ID”）;预览结果。头(diffTableIdentity)

ans = 8x6表ID Mean1 Mean2 Log2FoldChange PValue AdjustedPValue _____________ ______________ ______________ _________ ______________ "FBgn0000003" 0.51415 1.3808 1.4253 0.627 0.75892 "FBgn0000008" 135.9 129.48 -0.069799 0.48628 0.64516 "FBgn0000014" 16.838 13.384 -0.33119 0.44445 0.61806 "FBgn0000015" 3.1234 0.42413 -2.8806 0.05835 0.12584 "FBgn0000017" 6151.8 5117.4 -0.26559 2.864e-42 6.0233e-41 "FBgn0000018" 522.67 505.48 -0.048243 0.39015 0.5616 "FBgn0000024" 5.752 8.7189 0.60009 0.355110.52203 "FBgn0000028" 0 1.4893 Inf 0.252 0.39867

或者，通过将方差指定为shot noise项(即平均值)和常数乘以平均值的平方的总和，根据[5]中描述的分布对计数进行建模。常数项是使用数据中的所有行来估计的。

diffTableConstant = rnasedde (geneCountTable，[“untreated3”，“untreated4”]、[“treated2”，“treated3”), VarianceLink =“不变”IDColumns =“ID”）;预览结果。头(diffTableConstant)

ans = 8x6表ID Mean1 Mean2 Log2FoldChange PValue AdjustedPValue _____________ ______________ ______________ __________ ______________ "FBgn0000003" 0.51415 1.3808 1.4253 0.62769 0.7944 "FBgn0000008" 135.9 129.48 -0.069799 0.53367 0.72053 "FBgn0000014" 16.838 13.384 -0.33119 0.45592 0.68454 "FBgn0000015" 3.1234 0.42413 -2.8806 0.058924 0.16938 "FBgn0000017" 6151.8 5117.4 -0.26559 8.5529e-05 0.00077269 "FBgn0000018" 522.67 505.48 -0.048243 0.54834 0.73346 "FBgn0000024" 5.752 8.7189 0.600090.36131 0.58937 "FBgn0000028" 0 1.4893 Inf 0.2527 0.46047

最后，通过将方差视为镜头噪声项(即均值)和均值的局部回归非参数平滑函数的和，根据[3]中提出的分布对计数进行建模。

difftabllocalal = rnaseqde(geneCountTable，[“untreated3”，“untreated4”]、[“treated2”，“treated3”), VarianceLink =“本地”IDColumns =“ID”）;预览结果。头(diffTableLocal)

ans = 8x6表ID Mean1 Mean2 Log2FoldChange PValue AdjustedPValue _____________ ______________ ______________ _________ ______________ "FBgn0000003" 0.51415 1.3808 1.4253 11 "FBgn0000008" 135.9 129.48 -0.069799 0.67298 0.89231 "FBgn0000014" 16.838 13.384 -0.33119 0.6421 0.87234 "FBgn0000015" 3.1234 0.42413 -2.8806 0.22776 0.57215 "FBgn0000017" 6151.8 5117.4 -0.26559 0.0014429 0.014207 "FBgn0000018" 522.67 505.48 -0.048243 0.65307 0.88136 "FBgn0000024" 5.752 8.7189 0.60009 0.55154 0.81984"FBgn0000028" 0 1.4893 Inf 0.42929 0.7765

的输出rnaseqde包含一个p值向量。p值表示在原假设下发生与观察到的相同(甚至更强)的表达变化的概率，即条件对基因表达没有影响。在p值的直方图中，我们观察到低值的富集(由于差异表达基因)，而其他值则均匀分布(由于非差异表达基因)。值等于1的富集是由于基因计数非常低。

图;直方图(diffTableLocal.PValue 100)包含(“假定值”) ylabel (“频率”)标题(“假定值浓缩”）

过滤掉那些计数相对较低的基因，以观察非显著p值在范围内更均匀的分布(0,1)。注意，这并不影响显著p值的分布。

lowCountThreshold = 10;lowCountGenes = all(计数< lowCountThreshold, 2);直方图(diffTableLocal.PValue (~ lowCountGenes), 100)包含(“假定值”) ylabel (“频率”)标题(“没有低计数基因的p值富集”）

多次测试和调整p值

由于多重测试问题，用阈值p值来确定哪些折叠变化比其他变化更显著并不适合这种类型的数据分析。在执行大量同时进行的测试时，仅仅由于机会而获得显著结果的概率会随着测试的数量而增加。为了解释多重测试，对p值进行修正(或调整)，以使由于偶然性而观察到至少一个显著性结果的概率保持在所需的显著性水平以下。

Benjamini-Hochberg调整[6]是一种统计方法，它提供了一个调整后的p值，回答以下问题:如果所有调整后p值低于给定阈值的基因都被认为是显著的，那么假阳性的比例是多少?

的输出rnaseqde在“adjuststedpvalue”字段中包含一个调整p值的向量。默认情况下，p值使用Benjamini-Hochberg调整。中的“FDRMethod”名称-值参数rnaseqde可以设置为“storey”来执行storey的程序[7]。

将调整后的p值设置为0.1的阈值，相当于认为10%的假阳性是可以接受的，并通过考虑调整后p值低于此阈值的所有基因来识别显著表达的基因。

%创建一个含有重要基因的表sig = difftabllocalal。adjuststedpvalue < 0.1;difftablellocalsig = difftabllocalsig (sig，:);diffTableLocalSig = sortrows(diffTableLocalSig，“AdjustedPValue”）;numberSigGenes = size(diffTableLocalSig,1)

numberSigGenes = 1904

识别上调和下调最多的基因

现在，你可以通过考虑所选截止点以上的绝对折叠变化来确定调控最多或下调最多的基因。例如，在log2尺度中，1的截止值产生了2倍变化的上调基因列表。

发现上调基因up = diffTableLocalSig。Log2FoldChange > 1;upGenes = sortrows(diffTableLocalSig(up，:)，)，“Log2FoldChange”，“下”）;numberSigGenesUp = sum(up)%显示前10位上调基因top10GenesUp = upGenes(1:10，:)找到下调的基因down = diffTableLocalSig。Log2FoldChange < -1;downGenes = sortrows(diffTableLocalSig(down，:))，“Log2FoldChange”，“提升”）;numberSigGenesDown = sum(down)找到前10个下调基因top10GenesDown = downGenes(1:10，:)

numberSigGenesUp = 129 top10GenesUp = 10x6表ID Mean1 Mean2 Log2FoldChange PValue AdjustedPValue _____________ ____________ ______________ __________ ______________ "FBgn0030173" 0 6.7957 Inf 0.0063115 0.047764 "FBgn0036822" 0 6.2729 Inf 0.012203 0.079274 "FBgn0052548" 1.0476 15.269 3.8654 0.00016945 0.0022662 "FBgn0050495" 1.0283 12.635 3.6191 0.0018949 0.017972 "FBgn0063667" 3.1042 38.042 3.6153 8.5037e-08 2.3845e-06 "FBgn0033764" 16.324 167.61 3.3601 1.8345e-25 2.9174e-23 "FBgn0037290"16.228 155.46 3.26 3.5583e-23 4.6941e-21 "FBgn0033733" 1.5424 13.384 3.1172 0.0027276 0.024283 "FBgn0037191" 1.6003 12.753 2.9945 0.0047803 0.038193 "FBgn0033943" 1.581 12.319 2.962 0.0053635 0.041986 numberSigGenesDown = 181 top10GenesDown = 10x6表ID Mean1 Mean2 Log2FoldChange PValue AdjustedPValue _____________ _____________ ______________ __________ ______________ "FBgn0053498" 30.938 0 -Inf 9.8404e-11 4.345e-09 "FBgn0259236" 13.618 0 -Inf 1.5526e-05 0.00027393 "FBgn0052500" 8.7405 0-Inf 0.00066783 0.0075343“FBgn0039331”7.3908 0 -Inf 0.0019558 0.018474“FBgn0040697”6.8381 0 -Inf 0.0027378 0.024336“FBgn0034972”5.8291 0 -Inf 0.0068564 0.05073“FBgn0040967”5.2764 0 -Inf 0.0096039 0.065972“fbgn000031923”4.7429 0 -Inf 0.016164 0.098762“FBgn0085359”121.97 2.9786 -5.3557 5.5813e-33 1.5068e-30“FBgn0004854”14.402 0.53259 -4.7571 8.1587e-05 0.0012034

通过在对数尺度上绘制折叠变化与平均值的关系，并根据调整后的p值对数据点进行着色，可以很好地可视化基因表达及其意义。

图散射(log2 (geneTable.meanBase) diffTableLocal.Log2FoldChange 3, diffTableLocal.PValue,“o”) colormap(flipud(cool(256))) colorbar;ylabel (“log2(褶皱变化)”)包含(log2(归一化计数的平均值))标题(“罗斯福带来的巨大变化”）