RNA-seq miRNA-seq聯合分析

背景知識

肝癌細胞經常會入侵門靜脈系統，從而導致門靜脈癌栓，但是還沒有一個詳盡的研究來討論其中的作用機制，因此需要對肝癌組織(tumor)，門靜脈組織(PVTT)，癌旁組織(normal)進行取樣分析。

資料來源

資料來源於2017年5月24日清華大學更新的miRNA-seq,DNA methylation, CNV, RNA-seq

專案標題：The molecularlandscape of hepatocellular carcinoma with portal vein tumor thrombosis

實驗設計：

提取了來自20箇中國肝癌患者的腫瘤組織，門靜脈組織和癌旁組織，共計60個樣本，分別對其進行miRNA-seq，甲基化分析，拷貝數變異分析和RNA-seq分析。

RNA-seq資料分析

資料預處理

由於此資料原始資料sra太大，沒有表達矩陣，提供了測序序列reads經過標準化以後，在每個基因上的數目(normalized_count)，將各個樣本reads count檔案合併就可以得到表達矩陣。

差異表達基因篩選

根據文獻所述，使用R包DESeq2篩選差異表達基因，DESeq2使用負二項分佈產生的線性模型，具體原理可見如下網址

分組方式：源資料為肝癌組織(tumor)，門靜脈組織(PVTT)，癌旁組織(normal)，然而由於門靜脈組織也屬於病變組織的一種，可以和tumor劃分為一類

聚類熱圖

對前20個差異表達基因繪製聚類熱圖，可以發現normal和tumor明顯分開，這說明DESeq2找出來的差異表達基因還是蠻不錯的。

圖表 1聚類熱圖

深度分析

為了進一步探索資料和結果，繪製MA-plot，橫座標為每個基因上reads的數目(標準化後)；縱座標為log2fold change，即變化的程度；每個點就是一個基因，紅色小點為pvalue<0.001的基因；只繪製了log2foldchange在(-3,+3)以內的基因，即改變程度在(0.125,8)倍的基因，對於不在此範圍內的基因，用三角形的標誌畫在邊界線上。

圖表 2MA-plot

可以從圖中看出來，黑色部分大致形成一個三角形，而紅色部分（差異表達基因）包裹在黑色三角形外圍。這說明用DESeq2的負二項分佈模型找出來的差異表達基因，大部分都是reads數目多（測序深度高），且表達量差異很大的基因。

接下來繪製某一基因在不同組織的表達量。選取p值最小的那個基因

圖表 3某一基因的表達量

PVTT和tumor在差異表達基因篩選的時候合為兩組，此時繪圖的時候仍然將它們分開。可以看到ENSG00000077152在normal和PVTT+tumor間表達量明顯不同。

再接下來可以進行主成分分析，對整個表達矩陣計算主成分，然後選取前面兩個主成分繪製PCA圖，可以看見PCA1代表了原本36%的資訊量，PCA2代表了原本10%的資訊量，然後normal和其他兩類比較能分得開，比起之前那次作業晶片資料，這次的緊緻性要好得多。

圖表 4PCA圖

miRNA-seq資料分析

資料處理過程和上面的RNA-seq一樣，把程式碼切換一下目錄就成。

在2578個miRNA中，共有199個差異表達(pvalue<0.001)，繪製MA-plot發現上調的居多，

圖表 5MA-plot

接下來，也對差異表達的部分做了聚類熱圖，發現對於差異表達的部分，兩組確實分得挺開的。

圖表 6聚類熱圖

接下來也挑了p值最小的miRNA繪製reads count圖，發現兩組之間的差異確實蠻明顯的。

圖表 7p值最小miRNA

最後，進行了主成分分析，繪製PCA圖，緊緻性不如上面的RNA-seq，應該是前兩個PCA代表的資訊太少的緣故，第一主成分只有代表源資料19%的資訊，第二主成分代表17%的資訊，倆主成分加起來才有剛剛一個主成分那麼多資訊（RNA-seq第一主成分就有36%）。

圖表 8PCA圖

聯合分析

MAGIA（miRNA和基因整合分析）是一個進行靶預測、miRNA和基因表達資料整合分析的新的網路工具。接下來，使用magia進行miRNA與基因相互作用的聯合分析。

網址：http://gencomp.bio.unipd.it/magia/analysis/

Step1

由於miRNA-seq和RNA-seq是來源相同的配對資料，而且樣本數有60個。聯合分析演算法選擇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

接下來的預測方式選擇Pita和miRanda的交集

Pita score filter:-10 Miranda score filter:500（都是預設值）

Step3

接下來將上面分析出來的差異表達矩陣分別上傳，分析即可。下面就是繪製出來的相互作用網路圖。

圖表 9相互作用網路

紅色三角形為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指互資訊，是資訊理論裡一種有用的資訊度量，它可以看成是一個隨機變數中包含的關於另一個隨機變數的資訊量，或者說是一個隨機變數由於已知另一個隨機變數而減少的不肯定性。

也就是說，在這裡MutualInformation就可以看做兩者的相關程度。

就比如在下圖表的截圖中可以看出來，hsa-mir-1303和其對應的靶基因DBF4B、hsa-mir-501-5p和其對應的靶基因KIF2C就有很強的相關性。

圖表 10相互作用表

GO註釋

使用Gene Ontology官網上的線上註釋功能即可，輸入剛剛相互作用網路interactions.tsv檔案中的基因名，進行biologicalprocess(生化反應),molecularfunction(分子功能),cellularcomponent(細胞定位)三方面的富集分析，通過富集分析可以找出在統計上顯著富集的GO Term，這些富集的條目有可能與研究的目前有關。

圖表 11biological process

圖表 12molecular function

圖表 13cellular component

看上去確實有一些相關的富集條目，比如分子功能:染色體繫結(chromatin binding)；生化過程:有絲分裂過程(mitotic cell cycle)；細胞定位:染色體部位(chromosomal part)，這些都和癌症細胞的產生有著重要關係。

結語

本次實驗使用的是配對的miRNA和mRNA表達譜檔案，這給了我們一個通過生物資訊學工具構建miRNA-mRNA相互作用網路的好機會，在系統層次的分析表明，我們找到了許多的重要miRNA和mRNA，這些對於肝癌起始和發展的過程中起著重要作用。這個全域性的“miRNA-mRNA相互作用網路”對於篩選miRNA靶基因和發現新的治療靶標有著重要意義。