技術(shù)服務(wù)描述

miRNA結(jié)題報(bào)告

生信部

2021年01月28日

1. 建庫測序工作流程

1.1. 建庫測序流程

??MicroRNA (miRNA) 是一類長度約為20-24個(gè)核苷酸長度的具有調(diào)控功能的非編碼RNA，miRNA 是生物體內(nèi)一類重要的特殊分子，誘導(dǎo)基因沉默，參與細(xì)胞生長、發(fā)育、基因轉(zhuǎn)錄和翻譯等諸多生命活動(dòng)的調(diào)控過程。

??從RNA樣品到最終分析結(jié)果數(shù)據(jù)，需要經(jīng)過樣品檢測、建庫、測序和信息分析等過程，其中測序過程我們采用高通量測序illumina HiSeq/MiSeq測序平臺(tái)；illumina測序得到的原始圖像數(shù)據(jù)文件經(jīng)堿基識(shí)別（Base Calling）分析轉(zhuǎn)化為原始測序序列（Sequenced Reads），我們稱之為Raw Data或Raw Reads，結(jié)果以FASTQ（簡稱為fq）文件格式存儲(chǔ)，其中包含測序序列（reads）的序列信息以及其對(duì)應(yīng)的測序質(zhì)量信息。

??樣品質(zhì)量檢測：利用NanoDrop 2000分光光度計(jì)對(duì)RNA樣品進(jìn)行初步定量，Aglient 2100對(duì)樣品濃度精確定量。RNA樣品總量、濃度、完整性（RIN值，RNA Integrity Number）、純度（OD260/OD280比值）符合收樣立項(xiàng)標(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方法對(duì)文庫進(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)計(jì)序列長度分布等過程完成初級(jí)分析；將初級(jí)分析得到的clean reads與ncRNA庫比對(duì)，分類注釋ncRNA，并去除reads中的rNRA、tRNA、snRNA、snoRNA等；然后再將reads比對(duì)到參考基因組、miRBase上，計(jì)算miRNA表達(dá)量，注釋已知miRNA、預(yù)測新的miRNA；然后進(jìn)行差異miRNA、靶基因預(yù)測、功能注釋等后續(xù)分析。

miRNA生物信息分析流程圖

2.1. miRNA-seq 測序數(shù)據(jù)質(zhì)量評(píng)估

2.1.1. 數(shù)據(jù)質(zhì)量評(píng)估 及 數(shù)據(jù)過濾

??對(duì)于測序得到的rawdata，我們首先使用fastqc進(jìn)行數(shù)據(jù)質(zhì)量評(píng)估，通過結(jié)果，我們可以了解到測序數(shù)據(jù)reads質(zhì)量、長度、堿基分布，序列重復(fù)頻次等基本信息。

??測序得到的raw data，里面含有帶接頭的、低質(zhì)量的reads，為了保證后續(xù)信息分析的質(zhì)量，必須對(duì)raw reads進(jìn)行處理，得到clean reads。

??raw data的過濾步驟如下： 去除低質(zhì)量 reads; 去除有5’接頭污染的reads; 去除沒有3’接頭序列; 去除沒有插入片段的reads; 去除長度小于16nt的reads;

??數(shù)據(jù)經(jīng)過以上過濾后，我們對(duì)數(shù)據(jù)做一些基本統(tǒng)計(jì)。

結(jié)果：

2.2. 數(shù)據(jù)長度分布

??一般來說，小RNA的長度區(qū)間為18～30nt，數(shù)據(jù)過濾完成后，我們統(tǒng)計(jì)了clean reads數(shù)目及各自占總reads的百分比，并繪制reads長度分布圖。reads長度分布圖能幫助我們判斷小RNA的種類，如miRNA集中在21或22nt，siRNA集中在24nt，piRNA集中在30nt。

Clean Data Reads 長度分布統(tǒng)計(jì)直方圖：

Figure 2.2: Small RNA Clean Reads長度分布圖

Clean Reads長度分布直方圖，橫坐標(biāo)為Reads的長度，縱坐標(biāo)為樣品中對(duì)應(yīng)該長度的total Reads數(shù)量。

2.3. miRNA基因組比對(duì)分析

??我們將Clean Reads 與 mature miRNA, miRNA hairpin, mRNA, mature & primary tRNA, snoRNA, rRNA, other non-coding RNA, and (optional) known RNA 使用bowtie進(jìn)行比對(duì)搜索，然后統(tǒng)計(jì)reads的比對(duì)情況，分析樣本中各類ncRNA的數(shù)目及比例情況，統(tǒng)計(jì)結(jié)果如下。過濾出待分析的miRNA以便后續(xù)分析，將clean reads中ncRNA盡可能地去除，便于后續(xù) miRNA 檢測及預(yù)測。

Reads類型分布圖

圖顯示的是各類型的RNA reads和miRNA reads占總distinct reads的比例。

與參考基因組比對(duì)質(zhì)控結(jié)果：

2.4. 已知 miRNA 的檢測以及 novel miRNA 預(yù)測分析

??首先用去除了ncRNA后的Reads比對(duì)miRbase數(shù)據(jù)庫，比對(duì)上miRbase中 reads 用于檢測已知成熟的miRNA，并統(tǒng)計(jì)其表達(dá)量和RPM值。我們采用miRge軟件來進(jìn)行 miRNA reads 檢測注釋、novel miRNA的預(yù)測、miRNA表達(dá)量的統(tǒng)計(jì)。

miRNA分析結(jié)果匯總：

2.4.1. miRNA異構(gòu)體分析

??隨著高通量測序技術(shù)的發(fā)展，早期認(rèn)為的miRNA loci僅產(chǎn)生1條miRNA成熟體序列這一觀點(diǎn)被顛覆。對(duì)于某一miRNA來說，它并不是單一的序列，而是由一系列長度/序列及表達(dá)不同的異構(gòu)體 (isomiRs) 組成。這些isomiRs表達(dá)多樣且序列多樣，甚至引入多樣的5'端及種子區(qū)域。特定miRNA位點(diǎn)在疾病組織中可具有異常的表達(dá)模式，現(xiàn)已證實(shí)，部分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比對(duì)的情況，將各個(gè)miRNA Reads進(jìn)行isomiR識(shí)別，并同時(shí)匯總至miRNA reads 注釋結(jié)果中，分類統(tǒng)計(jì)，并選取表達(dá)量較高的前20個(gè)mirna進(jìn)行各類IsomiRs表達(dá)豐度統(tǒng)計(jì)。

圖2.4.1 所有樣本的變體類型分布（isomiRs）

圖2.4.2 各個(gè)樣本的變體分布TOP20的豐度統(tǒng)計(jì)

2.4.2. miRNA 堿基編輯

??成熟miRNA序列的第2-8個(gè)堿基被稱作“種子”序列，保守性很高。若在這一區(qū)域發(fā)生堿基突變，則可能改變miRNA的靶基因作用位點(diǎn)。我們通過將miRNA序列與miRBase中已知miRNA前體以及成熟miRNA序列進(jìn)行比對(duì)(只允許一個(gè)位點(diǎn)的錯(cuò)配)找出發(fā)生了堿基突變的miRNA，統(tǒng)計(jì)匯總在該miRNA的堿基突變位點(diǎn)，以及發(fā)生堿基編輯的Reads數(shù)量及百分比，并對(duì)統(tǒng)計(jì)概況結(jié)果進(jìn)行可視化。

圖2.4.3 鑒定為已存在mirna堿基編輯的在各個(gè)樣本中的相應(yīng)readCount占比統(tǒng)計(jì)圖

miRNA堿基編輯分析結(jié)果詳細(xì)說明：

1. A2I在各個(gè)樣本中統(tǒng)計(jì)匯總表

• 結(jié)果文件：02.miRNA/2.a2IEditing/a2IEditing.report.csv

• 表頭說明：

  結(jié)果說明	結(jié)果文件
  基本注釋結(jié)果：
  miRNA reads 變體類型注釋結(jié)果	02.miRNA/1.isomiR/sample_miRge3.gff
  miRNA reads 變體類型可視化統(tǒng)計(jì)	02.miRNA/1.isomiR/*svg
  miRNA堿基編輯分析結(jié)果：
  A2I在各個(gè)樣本中統(tǒng)計(jì)匯總表	02.miRNA/2.a2IEditing/a2IEditing.report.csv
  A2I在各個(gè)樣本中統(tǒng)計(jì)匯總表的篩選及簡化信息 ( 即匯總表的*.RPM.mismatch、*.AtoI.percentage列 )	02.miRNA/2.a2IEditing/a2IEditing.report.newform.csv
  A2I分析的序列及統(tǒng)計(jì)詳情信息	02.miRNA/2.a2IEditing/a2IEditing.detail.txt
  二級(jí)結(jié)構(gòu)預(yù)測結(jié)果：
  novel miRNA 二級(jí)結(jié)構(gòu)預(yù)測結(jié)果	02.miRNA/3.novel_predict/*.R1_novel_miRNAs/
  novel miRNA 二級(jí)結(jié)構(gòu)預(yù)測詳細(xì)信息表格匯總	02.miRNA/3.novel_predict/*.R1_novel_miRNAs/*_novel_miRNAs_report.csv

2. A2I在各個(gè)樣本中統(tǒng)計(jì)匯總表的篩選及簡化信息
   ( 即匯總表的*.RPM.mismatch、*.AtoI.percentage列，該表對(duì)應(yīng)上面的可視化結(jié)果 )

• 結(jié)果文件：02.miRNA/2.a2IEditing/a2IEditing.report.newform.csv

• 表頭說明：




結(jié)果說明	結(jié)果文件
基本注釋結(jié)果：
miRNA reads 變體類型注釋結(jié)果	02.miRNA/1.isomiR/sample_miRge3.gff
miRNA reads 變體類型可視化統(tǒng)計(jì)	02.miRNA/1.isomiR/*svg
miRNA堿基編輯分析結(jié)果：
A2I在各個(gè)樣本中統(tǒng)計(jì)匯總表	02.miRNA/2.a2IEditing/a2IEditing.report.csv
A2I在各個(gè)樣本中統(tǒng)計(jì)匯總表的篩選及簡化信息 ( 即匯總表的*.RPM.mismatch、*.AtoI.percentage列 )	02.miRNA/2.a2IEditing/a2IEditing.report.newform.csv
A2I分析的序列及統(tǒng)計(jì)詳情信息	02.miRNA/2.a2IEditing/a2IEditing.detail.txt
二級(jí)結(jié)構(gòu)預(yù)測結(jié)果：
novel miRNA 二級(jí)結(jié)構(gòu)預(yù)測結(jié)果	02.miRNA/3.novel_predict/*.R1_novel_miRNAs/
novel miRNA 二級(jí)結(jié)構(gòu)預(yù)測詳細(xì)信息表格匯總	02.miRNA/3.novel_predict/*.R1_novel_miRNAs/*_novel_miRNAs_report.csv




3. A2I分析的序列及統(tǒng)計(jì)詳情信息

• 結(jié)果文件：02.miRNA/2.a2IEditing/a2IEditing.detail.txt

• 結(jié)果說明：

  ①每一個(gè)miRNA塊，包含3小塊內(nèi)容：原始序列、篩選后的序列統(tǒng)計(jì)、篩選后的序列。重復(fù)的miRNA塊代表不同樣本，順序?yàn)?A2I在各個(gè)樣本中統(tǒng)計(jì)匯總表 的樣本信息。

  ②每一行序列包含3個(gè)信息， miRNA sequence（reads序列），*.readCount（Count計(jì)數(shù)），Filter（是否篩選到后續(xù)進(jìn)行堿基編輯分析）

2.4.3. novel miRNA 預(yù)測及二級(jí)結(jié)構(gòu)分析

??miRge根據(jù)序列特征，堿基配對(duì)原則，最小自由能等特征模型，進(jìn)行新穎的miRNA及二級(jí)結(jié)構(gòu)預(yù)測。預(yù)測示意圖如下。結(jié)果見：02.miRNA/3.novel_predict/

圖2.4.4 新miRNA的二級(jí)結(jié)構(gòu)預(yù)測示意圖

2.5. miRNA 定量分析

2.5.1. 樣本間相關(guān)性分析

??生物學(xué)重復(fù)通常是任何生物學(xué)實(shí)驗(yàn)所必須的，目前主流期刊也基本要求生物學(xué)重復(fù)。生物學(xué)重復(fù)主要有兩個(gè)用途：一個(gè)是證明所涉及的生物學(xué)實(shí)驗(yàn)操作不是偶然，而是可重復(fù)的。另一個(gè)是為了確保后續(xù)的差異分析得到更可靠的結(jié)果。樣品間基因表達(dá)水平相關(guān)性是檢驗(yàn)實(shí)驗(yàn)可靠性和樣本選擇是否合理的重要指標(biāo)。相關(guān)系數(shù)越接近1，表明樣品之間表達(dá)模式的相似度越高。Encode計(jì)劃建議皮爾遜相關(guān)系數(shù)的平方(R2)大于0.92(理想的取樣和實(shí)驗(yàn)條件下)。具體的項(xiàng)目操作中，我們要求生物學(xué)重復(fù)樣品間R2至少要大于0.8，否則需要對(duì)樣品做出合適的解釋，或者重新進(jìn)行實(shí)驗(yàn)。根據(jù)各樣本所有基因的表達(dá)值計(jì)算組內(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）也常用來評(píng)估組間差異及組內(nèi)樣本重復(fù)情況，PCA采用線性代數(shù)的計(jì)算方法，對(duì)數(shù)以萬計(jì)的基因變量進(jìn)行降維及主成分提取。我們對(duì)所有樣本的基因表達(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á)定量完成后，需要對(duì)其表達(dá)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)學(xué)分析，篩選樣本在不同狀態(tài)下表達(dá)水平顯著差異的基因。差異分析主要分為三個(gè)步驟。

• 首先對(duì)原始的readcount進(jìn)行標(biāo)準(zhǔn)化（normalization），主要是對(duì)測序深度的校正。

• 然后統(tǒng)計(jì)學(xué)模型進(jìn)行假設(shè)檢驗(yàn)概率（pvalue）的計(jì)算

• 最后進(jìn)行多重假設(shè)檢驗(yàn)校正，得到FDR值（錯(cuò)誤發(fā)現(xiàn)率，padj是其常見形式)[1-2]。

??針對(duì)不同的實(shí)驗(yàn)情況，我們選用合適的軟件進(jìn)行基因表達(dá)差異顯著性分析，具體如下表所示。

??若按照以上標(biāo)準(zhǔn)篩選得到的差異基因過少（低于100），很有可能導(dǎo)致后面的功能富集分析沒有顯著性結(jié)果，所以，我們會(huì)根據(jù)項(xiàng)目的具體情況，適當(dāng)?shù)亟档秃Y選差異基因的閾值標(biāo)準(zhǔn)。若項(xiàng)目實(shí)驗(yàn)只關(guān)注某幾個(gè)基因的表達(dá)情況（如基因敲除），不在意富集結(jié)果，從下面的差異分析表格中篩選關(guān)注的那幾個(gè)基因即可。

??一般來說，如果一個(gè)基因在兩組樣品中的表達(dá)量差異達(dá)到兩倍以上，我們認(rèn)為這樣的基因是具有表達(dá)差異的。為了判斷兩個(gè)樣品之間的表達(dá)量差異究竟是由于各種誤差導(dǎo)致的還是本質(zhì)差異，我們需要對(duì)所有基因在這兩個(gè)樣本中的表達(dá)量數(shù)據(jù)進(jìn)行假設(shè)檢驗(yàn)。而轉(zhuǎn)錄組分析是針對(duì)成千上萬個(gè)基因進(jìn)行的，這樣會(huì)導(dǎo)致假陽性的累積，基因數(shù)目越多，假設(shè)檢驗(yàn)的假陽性累積程度會(huì)越高，所以引入padj對(duì)假設(shè)檢驗(yàn)的P-value進(jìn)行校正，從而控制假陽性的比例[3]。

??差異基因的篩選標(biāo)準(zhǔn)是非常重要的，我們給出的標(biāo)準(zhǔn)|log2(FoldChange)| > 1 & padj< 0.05是常用的經(jīng)驗(yàn)值，在實(shí)際項(xiàng)目中可以根據(jù)情況靈活選擇。例如，差異倍數(shù)可以選擇1.5倍，也可以選擇3倍，padj常用的閾值包括0.01、0.05、0.1等。若按照以上標(biāo)準(zhǔn)篩選得到的差異基因過少，很有可能導(dǎo)致后?的功能富集分析沒有顯著性結(jié)果。若項(xiàng)目實(shí)驗(yàn)只關(guān)注某幾個(gè)基因的表達(dá)情況（如基因敲除），不在意富集結(jié)果，從下面的差異分析表格中篩選關(guān)注的那幾個(gè)基因即可。反之，如果得到的差異基因數(shù)目過多，不利于后續(xù)目標(biāo)基因的篩選，這個(gè)時(shí)候可使用更嚴(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的篩選，差異計(jì)算結(jié)果如下。

結(jié)果文件：

04.DE/Allgene_anno.xls表頭

2.6.2. 差異miRNA的熱圖聚類

??將所有比較組的差異miRNA取并集之后作為差異miRNA集。兩組以上的實(shí)驗(yàn)，可對(duì)差異miRNA集進(jìn)行聚類分析，將表達(dá)模式相近的基因聚在一起。我們采用主流的層次聚類對(duì)基因的表達(dá)值進(jìn)行聚類分析，對(duì)行（row）進(jìn)行均一化處理（Z-score）。熱圖中表達(dá)模式相近的基因或樣本會(huì)被聚集在一起，每個(gè)方格中的顏色反映的不是基因表達(dá)值，而是表達(dá)數(shù)據(jù)的行進(jìn)行均一化處理后得到的數(shù)值（一般在-1到1之間），所以熱圖中的顏色只能橫向比較（同一基因在不同樣本中的表達(dá)情況），不能縱向比較（同一樣本不同基因的表達(dá)情況）。結(jié)果文件中既有組間的聚類，也有樣品間的聚類。結(jié)題報(bào)告展示了樣品間的聚類，具體如下圖所示。

圖 2.6.2. 差異表達(dá)基因聚類熱圖

圖中橫坐標(biāo)為樣品名，縱坐標(biāo)為差異miRNA歸一化后的數(shù)值，顏色越紅，表達(dá)量越高，越藍(lán)，表達(dá)量越低。

結(jié)果文件：

差異miRNA的熱圖結(jié)果：04.DE/heatmap/

2.6.3. 差異miRNA的火山圖分布

??火山圖可直觀展示每個(gè)比較組合的差異miRNA分布情況，如下圖所示。圖中橫坐標(biāo)表示基因在處理和對(duì)照兩組中的表達(dá)倍數(shù)變化(log2FoldChange)，縱坐標(biāo)表示基因在處理和對(duì)照兩組中表達(dá)差異的顯著性水平(-log10padj或-log10pvalue)。為上調(diào)基因用紅色點(diǎn)表示，下調(diào)基因用藍(lán)色點(diǎ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á)。

all

down

up

圖 2.7.1. TargetScan vs miRDB 的靶基因交集veen圖