本發(fā)明涉及生物信息學(xué)與第二代測(cè)序技術(shù)領(lǐng)域，特別是關(guān)于不同實(shí)驗(yàn)室或平臺(tái)干擾下引入的批次效應(yīng)的去除方法，具體為一種去除測(cè)序數(shù)據(jù)噪聲的方法。

背景技術(shù)：

目前已有的去除測(cè)序數(shù)據(jù)噪聲算法有兩種，一種是去除不需要的變量方法。對(duì)于m個(gè)樣本和n組基因，基于對(duì)數(shù)線性模型，觀測(cè)的測(cè)序讀取計(jì)數(shù)在感興趣的已知協(xié)變量和不需要的變量的未知因素上回歸，利用數(shù)據(jù)的子集來估計(jì)不需要的變量并調(diào)整他們。另一種為替代變量分析方法。該算法結(jié)合奇異值分解和線性模型分析，通過線性模型刪除生物變量引入的差異后，對(duì)殘差矩陣通過奇異值分解估計(jì)特征值并確定重要的混雜變量后對(duì)其進(jìn)行移除。

以上兩種算法的缺陷具有以下三點(diǎn)缺陷：

(1)目前算法用于微陣列數(shù)據(jù)，不適用于第二代測(cè)序數(shù)據(jù)。

(2)沒有考慮模型中存在的異方差問題，數(shù)據(jù)處理的精度不準(zhǔn)確。

(3)目前算法對(duì)對(duì)計(jì)數(shù)值進(jìn)行對(duì)數(shù)變換，然而對(duì)數(shù)變化后形成的小計(jì)數(shù)值會(huì)存在內(nèi)在的噪聲，并且對(duì)數(shù)變換會(huì)加大較小計(jì)數(shù)值之間的差異，這些低計(jì)數(shù)值顯示樣本之間強(qiáng)大的相對(duì)差異。

技術(shù)實(shí)現(xiàn)要素：

根據(jù)現(xiàn)有技術(shù)存在的問題，本發(fā)明公開了一種去除測(cè)序數(shù)據(jù)噪聲的方法。

其采用如下技術(shù)方案：一種去除測(cè)序數(shù)據(jù)噪聲的方法，包括以下步驟：

S1：對(duì)原始數(shù)據(jù)集進(jìn)行過濾，應(yīng)用最小過濾原則，去除原始數(shù)據(jù)矩陣中數(shù)值為0或表達(dá)量極低的行；

S2：對(duì)過濾后的數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，去除系統(tǒng)偏差，將不同的結(jié)果根據(jù)全局?jǐn)?shù)值進(jìn)行調(diào)整，使個(gè)體之間的數(shù)據(jù)具有可比性；

S3：計(jì)算標(biāo)準(zhǔn)化數(shù)據(jù)后組變量和批次噪聲之間的相關(guān)性，進(jìn)行相關(guān)系數(shù)顯著性檢驗(yàn)，求出P值，若P<0.05，則有顯著相關(guān)性，否則沒有顯著相關(guān)性；

S4：若兩者相關(guān)，則對(duì)數(shù)據(jù)進(jìn)行r log變換；若兩：若兩者不相關(guān)，則先對(duì)數(shù)據(jù)進(jìn)行Z-Score處理，然后對(duì)處理后的數(shù)據(jù)進(jìn)行l(wèi)og變換，Z-Score的模型為：

其中：Y為處理前的數(shù)據(jù)，X為處理后的數(shù)據(jù)，μ為均值，ν為方差；

S5：確定替代變量。

進(jìn)一步的，所述S5步驟包括以下步驟：

S51：檢測(cè)批次效應(yīng)；

S52：計(jì)算替代變量。

進(jìn)一步的，所述S51步驟包括以下步驟：

①通過擬合模型x_ij＝μ_i+b_iy_j+e_ij來估計(jì)和并通過加權(quán)最小二乘法計(jì)算殘差形成m×n的殘差矩陣R；

②計(jì)算殘差矩陣的奇異值分解，即R＝UDV^T，其中U和V是A的特征向量，D表示A的特征值；

③使d_l為第l個(gè)特征值，其是D的第l個(gè)對(duì)角元素，l＝1,2,...,n，如果df是模型擬合的自由度，然后通過構(gòu)造最后的df特征值正好為零，將其刪除；對(duì)于特征基因k＝1,2,...,n-df設(shè)置觀察到的統(tǒng)計(jì)量為：

④通過置換R的每一行以形成矩陣R^*；

⑤擬合模型并計(jì)算殘差來形成m×n的空矩陣模型

⑥計(jì)算R⁰矩陣的奇異值分解

⑦對(duì)于R⁰中的特征基因k

⑧迭代4-7步驟共B次，得到空統(tǒng)計(jì)b＝1,2,...,B和k＝1,2,...,n-df；

⑨計(jì)算特征基因k的p值：

⑩對(duì)于用戶選擇的顯著性水平0≤α≤1，如果p_k≤α，則特征基因k為顯著性特征；否則，這些特征基因不顯著。

進(jìn)一步的，所述S52步驟包括以下步驟：

①通過擬合模型x_ij＝μ_i+b_iy_j+e_ij來估計(jì)和并通過加權(quán)最小二乘法計(jì)算殘差形成m×n的殘差矩陣R；

②計(jì)算殘差矩陣的奇異值分解R＝UDV^T，令e_k＝(e_k1,.....e_kn)^T是特征向量V的第k列，表示殘差特征基因，并且代表與主變量導(dǎo)致的信號(hào)無關(guān)的正交殘差信號(hào)；

設(shè)置為算法所確定的顯著特征基因數(shù)；

③在x_i(i＝1,2,...m)上回歸e_k并計(jì)算p值來檢測(cè)殘差特征基因和每個(gè)基因表達(dá)之間的關(guān)聯(lián)，p值測(cè)量殘差特征基因e_k和基因i的表達(dá)之間關(guān)聯(lián)的強(qiáng)度；

④令π₀是與e_k不相關(guān)的表達(dá)基因的比例，估計(jì)并估計(jì)與殘差特征基因相關(guān)的基因的數(shù)量為

⑤形成的簡(jiǎn)化矩陣

為與殘差特征基因k相關(guān)的基因數(shù)量的估計(jì)，計(jì)算X_r的特征基因，并用表示，j＝1,...,n；

⑥令即j^*是使得e_k和之間的相關(guān)性達(dá)到最大值所對(duì)應(yīng)的變量，并設(shè)置將替代變量的估計(jì)設(shè)置為與相應(yīng)的殘差特征基因最相關(guān)的簡(jiǎn)化矩陣的特征基因；

⑦在后續(xù)分析中，應(yīng)用模型

本發(fā)明具有以下有益效果：

(1)本發(fā)明的一種去除測(cè)序數(shù)據(jù)噪聲的方法，對(duì)原始第二代測(cè)序數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，降低了數(shù)據(jù)中的技術(shù)噪聲；

(2)本發(fā)明的一種去除測(cè)序數(shù)據(jù)噪聲的方法，當(dāng)組變量和批次相關(guān)時(shí)，使用正則對(duì)數(shù)變換后更穩(wěn)定，有助于多變量可視化和排序，效果更好；

(3)本發(fā)明的一種去除測(cè)序數(shù)據(jù)噪聲的方法，當(dāng)組變量和批次不相關(guān)時(shí)，用Z-Score處理后，提高了算法的精度；

(4)本發(fā)明的一種去除測(cè)序數(shù)據(jù)噪聲的方法，利用帶權(quán)重的最小二乘法求殘差矩陣，解決了模型中存在的異方差問題。

附圖說明

為了更清楚地說明本申請(qǐng)實(shí)施例或現(xiàn)有技術(shù)中的技術(shù)方案，下面將對(duì)實(shí)施例或現(xiàn)有技術(shù)描述中所需要使用的附圖作簡(jiǎn)單地介紹。顯而易見地，下面描述中的附圖僅僅是本申請(qǐng)中記載的一些實(shí)施例，對(duì)于本領(lǐng)域普通技術(shù)人員來講，在不付出創(chuàng)造性勞動(dòng)的前提下，還可以根據(jù)這些附圖獲得其他的附圖；

圖1為本發(fā)明所述去除測(cè)序數(shù)據(jù)噪聲的方法邏輯結(jié)構(gòu)示意圖；

圖2為本發(fā)明所述去除測(cè)序數(shù)據(jù)噪聲的方法步驟S5邏輯結(jié)構(gòu)示意圖；

圖3為組變量和批次效應(yīng)不相關(guān)數(shù)據(jù)差異表達(dá)結(jié)果的比較；

圖4為組變量和批次效應(yīng)相關(guān)數(shù)據(jù)差異表達(dá)結(jié)果的比較。

具體實(shí)施方式

為使本發(fā)明的技術(shù)方案和優(yōu)點(diǎn)更加清楚，下面結(jié)合本發(fā)明實(shí)施例中的附圖，對(duì)本發(fā)明實(shí)施例中的技術(shù)方案進(jìn)行清楚完整的描述。

實(shí)施例1

如圖1所示，一種去除測(cè)序數(shù)據(jù)噪聲的方法，包括以下步驟：

S1：對(duì)原始數(shù)據(jù)集進(jìn)行過濾，應(yīng)用最小過濾原則，去除原始數(shù)據(jù)矩陣中數(shù)值為0或表達(dá)量極低的行；

S2：對(duì)過濾后的數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，去除系統(tǒng)偏差，將不同的結(jié)果根據(jù)全局?jǐn)?shù)值進(jìn)行調(diào)整，使個(gè)體之間的數(shù)據(jù)具有可比性；

S3：計(jì)算標(biāo)準(zhǔn)化數(shù)據(jù)后組變量和批次噪聲之間的相關(guān)性，進(jìn)行相關(guān)系數(shù)顯著性檢驗(yàn)，求出P值，若P<0.05，則有顯著相關(guān)性，否則沒有顯著相關(guān)性；

S4：若兩者相關(guān)，則對(duì)數(shù)據(jù)進(jìn)行r log變換；若兩：若兩者不相關(guān)，則先對(duì)數(shù)據(jù)進(jìn)行Z-Score處理，然后對(duì)處理后的數(shù)據(jù)進(jìn)行l(wèi)og變換，Z-Score的模型為：

其中：Y為處理前的數(shù)據(jù)，X為處理后的數(shù)據(jù)，μ為均值，ν為方差；

S5：確定替代變量。

所述S1步驟中，由于計(jì)數(shù)數(shù)據(jù)矩陣中許多行僅包含0或者表達(dá)量極低，去除那些沒有含有或僅含有很少信息的行，能夠降低對(duì)象的大小，并增加處理的速度。

所述rlog變換為正則對(duì)數(shù)變換，所述log變換為對(duì)數(shù)變換，所述Z-Score處理為根據(jù)數(shù)據(jù)的均值和標(biāo)準(zhǔn)差進(jìn)行歸一化，經(jīng)過處理的數(shù)據(jù)符合正態(tài)分布。

如圖2所示，進(jìn)一步的，所述S5步驟包括以下步驟：

S51：檢測(cè)批次效應(yīng)；

S52：計(jì)算替代變量。

進(jìn)一步的，所述S51步驟包括以下步驟：

①通過擬合模型x_ij＝μ_i+b_iy_j+e_ij來估計(jì)和并通過加權(quán)最小二乘法計(jì)算殘差形成m×n的殘差矩陣R；

②計(jì)算殘差矩陣的奇異值分解，即R＝UDV^T，其中U和V是A的特征向量，D表示A的特征值；

③使d_l為第l個(gè)特征值，其是D的第l個(gè)對(duì)角元素，l＝1,2,...,n，如果df是模型擬合的自由度，然后通過構(gòu)造最后的df特征值正好為零，將其刪除；對(duì)于特征基因k＝1,2,...,n-df設(shè)置觀察到的統(tǒng)計(jì)量為：

④通過置換R的每一行以形成矩陣R^*；

⑤擬合模型并計(jì)算殘差來形成m×n的空矩陣模型

⑥計(jì)算R⁰矩陣的奇異值分解

⑦對(duì)于R⁰中的特征基因k

⑧迭代4-7步驟共B次，得到空統(tǒng)計(jì)b＝1,2,...,B和k＝1,2,...,n-df；

⑨計(jì)算特征基因k的p值：

⑩對(duì)于用戶選擇的顯著性水平0≤α≤1，如果p_k≤α，則特征基因k為顯著性特征；否則，這些特征基因不顯著。

進(jìn)一步的，所述S52步驟包括以下步驟：

①通過擬合模型x_ij＝μ_i+b_iy_j+e_ij來估計(jì)和并通過加權(quán)最小二乘法計(jì)算殘差形成m×n的殘差矩陣R；

②計(jì)算殘差矩陣的奇異值分解R＝UDV^T，令e_k＝(e_k1,.....e_kn)^T是特征向量V的第k列，表示殘差特征基因，并且代表與主變量導(dǎo)致的信號(hào)無關(guān)的正交殘差信號(hào)；

設(shè)置為算法所確定的顯著特征基因數(shù)；

③在x_i(i＝1,2,...m)上回歸e_k并計(jì)算p值來檢測(cè)殘差特征基因和每個(gè)基因表達(dá)之間的關(guān)聯(lián)，p值測(cè)量殘差特征基因e_k和基因i的表達(dá)之間關(guān)聯(lián)的強(qiáng)度；

④令π₀是與e_k不相關(guān)的表達(dá)基因的比例，估計(jì)并估計(jì)與殘差特征基因相關(guān)的基因的數(shù)量為

⑤形成的簡(jiǎn)化矩陣

為與殘差特征基因k相關(guān)的基因數(shù)量的估計(jì)，計(jì)算X_r的特征基因，并用表示，j＝1,...,n；

⑥令即j^*是使得e_k和之間的相關(guān)性達(dá)到最大值所對(duì)應(yīng)的變量，并設(shè)置將替代變量的估計(jì)設(shè)置為與相應(yīng)的殘差特征基因最相關(guān)的簡(jiǎn)化矩陣的特征基因；

⑦在后續(xù)分析中，應(yīng)用模型

實(shí)驗(yàn)結(jié)果：

以下實(shí)驗(yàn)為使用數(shù)據(jù)來比較算法的優(yōu)劣：

從ReCount網(wǎng)站下載了Pickrell和Montgomery研究的計(jì)數(shù)數(shù)據(jù)，并從Hapmap網(wǎng)站下載了譜系信息進(jìn)行分析。Montgomery研究為對(duì)具有北歐或西歐血統(tǒng)的猶他州居民進(jìn)行測(cè)序(HapMap中的CEU人群)，Pickrell為對(duì)尼日利亞伊巴丹的約魯巴人進(jìn)行測(cè)序(HapMap中的YRI人群)。通過把兩個(gè)不同群體的基因表達(dá)研究的數(shù)據(jù)組合起來，產(chǎn)生人為的批次效應(yīng)。在分析中把性別作為結(jié)果變量，然后使用不同的方法去除批次效應(yīng)。原始數(shù)據(jù)中，組變量和批次效應(yīng)幾乎完全正交，即兩者不相關(guān)，用不同的算法得到的結(jié)果如圖3所示。然后重新采樣數(shù)據(jù)，模擬了兩者相關(guān)的情況，得到的結(jié)果如圖4所示。

圖3和圖4中，算法1至算法5分別表示使用一種算法得到的結(jié)果，算法1為使用本發(fā)明的一種去除測(cè)序數(shù)據(jù)噪聲的方法得到的結(jié)果，算法2為使用RUVEmp去除批次效應(yīng)算法得到的結(jié)果，算法3為使用RUVRes去除批次效應(yīng)算法得到的結(jié)果，算法4為使用svaseq去除批次效應(yīng)算法得到的結(jié)果，算法5為使用Noadjustment去除批次效應(yīng)算法得到的結(jié)果。

從圖3中看到算法1的結(jié)果線條在其他算法的結(jié)果線條的上方，即本發(fā)明的算法去除批次效應(yīng)后的差異表達(dá)結(jié)果高于其他方法，因此說明該算法優(yōu)于其他算法。

圖4表示組變量和批次效應(yīng)相關(guān)時(shí)去除批次效應(yīng)后得到的差異表達(dá)結(jié)果，從圖4中看到算法1的結(jié)果線條在其他算法的結(jié)果線條的上方，即本發(fā)明的算法得到的差異表達(dá)分?jǐn)?shù)明顯高于其他方法，因此可以說明本算法明顯優(yōu)于其他算法。

綜合這兩幅圖可以得到，本發(fā)明的算法在組變量和批次效應(yīng)相關(guān)和不相關(guān)這兩種情況下，表現(xiàn)都優(yōu)于其他常用算法。

由于采用了上述技術(shù)方案，本發(fā)明提供的一種去除測(cè)序數(shù)據(jù)噪聲的方法，對(duì)原始第二代測(cè)序數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，降低了數(shù)據(jù)中的技術(shù)噪聲，當(dāng)組變量和批次相關(guān)時(shí)，使用正則對(duì)數(shù)變換后更穩(wěn)定，有助于多變量可視化和排序，效果更好，當(dāng)組變量和批次不相關(guān)時(shí)，用Z-Score處理后，提高了算法的精度，利用帶權(quán)重的最小二乘法求殘差矩陣，解決了模型中存在的異方差問題。

以上所述，僅為本發(fā)明較佳的具體實(shí)施方式，但本發(fā)明的保護(hù)范圍并不局限于此，任何熟悉本技術(shù)領(lǐng)域的技術(shù)人員在本發(fā)明揭露的技術(shù)范圍內(nèi)，根據(jù)本發(fā)明的技術(shù)方案及其發(fā)明構(gòu)思加以等同替換或改變，都應(yīng)涵蓋在本發(fā)明的保護(hù)范圍之內(nèi)。