廣州市黃埔區(qū)學(xué)大道攬月路廣州企業(yè)孵化器B座402
電話:020-85625352
手機(jī):18102256923、18102253682
Email:servers@gzscbio.com
Fax:020-85625352
QQ:386244141
項目名稱:micRNA測序與分析報告
所屬分類:生物信息學(xué)分析-報告解讀
聯(lián)系電話:020-85625352
QQ:386244141
Email:servers@gzscbio.com
技術(shù)服務(wù)描述
micRNA測序與分析報告
1. 建庫測序工作流程
1.1. 建庫測序流程
??MicroRNA (miRNA) 是一類長度約為20-24個核苷酸長度的具有調(diào)控功能的非編碼RNA,miRNA 是生物體內(nèi)一類重要的特殊分子,誘導(dǎo)基因沉默,參與細(xì)胞生長、發(fā)育、基因轉(zhuǎn)錄和翻譯等諸多生命活動的調(diào)控過程。
??從RNA樣品到最終分析結(jié)果數(shù)據(jù),需要經(jīng)過樣品檢測、建庫、測序和信息分析等過程,其中測序過程我們采用高通量測序illumina HiSeq/MiSeq測序平臺;illumina測序得到的原始圖像數(shù)據(jù)文件經(jīng)堿基識別(Base Calling)分析轉(zhuǎn)化為原始測序序列(Sequenced Reads),我們稱之為Raw Data或Raw Reads,結(jié)果以FASTQ(簡稱為fq)文件格式存儲,其中包含測序序列(reads)的序列信息以及其對應(yīng)的測序質(zhì)量信息。
??樣品質(zhì)量檢測:利用NanoDrop 2000分光光度計對RNA樣品進(jìn)行初步定量,Aglient 2100對樣品濃度精確定量。RNA樣品總量、濃度、完整性(RIN值,RNA Integrity Number)、純度(OD260/OD280比值)符合收樣立項標(biāo)準(zhǔn)后,進(jìn)入下一步的文庫構(gòu)建。
??文庫構(gòu)建:Total RNA PEG (polyethylene glycol)沉淀,連接3’接頭,PAGE (polyacrylamide gel electrophoresis) 膠分離回收,連接5’接頭,反轉(zhuǎn)錄,PCR擴(kuò)增,文庫切膠回收。
??文庫質(zhì)檢:Q-PCR方法對文庫進(jìn)行精確定量,文庫有效濃度>2nM即可上機(jī)測序。
??上機(jī)測序:將各文庫按照有效濃度及目標(biāo)下機(jī)數(shù)據(jù)量的需求pooling后,進(jìn)行SE75測序。
圖1: Small RNA 文庫構(gòu)建 Workflow
??Small RNA測序的接頭信息:
??RNA 5‘Adapter (RA5), part:
??5’-GTTCAGAGTTCTACAGTCCGACGATC-3’
??RNA 3’Adapter (RA3), part:
??5’-AGATCGGAAGAGCACACGTCT-3’ (1ug流程)
??5’-CTGTAGGCACCATCAATCAGATCGGAAGAGCACACGTCTG-3’ (10ug流程)
2. 信息分析流程
??Raw Data通過去接頭、去低質(zhì)量、去污染、統(tǒng)計序列長度分布等過程完成初級分析;將初級分析得到的clean reads與ncRNA庫比對,分類注釋ncRNA,并去除reads中的rNRA、tRNA、snRNA、snoRNA等;然后再將reads比對到參考基因組、miRBase上,計算miRNA表達(dá)量,注釋已知miRNA、預(yù)測新的miRNA;然后進(jìn)行差異miRNA、靶基因預(yù)測、功能注釋等后續(xù)分析。
miRNA生物信息分析流程圖
2.1. miRNA-seq 測序數(shù)據(jù)質(zhì)量評估
2.1.1. 數(shù)據(jù)質(zhì)量評估 及 數(shù)據(jù)過濾
??對于測序得到的rawdata,我們首先使用fastqc進(jìn)行數(shù)據(jù)質(zhì)量評估,通過結(jié)果,我們可以了解到測序數(shù)據(jù)reads質(zhì)量、長度、堿基分布,序列重復(fù)頻次等基本信息。
??測序得到的raw data,里面含有帶接頭的、低質(zhì)量的reads,為了保證后續(xù)信息分析的質(zhì)量,必須對raw reads進(jìn)行處理,得到clean reads。
??raw data的過濾步驟如下: 去除低質(zhì)量 reads; 去除有5’接頭污染的reads; 去除沒有3’接頭序列; 去除沒有插入片段的reads; 去除長度小于16nt的reads;
??數(shù)據(jù)經(jīng)過以上過濾后,我們對數(shù)據(jù)做一些基本統(tǒng)計。
結(jié)果:
結(jié)果說明 | 結(jié)果文件 |
---|---|
原始數(shù)據(jù) RawData QC 結(jié)果 | 00.QC/qc_rawdata/ |
過濾數(shù)據(jù) CleanData QC 結(jié)果 | 00.QC/qc_cleandata |
Trimmed Reads (unique) 結(jié)果文件 | 00.QC/*.trim.collapse.fa |
2.2. 數(shù)據(jù)長度分布
??一般來說,小RNA的長度區(qū)間為18~30nt,數(shù)據(jù)過濾完成后,我們統(tǒng)計了clean reads數(shù)目及各自占總reads的百分比,并繪制reads長度分布圖。reads長度分布圖能幫助我們判斷小RNA的種類,如miRNA集中在21或22nt,siRNA集中在24nt,piRNA集中在30nt。
Clean Data Reads 長度分布統(tǒng)計直方圖:
Figure 2.2: Small RNA Clean Reads長度分布圖
Clean Reads長度分布直方圖,橫坐標(biāo)為Reads的長度,縱坐標(biāo)為樣品中對應(yīng)該長度的total Reads數(shù)量。
2.3. miRNA基因組比對分析
??我們將Clean Reads 與 mature miRNA, miRNA hairpin, mRNA, mature & primary tRNA, snoRNA, rRNA, other non-coding RNA, and (optional) known RNA 使用bowtie進(jìn)行比對搜索,然后統(tǒng)計reads的比對情況,分析樣本中各類ncRNA的數(shù)目及比例情況,統(tǒng)計結(jié)果如下。過濾出待分析的miRNA以便后續(xù)分析,將clean reads中ncRNA盡可能地去除,便于后續(xù) miRNA 檢測及預(yù)測。
Reads類型分布圖
圖顯示的是各類型的RNA reads和miRNA reads占總distinct reads的比例。
與參考基因組比對質(zhì)控結(jié)果:
結(jié)果說明 | 結(jié)果文件 |
---|---|
Trimmed Reads (unique) 與參考基因組的比對結(jié)果(bam) | 01.Mapping/*.bam |
Trimmed Reads (unique) 與參考基因組的比對結(jié)果(bw) | 01.Mapping/*.bw |
各類型RNA統(tǒng)計情況 | 00.QC/annotation.report.csv |
與各個smallRNA比對的詳細(xì)結(jié)果 | 01.Mapping/mapped.csv |
2.4. 已知 miRNA 的檢測以及 novel miRNA 預(yù)測分析
??首先用去除了ncRNA后的Reads比對miRbase數(shù)據(jù)庫,比對上miRbase中 reads 用于檢測已知成熟的miRNA,并統(tǒng)計其表達(dá)量和RPM值。我們采用miRge軟件來進(jìn)行 miRNA reads 檢測注釋、novel miRNA的預(yù)測、miRNA表達(dá)量的統(tǒng)計。
miRNA分析結(jié)果匯總:
結(jié)果說明 | 結(jié)果文件 |
---|---|
基本注釋結(jié)果: | |
miRNA reads 變體類型注釋結(jié)果 | 02.miRNA/1.isomiR/sample_miRge3.gff |
miRNA reads 變體類型可視化統(tǒng)計 | 02.miRNA/1.isomiR/*svg |
miRNA堿基編輯分析結(jié)果: | |
A2I在各個樣本中統(tǒng)計匯總表 | 02.miRNA/2.a2IEditing/a2IEditing.report.csv |
A2I在各個樣本中統(tǒng)計匯總表的篩選及簡化信息 ( 即匯總表的 *.RPM.mismatch 、*.AtoI.percentage 列 ) | 02.miRNA/2.a2IEditing/a2IEditing.report.newform.csv |
A2I分析的序列及統(tǒng)計詳情信息 | 02.miRNA/2.a2IEditing/a2IEditing.detail.txt |
二級結(jié)構(gòu)預(yù)測結(jié)果: | |
novel miRNA 二級結(jié)構(gòu)預(yù)測結(jié)果 | 02.miRNA/3.novel_predict/*.R1_novel_miRNAs/ |
novel miRNA 二級結(jié)構(gòu)預(yù)測詳細(xì)信息表格匯總 | 02.miRNA/3.novel_predict/*.R1_novel_miRNAs/*_novel_miRNAs_report.csv |
以上miRNA各圖的可視化網(wǎng)頁:02.miRNA/miR_visualization.html
2.4.1. miRNA異構(gòu)體分析
??隨著高通量測序技術(shù)的發(fā)展,早期認(rèn)為的miRNA loci僅產(chǎn)生1條miRNA成熟體序列這一觀點被顛覆。對于某一miRNA來說,它并不是單一的序列,而是由一系列長度/序列及表達(dá)不同的異構(gòu)體 (isomiRs) 組成。這些isomiRs表達(dá)多樣且序列多樣,甚至引入多樣的5'端及種子區(qū)域。特定miRNA位點在疾病組織中可具有異常的表達(dá)模式,現(xiàn)已證實,部分isomiRs具有重要的生物學(xué)功能。
??IsomiRs目前主要分為三類:5′isomiRs
, 3′isomiRs
, polymorphic isomiRs
。IsomiRs的產(chǎn)生機(jī)制主要有:miRNA加工和成熟過程中Darsha和Dicer酶的不精確或選擇性剪切;3'端核苷酸添加;RNA編輯和單核苷酸多態(tài)性( single nucleotide polymorphism,SNP)。主要表現(xiàn)為:5'端修剪;3' 端修剪;3'端核苷酸附加和堿基替換。其中,5'端修剪和堿基替換可發(fā)生在種子區(qū)域內(nèi)部,產(chǎn)生“種子轉(zhuǎn)移”(seed shifting)。
??我們通過序列與miRbase比對的情況,將各個miRNA Reads進(jìn)行isomiR識別,并同時匯總至miRNA reads 注釋結(jié)果中,分類統(tǒng)計,并選取表達(dá)量較高的前20個mirna進(jìn)行各類IsomiRs表達(dá)豐度統(tǒng)計。
圖2.4.1 所有樣本的變體類型分布(isomiRs)
圖2.4.2 各個樣本的變體分布TOP20的豐度統(tǒng)計
2.4.2. miRNA 堿基編輯
??成熟miRNA序列的第2-8個堿基被稱作“種子”序列,保守性很高。若在這一區(qū)域發(fā)生堿基突變,則可能改變miRNA的靶基因作用位點。我們通過將miRNA序列與miRBase中已知miRNA前體以及成熟miRNA序列進(jìn)行比對(只允許一個位點的錯配)找出發(fā)生了堿基突變的miRNA,統(tǒng)計匯總在該miRNA的堿基突變位點,以及發(fā)生堿基編輯的Reads數(shù)量及百分比,并對統(tǒng)計概況結(jié)果進(jìn)行可視化。
圖2.4.3 鑒定為已存在mirna堿基編輯的在各個樣本中的相應(yīng)readCount占比統(tǒng)計圖
miRNA堿基編輯分析結(jié)果詳細(xì)說明:
A2I在各個樣本中統(tǒng)計匯總表
表頭 說明 miRNA
miRNA的名稱 A-to-I position in the miRNA
A2I檢測位置,即序列的第幾個堿基 miRNA sequence
miRNA的序列 Sample.readCount
該樣本miRNA的所有序列的readCount總和 Sample.readCount.canonical
該樣本miRNA的所有序列的readCount總和(篩選低質(zhì)量后) Sample.RPM.canonical
該樣本miRNA的所有序列的readCount總和(篩選低質(zhì)量后)的RPM值 Sample.readCount.mismatch
A2I檢測的序列readCount總和 Sample.RPM.mismatch
A2I檢測的序列readCount總和的RPM值 Sample.AtoI.percentage
A2I檢測的序列readCount占比,計算方法為: Sample.readCount.mismatch/Sample.readCount.canonical*100%
Sample.AtoI.adjusted.pValue
A2I檢測的置信度計算 表頭說明:
A2I在各個樣本中統(tǒng)計匯總表的篩選及簡化信息
( 即匯總表的*.RPM.mismatch
、*.AtoI.percentage
列,該表對應(yīng)上面的可視化結(jié)果 )表頭 說明 miRNA:position
miRNA和position位置信息 sample
樣本名 A-to-I percentage
A2I檢測的序列readCount占比,對應(yīng)第一個表中的 Sample.AtoI.percentage
log2RPM
A2I檢測的序列readCount總和的log2RPM值,對應(yīng)第一個表中的A2I檢測的序列
readCount總和的RPM值Sample.RPM.mismatch
,進(jìn)行l(wèi)og2計算結(jié)果文件:
02.miRNA/2.a2IEditing/a2IEditing.report.newform.csv
表頭說明:
A2I分析的序列及統(tǒng)計詳情信息
每一個miRNA塊,包含3小塊內(nèi)容:原始序列、篩選后的序列統(tǒng)計、篩選后的序列。重復(fù)的miRNA塊代表不同樣本,順序為 A2I在各個樣本中統(tǒng)計匯總表 的樣本信息。
每一行序列包含3個信息,
miRNA sequence
(reads序列),*.readCount
(Count計數(shù)),Filter
(是否篩選到后續(xù)進(jìn)行堿基編輯分析)結(jié)果說明:
2.4.3. novel miRNA 預(yù)測及二級結(jié)構(gòu)分析
??miRge根據(jù)序列特征,堿基配對原則,最小自由能等特征模型,進(jìn)行新穎的miRNA及二級結(jié)構(gòu)預(yù)測。預(yù)測示意圖如下。結(jié)果見:02.miRNA/3.novel_predict/
圖2.4.4 新miRNA的二級結(jié)構(gòu)預(yù)測示意圖
2.5. miRNA 定量分析
2.5.1. 樣本間相關(guān)性分析
??生物學(xué)重復(fù)通常是任何生物學(xué)實驗所必須的,目前主流期刊也基本要求生物學(xué)重復(fù)。生物學(xué)重復(fù)主要有兩個用途:一個是證明所涉及的生物學(xué)實驗操作不是偶然,而是可重復(fù)的。另一個是為了確保后續(xù)的差異分析得到更可靠的結(jié)果。樣品間基因表達(dá)水平相關(guān)性是檢驗實驗可靠性和樣本選擇是否合理的重要指標(biāo)。相關(guān)系數(shù)越接近1,表明樣品之間表達(dá)模式的相似度越高。Encode計劃建議皮爾遜相關(guān)系數(shù)的平方(R2)大于0.92(理想的取樣和實驗條件下)。具體的項目操作中,我們要求生物學(xué)重復(fù)樣品間R2至少要大于0.8,否則需要對樣品做出合適的解釋,或者重新進(jìn)行實驗。根據(jù)各樣本所有基因的表達(dá)值計算組內(nèi)及組間樣本的相關(guān)性系數(shù),繪制成熱圖,可直觀顯示組間樣本差異及組內(nèi)樣本重復(fù)情況。樣本間相關(guān)性系數(shù)越高,其表達(dá)模式越為接近,樣本相關(guān)性熱圖如下圖所示。
圖 2.5.1. 樣本間相關(guān)性熱圖
圖中橫縱坐標(biāo)為各樣本相關(guān)系數(shù)的平方
結(jié)果文件:
樣本間相關(guān)性熱圖結(jié)果:
04.DE/Quant/cor_pheatmap*
2.5.2. 主成分分析
??主成分分析(PCA)也常用來評估組間差異及組內(nèi)樣本重復(fù)情況,PCA采用線性代數(shù)的計算方法,對數(shù)以萬計的基因變量進(jìn)行降維及主成分提取。我們對所有樣本的基因表達(dá)值進(jìn)行PCA分析,如下圖所示。理想條件下,PCA圖中,組間樣本應(yīng)該分散,組內(nèi)樣本應(yīng)該聚在一起。
圖 2.5.2. 主成分分析結(jié)果圖
圖中橫坐標(biāo)為第一主成分,縱坐標(biāo)為第二主成分
結(jié)果文件:
主成分分析結(jié)果:
04.DE/Quant/pca*
2.6. miRNA 差異分析
??基因表達(dá)定量完成后,需要對其表達(dá)數(shù)據(jù)進(jìn)行統(tǒng)計學(xué)分析,篩選樣本在不同狀態(tài)下表達(dá)水平顯著差異的基因。差異分析主要分為三個步驟。
首先對原始的readcount進(jìn)行標(biāo)準(zhǔn)化(normalization),主要是對測序深度的校正。
然后統(tǒng)計學(xué)模型進(jìn)行假設(shè)檢驗概率(pvalue)的計算
最后進(jìn)行多重假設(shè)檢驗校正,得到FDR值(錯誤發(fā)現(xiàn)率,padj是其常見形式)[1-2]。
??針對不同的實驗情況,我們選用合適的軟件進(jìn)行基因表達(dá)差異顯著性分析,具體如下表所示。
表1 表達(dá)差異分析所用軟件及差異基因篩選標(biāo)準(zhǔn)
類型 | 軟件 | 標(biāo)準(zhǔn)化方法 | pvalue計算模型 | FDR計算方法 | 差異基因篩選標(biāo)準(zhǔn) |
---|---|---|---|---|---|
有生物學(xué)重復(fù) | DESeq2(Anders et al, 2014) | DESeq | 負(fù)二項分布 | BH | |log2(FoldChange)| > 0 & padj < 0.05 |
無生物學(xué)重復(fù) | edgeR(Robinson et al, 2010) | TMM | 負(fù)二項分布 | BH | |log2(FoldChange)| > 1 & padj < 0.05 |
??若按照以上標(biāo)準(zhǔn)篩選得到的差異基因過少(低于100),很有可能導(dǎo)致后面的功能富集分析沒有顯著性結(jié)果,所以,我們會根據(jù)項目的具體情況,適當(dāng)?shù)亟档秃Y選差異基因的閾值標(biāo)準(zhǔn)。若項目實驗只關(guān)注某幾個基因的表達(dá)情況(如基因敲除),不在意富集結(jié)果,從下面的差異分析表格中篩選關(guān)注的那幾個基因即可。
??一般來說,如果一個基因在兩組樣品中的表達(dá)量差異達(dá)到兩倍以上,我們認(rèn)為這樣的基因是具有表達(dá)差異的。為了判斷兩個樣品之間的表達(dá)量差異究竟是由于各種誤差導(dǎo)致的還是本質(zhì)差異,我們需要對所有基因在這兩個樣本中的表達(dá)量數(shù)據(jù)進(jìn)行假設(shè)檢驗。而轉(zhuǎn)錄組分析是針對成千上萬個基因進(jìn)行的,這樣會導(dǎo)致假陽性的累積,基因數(shù)目越多,假設(shè)檢驗的假陽性累積程度會越高,所以引入padj對假設(shè)檢驗的P-value進(jìn)行校正,從而控制假陽性的比例[3]。
??差異基因的篩選標(biāo)準(zhǔn)是非常重要的,我們給出的標(biāo)準(zhǔn)|log2(FoldChange)| > 1 & padj< 0.05是常用的經(jīng)驗值,在實際項目中可以根據(jù)情況靈活選擇。例如,差異倍數(shù)可以選擇1.5倍,也可以選擇3倍,padj常用的閾值包括0.01、0.05、0.1等。若按照以上標(biāo)準(zhǔn)篩選得到的差異基因過少,很有可能導(dǎo)致后?的功能富集分析沒有顯著性結(jié)果。若項目實驗只關(guān)注某幾個基因的表達(dá)情況(如基因敲除),不在意富集結(jié)果,從下面的差異分析表格中篩選關(guān)注的那幾個基因即可。反之,如果得到的差異基因數(shù)目過多,不利于后續(xù)目標(biāo)基因的篩選,這個時候可使用更嚴(yán)格的閾值標(biāo)準(zhǔn)進(jìn)行篩選,則可以使用更嚴(yán)格的閾值標(biāo)準(zhǔn)進(jìn)行篩選。
2.6.1. 差異miRNA的篩選
??通過Deseq2進(jìn)行差異分析,我們通常采用 |log2FC|>1 & padj < 0.05
進(jìn)行差異miRNA的篩選,差異計算結(jié)果如下。
結(jié)果文件:
結(jié)果說明 | 結(jié)果文件 |
---|---|
定量分析表達(dá)矩陣結(jié)果 | 04.DE/miR.Counts.csv |
差異miRNA分析結(jié)果(ALL) | 04.DE/Allgene_anno_ALL.xls |
差異miRNA分析結(jié)果(篩選后) | 04.DE/Allgene_anno.xls |
04.DE/Allgene_anno.xls
表頭
表頭 | 說明 |
---|---|
ENSEMBL | 差異miRNA的ENSEMBL名 |
pvalue | 差異miRNA的置信度計算結(jié)果 |
padj | 差異miRNA的多重校驗FDR |
log2FC | Treat組 vs Control組 差異倍數(shù) 的log2標(biāo)準(zhǔn)化結(jié)果 |
FC | Treat組 vs Control組 差異倍數(shù) |
log2FC_abs | Treat組 vs Control組 差異倍數(shù) 的log2標(biāo)準(zhǔn)化結(jié)果的絕對值(此列便于篩選log2FC閾值) |
FC_HvsL | 高表達(dá)組 vs 低表達(dá)組 差異倍數(shù) (此列便于篩選FC閾值) |
change | 使用本次分析的閾值,對差異miRNA的上下調(diào)標(biāo)記 |
miRNA | 差異miRNA的miRNA名 |
ENTREZID | 差異miRNA的ENTREZID號 |
GENENAME | 差異miRNA的基本描述信息 |
baseMean | 差異miRNA的表達(dá)量標(biāo)準(zhǔn)化后的平均值 |
Samples* | 樣本的原始表達(dá)矩陣表達(dá)量結(jié)果 |
Samples*_normal | 樣本的表達(dá)矩陣標(biāo)準(zhǔn)化后的結(jié)果 |
2.6.2. 差異miRNA的熱圖聚類
??將所有比較組的差異miRNA取并集之后作為差異miRNA集。兩組以上的實驗,可對差異miRNA集進(jìn)行聚類分析,將表達(dá)模式相近的基因聚在一起。我們采用主流的層次聚類對基因的表達(dá)值進(jìn)行聚類分析,對行(row)進(jìn)行均一化處理(Z-score)。熱圖中表達(dá)模式相近的基因或樣本會被聚集在一起,每個方格中的顏色反映的不是基因表達(dá)值,而是表達(dá)數(shù)據(jù)的行進(jìn)行均一化處理后得到的數(shù)值(一般在-1到1之間),所以熱圖中的顏色只能橫向比較(同一基因在不同樣本中的表達(dá)情況),不能縱向比較(同一樣本不同基因的表達(dá)情況)。結(jié)果文件中既有組間的聚類,也有樣品間的聚類。結(jié)題報告展示了樣品間的聚類,具體如下圖所示。
圖 2.6.2. 差異表達(dá)基因聚類熱圖
圖中橫坐標(biāo)為樣品名,縱坐標(biāo)為差異miRNA歸一化后的數(shù)值,顏色越紅,表達(dá)量越高,越藍(lán),表達(dá)量越低。
結(jié)果文件:
差異miRNA的熱圖結(jié)果:
04.DE/heatmap/
2.6.3. 差異miRNA的火山圖分布
??火山圖可直觀展示每個比較組合的差異miRNA分布情況,如下圖所示。圖中橫坐標(biāo)表示基因在處理和對照兩組中的表達(dá)倍數(shù)變化(log2FoldChange),縱坐標(biāo)表示基因在處理和對照兩組中表達(dá)差異的顯著性水平(-log10padj或-log10pvalue)。為上調(diào)基因用紅色點表示,下調(diào)基因用藍(lán)色點表示。
圖 2.6.3. 差異miRNA火山圖
圖中橫坐標(biāo)為log2FoldChange值,縱坐標(biāo)為-log10padj或-log10pvalue,藍(lán)色的虛線表示差異基因篩選標(biāo)準(zhǔn)的閾值線
結(jié)果文件:
差異基因的火山圖結(jié)果:
04.DE/volcano/
2.7. 靶基因預(yù)測
??miRNA是一類在生物體內(nèi)起到重要調(diào)控作用的的小片段非編碼RNA。一般認(rèn)為miRNA通過和mRNA的結(jié)合,可以抑制mRNA的表達(dá),從而影響到Gene的表達(dá)。
2.7.1. 靶基因預(yù)測結(jié)果
??目前miRNA的靶基因預(yù)測軟件主要是依據(jù)miRNA與其靶位點的互補性、miRNA與其靶位點的調(diào)控關(guān)系、miRNA靶位點的保守性、miRNA-mRNA雙鏈之間的熱穩(wěn)定性及附近序列的二級結(jié)構(gòu)等原則設(shè)計的。 目前miRNA預(yù)測主要基于算法預(yù)測,目前也有一部分以實驗結(jié)果為收集對象的數(shù)據(jù)庫。使用算法預(yù)測的數(shù)據(jù)庫,有如 TargetScan、miRDB、RNAhybrid、DIANA-microT-CDS 等。以實驗實驗驗證的收集的靶標(biāo)相互作用數(shù)據(jù)庫,有如 miRTarBase、Tarbase v8 等。 考慮到研究時效性,我們選擇基于miRbase22為基礎(chǔ)研究的數(shù)據(jù)庫,在此基礎(chǔ)上篩選出近5年內(nèi)較新的數(shù)據(jù)庫。
??其中,miRDB
數(shù)據(jù)庫 是根據(jù) “miRNA通過下調(diào)其靶基因的表達(dá)” miRNA主要功能特征,采用機(jī)器學(xué)習(xí)來進(jìn)行模型預(yù)測得到的靶標(biāo)數(shù)據(jù)庫。它通過采用大量RNA數(shù)據(jù)集,使用 “過表達(dá)的miRNA“ 及 “下調(diào)的轉(zhuǎn)錄本” 進(jìn)一步量化miRNA靶向下調(diào)的特征,采用SVM算法開發(fā)了預(yù)測模型 MirTarget
。 該算法結(jié)合實驗數(shù)據(jù)同時比較 TargetScan
,DIANA-MicroT
,miRanda
和 PITA
另外四種預(yù)測算法,獲得了最優(yōu)表現(xiàn)。 該模型對于每個候選靶標(biāo)基因,它都會生成由SVM計算出的概率分?jǐn)?shù),該分?jǐn)?shù)反映了預(yù)測結(jié)果的統(tǒng)計置信度,分?jǐn)?shù)越高,置信度越高。所有預(yù)測目標(biāo)的目標(biāo)預(yù)測分?jǐn)?shù)在50到100之間,一般而言,預(yù)測得分 > 80 的預(yù)測目標(biāo)具有較高的可信度。如果分?jǐn)?shù)低于60,則需要謹(jǐn)慎選擇,一般需要提供其他證據(jù)支持。
??TargetScan
數(shù)據(jù)庫 則基于序列互補原則,找到比對到靶 3'UTR 的保守性 8 mer、7 mer 或 6 mer 位點(seed match 序列),進(jìn)一步根據(jù)熱力學(xué)穩(wěn)定性篩選得到 miRNA 的靶基因。 隨后,對于不具備保守性的 seed match 區(qū)域計算相應(yīng)的 context score
。 將保守和不保守區(qū)域的 context score
進(jìn)行排序即得到 context score percentile
。 一般考慮 context score percentile
> 90 為預(yù)測的可能具有功能的 miRNA 的靶基因。
??綜上,miRDB數(shù)據(jù)庫主要通過miRNA-mRNA調(diào)控關(guān)系來進(jìn)行量化計算,而TargetScan通過序列互補原則篩選后再進(jìn)行熱穩(wěn)定性量化計算,前者相較于后者更有說服力,而后者得到的數(shù)據(jù)較為全面。我們綜合兩者,通過標(biāo)準(zhǔn) miRDB Score >= 80 & TargetScan context score percentile >= 90
的靶基因作為預(yù)測結(jié)果。通過預(yù)測給出的miRNA-靶基因?qū)?,用于后續(xù)進(jìn)一步網(wǎng)絡(luò)圖或富集分析。
圖 2.7.1. TargetScan vs miRDB 的靶基因交集veen圖
結(jié)果文件:
miRDB/TargetScan預(yù)測結(jié)果:
結(jié)果說明 結(jié)果文件 miRDB預(yù)測結(jié)果(原始) 05.Target/*/Res_miRDB/Res_*.xls
miRDB預(yù)測結(jié)果( miRDB Score >= 80
過濾后)05.Target/*/Res_miRDB/Res_*.sig.xls
miRDB預(yù)測結(jié)果(過濾情況統(tǒng)計) 05.Target/*/Res_miRDB/Res_*.Summary.xls
TargetScan預(yù)測結(jié)果(原始) 05.Target/*/Res_TargetScan/Res_*.xls
TargetScan預(yù)測結(jié)果( TargetScan score >= 90
過濾后)05.Target/*/Res_TargetScan/Res_*.sig.xls
TargetScan預(yù)測結(jié)果(過濾情況統(tǒng)計) 05.Target/*/Res_TargetScan/Res_*.Summary.xls
miRDB/TargetScan預(yù)測結(jié)果過濾后的交集結(jié)果:
結(jié)果說明 結(jié)果文件 miRDB/TargetScan兩者過濾后預(yù)測結(jié)果交集統(tǒng)計(總表) 05.Target/*/Compare/Compare_diffMirna2_ALL.xls
miRDB/TargetScan兩者過濾后預(yù)測結(jié)果交集統(tǒng)計(僅共存的)
(當(dāng)兩者中出現(xiàn)僅有一個軟件有結(jié)果時,則按一個軟件篩選)05.Target/*/Compare/Compare_diffMirna2.xls
僅共存列表的各個miRNA的靶基因數(shù)量統(tǒng)計 05.Target/*/Compare/Compare_Summary_miRNA.xls
僅共存列表的各個靶基因的相應(yīng)miRNA數(shù)量統(tǒng)計 05.Target/*/Compare/Compare_Summary_SYMBOL.xls
交集veen統(tǒng)計文件夾 05.Target/*/Compare/veen
以上結(jié)果均位于文件夾: 05.Target/
Results/05.Target/Compare/veen/
交集veen統(tǒng)計文件夾說明
表頭 | 說明 |
---|---|
所有基因 | |
??比較組A和B的所有基因 | compare/union0_all_* |
??比較組A的基因 | compare/union1_A_* |
??比較組B的基因 | compare/union2_B_* |
共有基因 | |
??比較組的共有基因 | compare/com_all_* |
特有基因 | |
??比較組A和B的特有基因 | compare/diff0_all_* |
??比較組A的特有基因 | compare/diff1_A_* |
??比較組B的特有基因 | compare/diff2_B_* |
比較組A和B的veen圖及相關(guān)說明 | compare/venn_gene_* |
2.8. 富集分析
??我們根據(jù)基因表達(dá)量分析得到差異基因之后,必須進(jìn)一步落到基因的功能上來。對于轉(zhuǎn)錄組分析而言,往往涉及到成千上萬個基因,這會使分析變得很復(fù)雜。解決思路是將一個基因列表分成多個部分,從而減少分析的復(fù)雜度。為了解決怎么分成不同類,通常會對基因功能進(jìn)行富集分析, 期望發(fā)現(xiàn)在生物學(xué)過程中起關(guān)鍵作用的生物通路, 從而揭示和理解生物學(xué)過程的基本分子機(jī)制。功能富集分析可以將成百上千個基因、蛋白或者其他分子分到不同的通路中,以減少分析的復(fù)雜度。另外,在兩種不同實驗條件下,激活的通路顯然比簡單的基因或蛋白列表更有說服力?;蚬δ芨患治鍪紫纫獦?gòu)建基因集(gene set,如GO和KEGG數(shù)據(jù)庫等),也就是基因組注釋信息進(jìn)行分類。然后再把我們的目標(biāo)基因集(差異基因集或者其他基因集)映射到背景基因集上,注意區(qū)分注釋與富集。
??我們采用clusterProfiler軟件對差異基因集進(jìn)行GO功能富集分析,KEGG通路富集分析等。富集分析基于超幾何分布原理,其中差異基因集為差異顯著分析所得差異基因并注釋到GO或KEGG數(shù)據(jù)庫的基因集,背景基因集為所有進(jìn)行差異顯著分析的基因并注釋到GO或KEGG數(shù)據(jù)庫的基因集。富集分析結(jié)果是對每個差異比較組合的所有差異基因集、上調(diào)差異基因集、下調(diào)差異基因集進(jìn)行富集。本報告中展示的表格是選取某一個比較組合的富集分析結(jié)果,圖片是所有組合的富集分析結(jié)果。
圖 2.8.1. 基因富集分析原理圖
2.8.1. 富集分析結(jié)果文件
結(jié)果說明 | 結(jié)果路徑 |
---|---|
GO富集分析結(jié)果 | |
GO富集結(jié)果列表(所有結(jié)果) | 06.ENRICH/*/gene.ego_all-p.adjust1.00.csv |
GO富集結(jié)果列表(按p.adj<0.05篩選后) | 06.ENRICH/*/gene.ego_all-p.adjust0.05.csv |
GO富集結(jié)果列表(MF、BP、CC所有結(jié)果) | 06.ENRICH/*/gene.ego_ALL.csv |
GO富集分析柱狀圖 | 06.ENRICH/*/gene.GO-*-barplot.p* |
GO富集分析散點圖 | 06.ENRICH/*/gene.GO-*-dotplot.p* |
GO富集分析DAG圖 | 06.ENRICH/*/gene.GO-*-DAG.p* |
KEGG富集分析結(jié)果 | |
KEGG富集結(jié)果列表(所有) | 06.ENRICH/*/gene.KEGG.csv |
KEGG富集結(jié)果列表(按p.adj<0.05篩選后) | 06.ENRICH/*/gene.KEGG_significant.csv |
KEGG富集分析柱狀圖 | 06.ENRICH/*/gene.KEGG-*-barplot.p* |
KEGG富集分析散點圖 | 06.ENRICH/*/gene.KEGG-*-dotplot.p* |
結(jié)果文件夾:
分析結(jié)果文件夾:
06.ENRICH/
以上富集分析各圖的可視化網(wǎng)頁:
all-pdf.html
、up-pdf.html
、down-pdf.html
表頭說明: (Results/*enrich_*/gene.ego_*.csv
GO富集結(jié)果列表)
ID | 對應(yīng)GO數(shù)據(jù)庫中的ID |
---|---|
ONTOLOGY | 分子功能(Molecular Function),生物過程(biological process)和細(xì)胞組成(cellular component) |
Description | GO的描述 |
GeneRatio | 對應(yīng)GO 差異基因數(shù) / 能夠?qū)?yīng)到GO數(shù)據(jù)庫中同類型的差異基因數(shù) |
BgRatio | 對應(yīng)GO包含對應(yīng)物種的基因數(shù) / GO數(shù)據(jù)庫中包含對應(yīng)物種的基因數(shù) |
pvalue | 富集分析得到的p-value |
p.adjust | 校正后的p-value |
qvalue | 富集分析得到的qvalue |
Count | 富集基因數(shù)目 |
ENTREZID | 富集基因列表(ENTREZID) |
SYMBOL | 富集基因列表(SYMBOL) |
表頭說明: (Results/*enrich_*/gene.KEGG*.csv
KEGG富集結(jié)果列表)
ID | 對應(yīng)PATHWAY數(shù)據(jù)庫中的ID |
---|---|
Description | PATHWAY的描述 |
GeneRatio | 對應(yīng)PATHWAY 差異基因數(shù) / 能夠?qū)?yīng)到PATHWAY數(shù)據(jù)庫中的差異基因數(shù) |
BgRatio | 對應(yīng)PATHWAY包含對應(yīng)物種的基因數(shù) / PATHWAY數(shù)據(jù)庫中包含對應(yīng)物種的基因數(shù) |
pvalue | 富集分析得到的p-value |
p.adjust | 校正后的p-value |
qvalue | 富集分析得到的qvalue |
Count | 富集基因數(shù)目 |
ENTREZID | 富集基因列表(ENTREZID) |
SYMBOL | 富集基因列表(SYMBOL) |
2.8.2. GO功能富集分析
??GO(Gene Ontology)是描述基因功能的綜合性數(shù)據(jù)庫,可分為生物過程(biological process)和細(xì)胞組成(cellular component)分子功能(Molecular Function)三個部分。GO功能富集以padj小于0.05作為為顯著性富集的閾值,富集結(jié)果見結(jié)果文件。
??從GO富集分析結(jié)果中,選取最顯著的30個Term繪制柱狀圖進(jìn)行展示,若不足30個,則繪制所有Term,按生物過程、細(xì)胞組分和分子功能三大類別及差異基因上下調(diào)分類畫的柱狀圖。
??有向無環(huán)圖(Directed Acyclic Graph,DAG)為差異基因GO富集分析結(jié)果的圖形化展示方式。圖中,分支代表包含關(guān)系,從上至下所定義的功能范圍越來越小,選取每個差異比較組合的GO富集結(jié)果最顯著性前5位的GO Term作為有向無環(huán)圖的主節(jié)點,并通過包含關(guān)系,將相關(guān)聯(lián)的GO Term一起展示,顏色的深淺代表富集程度。我們的項目中分別繪制生物過程、分子功能和細(xì)胞組分的DAG圖。
圖 2.8.2.1. GO富集分析柱狀圖
圖中縱坐標(biāo)為GO Term,橫坐標(biāo)為GO Term富集的顯著性水平,數(shù)值越高越顯著
圖 2.8.2.2. GO富集分析散點圖
圖中橫坐標(biāo)為注釋到GO Term上的差異基因數(shù)與差異基因總數(shù)的比值,縱坐標(biāo)為GO Term
圖 2.8.2.3. GO富集分析DAG圖
每個節(jié)點代表一個GO術(shù)語,方框代表的是富集程度為TOP5的GO,顏色的深淺代表富集程度,顏色越深就表示富集程度越高,每個節(jié)點上展示了該TERM的名稱及富集分析的padj
2.8.3. KEGG通路富集分析
??KEGG(Kyoto Encyclopedia of Genes and Genomes)是整合了基因組、化學(xué)和系統(tǒng)功能信息的綜合性數(shù)據(jù)庫。KEGG通路富集以padj小于0.05作為顯著性富集的閾值,富集結(jié)果見結(jié)果文件。
??從KEGG富集結(jié)果中,選取最顯著的20個KEGG通路繪制柱狀圖進(jìn)行展示,若不足20個,則繪制所有通路,如下圖所示。圖中橫坐標(biāo)為通路富集的顯著性水平,數(shù)值越高越顯著,縱坐標(biāo)為KEGG通路。
??從KEGG富集結(jié)果中,選取最顯著的20個KEGG通路繪制散點圖進(jìn)行展示,若不足20個,則繪制所有通路,如下圖所示。圖中橫坐標(biāo)為注釋到KEGG通路上的差異基因數(shù)與差異基因總數(shù)的比值,縱坐標(biāo)為KEGG通路,點的大小代表注釋到KEGG通路上的基因數(shù),顏色從紅到紫代表富集的顯著性大小。
圖 2.8.3.1. KEGG富集分析柱狀圖
圖中橫坐標(biāo)為通路富集的顯著性水平,數(shù)值越高越顯著,縱坐標(biāo)為KEGG通路。
圖 2.8.3.2. KEGG富集散點圖
圖中橫坐標(biāo)為注釋到KEGG通路上的差異基因數(shù)與差異基因總數(shù)的比值,縱坐標(biāo)為KEGG通路