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

項(xiàng)目信息

合同編號(hào)：DEMO-2021-01-29-xx

客戶姓名：Client-name

客戶單位：Unit-address

1. 分析流程

1.1. 建庫測(cè)序流程

??從RNA樣品提取到最終數(shù)據(jù)獲得，樣品檢測(cè)、建庫、測(cè)序等每一環(huán)節(jié)都會(huì)直接影響數(shù)據(jù)的數(shù)量和質(zhì)量，從而影響后續(xù)數(shù)據(jù)分析的結(jié)果。為從源頭保證測(cè)序數(shù)據(jù)準(zhǔn)確可靠，在數(shù)據(jù)的所有生產(chǎn)環(huán)節(jié)都嚴(yán)格把關(guān)，從根源上確保高質(zhì)量數(shù)據(jù)的產(chǎn)出。建庫測(cè)序的流程：

1. Total RNA 樣本檢測(cè)

2. RNA 富集

3. 雙鏈cDNA合成

4. 末端修復(fù)、加A和接頭

5. 片段選擇和 PCR 擴(kuò)增

6. 文庫質(zhì)檢

7. Illumina測(cè)序

1.2. 信息分析流程

??RNA-seq的核心是基因表達(dá)差異的顯著性分析，使用統(tǒng)計(jì)學(xué)方法，比較兩個(gè)條件或多個(gè)條件下的基因表達(dá)差異，從中找出與條件相關(guān)的特異性基因，然后進(jìn)一步分析這些特異性基因的生物學(xué)意義，分析過程包括質(zhì)控、比對(duì)、定量、差異顯著性分析、功能富集等環(huán)節(jié)。信息分析流程如下圖所示：

2. 信息分析

2.2. 測(cè)序數(shù)據(jù)質(zhì)量控制

對(duì)原始測(cè)序數(shù)據(jù)及去除接頭后的可用數(shù)據(jù)進(jìn)行質(zhì)量評(píng)估。測(cè)序數(shù)據(jù)一般為雙端測(cè)序，因此，每個(gè)測(cè)序樣本會(huì)有兩個(gè)測(cè)序結(jié)果。

評(píng)估的具體內(nèi)容見：

RawData-fastqc 文件鏈接： /result/qc/qc_rawdata/*.html
CleanData-fastqc 文件鏈接： /result/qc/qc_cleandata/*.html
Fastqc 格式補(bǔ)充說明： /result/qc/qc_Supplement.html

2.3. 參考基因組比對(duì)

??測(cè)序片段（fragments）是mRNA隨機(jī)打斷的，為了確定這些一段由哪些基因轉(zhuǎn)錄來，需要將質(zhì)控后的clean reads比對(duì)到參考基因組上。使用HISAT2軟件將Clean Reads與參考基因組進(jìn)行快速精確的比對(duì)，獲取Reads在參考基因組上的定位信息^[4]。HISAT2軟件官方手冊(cè)。

??如果參考基因組組裝的較為完善，而且所測(cè)物種與參考基因組一致，且相關(guān)實(shí)驗(yàn)不存在污染，那么實(shí)驗(yàn)所產(chǎn)生的測(cè)序reads成功比對(duì)到基因組的比例會(huì)高于70% (Total Mapped Reads or Fragments)。本項(xiàng)目所用參考基因組為 hg38 ，下載鏈接：Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz 。

結(jié)果文件：

各個(gè)樣本的比對(duì)情況統(tǒng)計(jì)文件：/result/map_stat/*.flagstat.txt

2.4. 定量分析

2.4.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ù)的平方(R²)大于0.92(理想的取樣和實(shí)驗(yàn)條件下)。具體的項(xiàng)目操作中，我們要求生物學(xué)重復(fù)樣品間R²至少要大于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)性熱圖如下圖所示。

圖 1 樣本間相關(guān)性熱圖

圖中橫縱坐標(biāo)為各樣本相關(guān)系數(shù)的平方

結(jié)果文件：

樣本間相關(guān)性熱圖結(jié)果：Quant/cor_pheatmap*

2.4.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 主成分分析結(jié)果圖

圖中橫坐標(biāo)為第一主成分，縱坐標(biāo)為第二主成分

結(jié)果文件：

主成分分析結(jié)果：Quant/pca*

2.5. 差異分析

??基因表達(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ì)測(cè)序深度的校正。

• 然后統(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.5.1. 差異基因的篩選

??通過Deseq2進(jìn)行差異分析，我們通常采用 |log2FC|>1 & padj < 0.05 進(jìn)行差異基因的篩選，隨后對(duì)差異基因進(jìn)行注釋，得到包含注釋信息的差異基因列表。

結(jié)果文件：

差異基因列表及相關(guān)注釋信息（總的結(jié)果）：result/Enrichment/Allgene_anno_ALL.xls
差異基因列表及相關(guān)注釋信息（篩選結(jié)果）：result/Enrichment/Allgene_anno.xls

Differential/Allgene_anno*.xls表頭

2.5.2. 差異基因的熱圖聚類

??將所有比較組的差異基因取并集之后作為差異基因集。兩組以上的實(shí)驗(yàn)，可對(duì)差異基因集進(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)告展示了樣品間的聚類，具體如下圖所示。

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

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

結(jié)果文件：

差異基因的熱圖結(jié)果：Differential/heatmap/

2.5.3. 差異基因的火山圖分布

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

圖 4 差異基因火山圖

圖中橫坐標(biāo)為log2FoldChange值，縱坐標(biāo)為-log10padj或-log10pvalue，藍(lán)色的虛線表示差異基因篩選標(biāo)準(zhǔn)的閾值線

結(jié)果文件：

差異基因的火山圖結(jié)果：Differential/volcano/volcano.png

2.6. 富集分析

??我們根據(jù)基因表達(dá)量分析得到差異基因之后，必須進(jìn)一步落到基因的功能上來。對(duì)于轉(zhuǎn)錄組分析而言，往往涉及到成千上萬個(gè)基因，這會(huì)使分析變得很復(fù)雜。解決思路是將一個(gè)基因列表分成多個(gè)部分，從而減少分析的復(fù)雜度。為了解決怎么分成不同類，通常會(huì)對(duì)基因功能進(jìn)行富集分析, 期望發(fā)現(xiàn)在生物學(xué)過程中起關(guān)鍵作用的生物通路, 從而揭示和理解生物學(xué)過程的基本分子機(jī)制。功能富集分析可以將成百上千個(gè)基因、蛋白或者其他分子分到不同的通路中，以減少分析的復(fù)雜度。另外，在兩種不同實(shí)驗(yàn)條件下，激活的通路顯然比簡(jiǎn)單的基因或蛋白列表更有說服力?；蚬δ芨患治鍪紫纫獦?gòu)建基因集（gene set，如GO和KEGG數(shù)據(jù)庫等），也就是基因組注釋信息進(jìn)行分類。然后再把我們的目標(biāo)基因集（差異基因集或者其他基因集）映射到背景基因集上，注意區(qū)分注釋與富集。

??我們采用clusterProfiler軟件對(duì)差異基因集進(jìn)行GO功能富集分析，KEGG通路富集分析等。富集分析基于超幾何分布原理，其中差異基因集為差異顯著分析所得差異基因并注釋到GO或KEGG數(shù)據(jù)庫的基因集，背景基因集為所有進(jìn)行差異顯著分析的基因并注釋到GO或KEGG數(shù)據(jù)庫的基因集。富集分析結(jié)果是對(duì)每個(gè)差異比較組合的所有差異基因集、上調(diào)差異基因集、下調(diào)差異基因集進(jìn)行富集。本報(bào)告中展示的表格是選取某一個(gè)比較組合的富集分析結(jié)果，圖片是所有組合的富集分析結(jié)果。

圖 5 基因富集分析原理圖

2.6.1. 富集分析結(jié)果文件

結(jié)果文件夾：

• all 分析結(jié)果文件夾：result/Enrichment/all/

• up 分析結(jié)果文件夾：result/Enrichment/up/

• down 分析結(jié)果文件夾：result/Enrichment/down/

• all 網(wǎng)頁預(yù)覽圖：result/Enrichment/all-pdf.html

• up 網(wǎng)頁預(yù)覽圖：result/Enrichment/up-pdf.html

• down 網(wǎng)頁預(yù)覽圖：result/Enrichment/down-pdf.html

說明：

• all/up/down分別對(duì)應(yīng)總差異基因，上調(diào)差異基因，下調(diào)差異基因進(jìn)行對(duì)應(yīng)的富集分析。

表頭說明： （Results/*enrich_*/gene.ego_*.csv GO富集結(jié)果列表）

表頭說明： （Results/*enrich_*/gene.KEGG*.csv KEGG富集結(jié)果列表）

2.6.1. GO功能富集分析

??GO(Gene Ontology)是描述基因功能的綜合性數(shù)據(jù)庫，可分為生物過程（biological process）和細(xì)胞組成（cellular component）分子功能（Molecular Function）三個(gè)部分。GO功能富集以padj小于0.05作為為顯著性富集的閾值，富集結(jié)果見結(jié)果文件。

??從GO富集分析結(jié)果中，選取最顯著的30個(gè)Term繪制柱狀圖進(jìn)行展示，若不足30個(gè)，則繪制所有Term，按生物過程、細(xì)胞組分和分子功能三大類別及差異基因上下調(diào)分類畫的柱狀圖。

??有向無環(huán)圖(Directed Acyclic Graph，DAG)為差異基因GO富集分析結(jié)果的圖形化展示方式。圖中，分支代表包含關(guān)系，從上至下所定義的功能范圍越來越小，選取每個(gè)差異比較組合的GO富集結(jié)果最顯著性前5位的GO Term作為有向無環(huán)圖的主節(jié)點(diǎn)，并通過包含關(guān)系，將相關(guān)聯(lián)的GO Term一起展示，顏色的深淺代表富集程度。我們的項(xiàng)目中分別繪制生物過程、分子功能和細(xì)胞組分的DAG圖。

圖 6 GO富集分析柱狀圖

圖中縱坐標(biāo)為GO Term，橫坐標(biāo)為GO Term富集的顯著性水平，數(shù)值越高越顯著

圖 7 GO富集分析散點(diǎn)圖

圖中橫坐標(biāo)為注釋到GO Term上的差異基因數(shù)與差異基因總數(shù)的比值，縱坐標(biāo)為GO Term

圖 8 GO富集分析DAG圖

每個(gè)節(jié)點(diǎn)代表一個(gè)GO術(shù)語，方框代表的是富集程度為TOP5的GO，顏色的深淺代表富集程度，顏色越深就表示富集程度越高，每個(gè)節(jié)點(diǎn)上展示了該TERM的名稱及富集分析的padj

2.6.2. KEGG通路富集分析

??KEGG(Kyoto Encyclopedia of Genes and Genomes)是整合了基因組、化學(xué)和系統(tǒng)功能信息的綜合性數(shù)據(jù)庫。KEGG通路富集以padj小于0.05作為顯著性富集的閾值，富集結(jié)果見結(jié)果文件。

??從KEGG富集結(jié)果中，選取最顯著的20個(gè)KEGG通路繪制柱狀圖進(jìn)行展示，若不足20個(gè)，則繪制所有通路，如下圖所示。圖中橫坐標(biāo)為通路富集的顯著性水平，數(shù)值越高越顯著，縱坐標(biāo)為KEGG通路。

??從KEGG富集結(jié)果中，選取最顯著的20個(gè)KEGG通路繪制散點(diǎn)圖進(jìn)行展示，若不足20個(gè)，則繪制所有通路，如下圖所示。圖中橫坐標(biāo)為注釋到KEGG通路上的差異基因數(shù)與差異基因總數(shù)的比值，縱坐標(biāo)為KEGG通路，點(diǎn)的大小代表注釋到KEGG通路上的基因數(shù)，顏色從紅到紫代表富集的顯著性大小。

圖 9 KEGG富集分析柱狀圖

圖中橫坐標(biāo)為通路富集的顯著性水平，數(shù)值越高越顯著，縱坐標(biāo)為KEGG通路。

圖 10 KEGG富集散點(diǎn)圖

圖中橫坐標(biāo)為注釋到KEGG通路上的差異基因數(shù)與差異基因總數(shù)的比值，縱坐標(biāo)為KEGG通路