RNA-seq miRNA-seq聯(lián)合分析
背景知識
肝癌細(xì)胞經(jīng)常會入侵門靜脈系統(tǒng),從而導(dǎo)致門靜脈癌栓,但是還沒有一個詳盡的研究來討論其中的作用機制,因此需要對肝癌組織(tumor),門靜脈組織(PVTT),癌旁組織(normal)進行采樣分析。
數(shù)據(jù)來源
數(shù)據(jù)來源于2017年5月24日清華大學(xué)更新的miRNA-seq,DNA methylation, CNV, RNA-seq
項目標(biāo)題:The molecularlandscape of hepatocellular carcinoma with portal vein tumor thrombosis
實驗設(shè)計:
提取了來自20個中國肝癌患者的腫瘤組織,門靜脈組織和癌旁組織,共計60個樣本,分別對其進行miRNA-seq,甲基化分析,拷貝數(shù)變異分析和RNA-seq分析。
數(shù)據(jù)下載網(wǎng)址:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE77276
RNA-seq數(shù)據(jù)分析
數(shù)據(jù)預(yù)處理
由于此數(shù)據(jù)原始數(shù)據(jù)sra太大,沒有表達矩陣,提供了測序序列reads經(jīng)過標(biāo)準(zhǔn)化以后,在每個基因上的數(shù)目(normalized_count),將各個樣本reads count文件合并就可以得到表達矩陣。
差異表達基因篩選
根據(jù)文獻所述,使用R包DESeq2篩選差異表達基因,DESeq2使用負(fù)二項分布產(chǎn)生的線性模型,具體原理可見如下網(wǎng)址
分組方式:源數(shù)據(jù)為肝癌組織(tumor),門靜脈組織(PVTT),癌旁組織(normal),然而由于門靜脈組織也屬于病變組織的一種,可以和tumor劃分為一類
最終在pvalue<0.001的條件下篩選出5676個差異表達基因,具體可以參見文件condition_treated_results.txt。
聚類熱圖
對前20個差異表達基因繪制聚類熱圖,可以發(fā)現(xiàn)normal和tumor明顯分開,這說明DESeq2找出來的差異表達基因還是蠻不錯的。
圖表 1聚類熱圖
深度分析
為了進一步探索數(shù)據(jù)和結(jié)果,繪制MA-plot,橫坐標(biāo)為每個基因上reads的數(shù)目(標(biāo)準(zhǔn)化后);縱坐標(biāo)為log2fold change,即變化的程度;每個點就是一個基因,紅色小點為pvalue<0.001的基因;只繪制了log2foldchange在(-3,+3)以內(nèi)的基因,即改變程度在(0.125,8)倍的基因,對于不在此范圍內(nèi)的基因,用三角形的標(biāo)志畫在邊界線上。
圖表 2MA-plot
可以從圖中看出來,黑色部分大致形成一個三角形,而紅色部分(差異表達基因)包裹在黑色三角形外圍。這說明用DESeq2的負(fù)二項分布模型找出來的差異表達基因,大部分都是reads數(shù)目多(測序深度高),且表達量差異很大的基因。
接下來繪制某一基因在不同組織的表達量。選取p值最小的那個基因
圖表 3某一基因的表達量
PVTT和tumor在差異表達基因篩選的時候合為兩組,此時繪圖的時候仍然將它們分開??梢钥吹紼NSG00000077152在normal和PVTT+tumor間表達量明顯不同。
再接下來可以進行主成分分析,對整個表達矩陣計算主成分,然后選取前面兩個主成分繪制PCA圖,可以看見PCA1代表了原本36%的信息量,PCA2代表了原本10%的信息量,然后normal和其他兩類比較能分得開,比起之前那次作業(yè)芯片數(shù)據(jù),這次的緊致性要好得多。
圖表 4PCA圖
miRNA-seq數(shù)據(jù)分析
數(shù)據(jù)處理過程和上面的RNA-seq一樣,把代碼切換一下目錄就成。
在2578個miRNA中,共有199個差異表達(pvalue<0.001),繪制MA-plot發(fā)現(xiàn)上調(diào)的居多,
圖表 5MA-plot
接下來,也對差異表達的部分做了聚類熱圖,發(fā)現(xiàn)對于差異表達的部分,兩組確實分得挺開的。
圖表 6聚類熱圖
接下來也挑了p值最小的miRNA繪制reads count圖,發(fā)現(xiàn)兩組之間的差異確實蠻明顯的。
圖表 7p值最小miRNA
最后,進行了主成分分析,繪制PCA圖,緊致性不如上面的RNA-seq,應(yīng)該是前兩個PCA代表的信息太少的緣故,第一主成分只有代表源數(shù)據(jù)19%的信息,第二主成分代表17%的信息,倆主成分加起來才有剛剛一個主成分那么多信息(RNA-seq第一主成分就有36%)。
圖表 8PCA圖
聯(lián)合分析
MAGIA(miRNA和基因整合分析)是一個進行靶預(yù)測、miRNA和基因表達數(shù)據(jù)整合分析的新的網(wǎng)絡(luò)工具。接下來,使用magia進行miRNA與基因相互作用的聯(lián)合分析。
網(wǎng)址:http://gencomp.bio.unipd.it/magia/analysis/
Step1
由于miRNA-seq和RNA-seq是來源相同的配對數(shù)據(jù),而且樣本數(shù)有60個。聯(lián)合分析算法選擇MATCHED:Mutual Information
MATCHED: Mutual Information: a classicinformation measure quantifying the mutual dependence of variables, includingnon-linear relationships. Suitable for large sample size (>20 needed).
Step2
接下來的預(yù)測方式選擇Pita和miRanda的交集
Pita score filter:-10 Miranda score filter:500(都是默認(rèn)值)
Step3
接下來將上面分析出來的差異表達矩陣分別上傳,分析即可。下面就是繪制出來的相互作用網(wǎng)絡(luò)圖。
圖表 9相互作用網(wǎng)絡(luò)
紅色三角形為miRNA,綠色圓形為基因。
紅色圈圈是看上去連線比較多的幾個miRNA,比較重要,名字分別是:hsa-miR-760、hsa-miR-1303 、hsa-miR-671-5p、hsa-miR-324-3p、hsa-miR-423-3p
還能做出來相互作用(interaction)的程度,下載為tsv文件,
就是一張包含了MicroRNA、Gene Symbol、MutualInformation的表,Mutual Information指互信息,是信息論里一種有用的信息度量,它可以看成是一個隨機變量中包含的關(guān)于另一個隨機變量的信息量,或者說是一個隨機變量由于已知另一個隨機變量而減少的不肯定性。
也就是說,在這里MutualInformation就可以看做兩者的相關(guān)程度。
就比如在下圖表的截圖中可以看出來,hsa-mir-1303和其對應(yīng)的靶基因DBF4B、hsa-mir-501-5p和其對應(yīng)的靶基因KIF2C就有很強的相關(guān)性。
圖表 10相互作用表
GO注釋
使用Gene Ontology官網(wǎng)上的在線注釋功能即可,輸入剛剛相互作用網(wǎng)絡(luò)interactions.tsv文件中的基因名,進行biologicalprocess(生化反應(yīng)),molecularfunction(分子功能),cellularcomponent(細(xì)胞定位)三方面的富集分析,通過富集分析可以找出在統(tǒng)計上顯著富集的GO Term,這些富集的條目有可能與研究的目前有關(guān)。
圖表 11biological process
圖表 12molecular function
圖表 13cellular component
看上去確實有一些相關(guān)的富集條目,比如分子功能:染色體綁定(chromatin binding);生化過程:有絲分裂過程(mitotic cell cycle);細(xì)胞定位:染色體部位(chromosomal part),這些都和癌癥細(xì)胞的產(chǎn)生有著重要關(guān)系。
結(jié)語
本次實驗使用的是配對的miRNA和mRNA表達譜文件,這給了我們一個通過生物信息學(xué)工具構(gòu)建miRNA-mRNA相互作用網(wǎng)絡(luò)的好機會,在系統(tǒng)層次的分析表明,我們找到了許多的重要miRNA和mRNA,這些對于肝癌起始和發(fā)展的過程中起著重要作用。這個全局的“miRNA-mRNA相互作用網(wǎng)絡(luò)”對于篩選miRNA靶基因和發(fā)現(xiàn)新的治療靶標(biāo)有著重要意義。
聯(lián)系客服