學習目標

• 瞭解 RNA-seq 和差異表達基因的分析流程
• 瞭解如何設計實驗
• 瞭解如何使用 R 語言進行資料分析

1. 簡介

在過去的十年中，RNA-seq 已成為轉錄組差異表達基因和 mRNA 可變剪下分析不可或缺的技術。正確識別哪些基因或轉錄本在特定條件下的表達情況，是理解生物反應過程的關鍵。

在本教程中，將藉助許多R包，帶你進行一個完整的 RNA-seq 分析過程。將從讀取資料開始，將偽計數轉換為計數，執行資料分析以進行質量評估並探索樣本之間的關係，執行差異表達分析，並在執行下游功能分析之前直觀地檢視結果。下面是流程圖。

2. 資料集

本教程將使用Kenny PJ et al, Cell Rep 2014

重複的情況如下圖：

使用這些資料，我們將探究 MOV10 基因表達出現不同的轉錄模式。此資料集的原始序列是從 Sequence Read Archive (SRA) 獲得的。然後經過之前介紹的工具，進行處理。所有的操作都是在 (Linux/Unix) 環境進行的。

MOV10：是一種 RNA 解旋酶，在 microRNA 通路的背景下與性別相關發育有關。

3. 問題

1. MOV10 基因的表達變化會產生什麼影響？
2. 變化之間是否有共同的特徵？

4. 配置

開啟 RStudio 併為此分析建立一個新專案。

1. 轉到 File 選單並選擇 New Project。

2. 選擇 New Directory ，然後建立 DEanalysis目錄。

3. RStudio 會自動開啟該專案。

使用 getwd()，檢查是否在正確的工作目錄中。返回的結果應該是：path/DEanalysis

（考慮到每個人的路徑不同，因此只需要最後是/DEanalysis即可）。在您的工作目錄中，建立兩個新目錄：meta

現在我們需要獲取用於分析的檔案：Mov10，點選即可下載（不能下載的，可以在文末連結獲取）。下載 zip 檔案後，您需要解壓它。將建立一個 data 目錄，其中的子目錄對應於我們資料集中的每個樣本。

接下來，我們將下載annotation file 用於將轉錄本識別符號轉換為基因名稱（如下圖）。此檔案是從 R 包 AnnotationHub 得到的（後續將介紹如何獲取過程）。

然後用 RStudio 開啟之前的 DEanalysis目錄，建立一個 de_script.R 檔案，寫入下面的註釋，並儲存。

## Gene-level differential expression analysis using DESeq2
## 使用 DESeq2 進行差異表達基因分析

完成以上步驟後，最後的工作目錄如下圖：

5. 載入包

分析將使用幾個 R 包，一些是從 CRAN 安裝的，另一些是從 Bioconductor 安裝的。要使用這些包，需要載入包。將以下內容新增到指令碼中。

library(DESeq2)
library(tidyverse)
library(RColorBrewer)
library(pheatmap)
library(DEGreport)
library(tximport)
library(ggplot2)
library(ggrepel)

6. 資料匯入

Salmon 的主要輸出是一個 quant.sf 檔案，資料集中的每個樣本都有一個這樣的檔案。該檔案如下圖所示：

• 內容如下

1. 轉錄本識別符號
2. 轉錄本長度
3. 有效長度（effective length）
4. TPM（transcripts per million），使用有效長度計算
5. 估計的計數

effective length：

The sequence composition of a transcript affects how many reads are sampled from it. While two transcripts might be of identical actual length, depending on the sequence composition we are more likely to generate fragments from one versus the other. The transcript that has a higer likelihood of being sampled, will end up with the larger effective length. The effective length is transcript length which has been “corrected” to include factors due to sequence-specific and GC biases.

將使用 tximport 包來為 DESeq2 準備 quant.sf檔案。需要做的第一件事是建立一個變數，其中包含每個 quant.sf 檔案的路徑。然後將名稱新增到我們的 quant 檔案中，這將使我們能夠輕鬆區分最終輸出矩陣中的樣本。

## 列出所有檔案
samples <- list.files(path = "./data", full.names = T, pattern="salmon$")

## 獲取檔名和路徑的向量
files <- file.path(samples, "quant.sf")

## 為每個檔案單獨命名
names(files) <- str_replace(samples, "./data/", "") %>% 
                str_replace(".salmon", "")

Salmon 索引是使用 Ensembl ID 列出的轉錄序列生成的，但 tximport 需要知道這些轉錄本來自哪些基因。我們將使用下載的 annotation file 來提取轉錄本的基因資訊。

# 載入註釋檔案
tx2gene <- read.delim("tx2gene_grch38_ens94.txt")

# 檢視
tx2gene %>% View()

tx2gene 是一個3列的 data frame，包含：轉錄本ID；基因ID；基因symbol。

tximport() 函式從各種外部軟體（例如 Salmon、Kallisto）匯入轉錄水平計數，並彙總到基因水平或輸出轉錄水平矩陣。有可選的引數來使用出現在 quant.sf 檔案中的丰度估計值或計算替代值。

對於我們的分析，需要基因水平的非標準化或“原始”計數估計來執行 DESeq2 分析。由於基因計數矩陣是預設值，因此我們只需修改一個附加引數即可指定獲取“原始”計數值。 countsFromAbundance 的選項如下：

• no（預設）：這將採用 TPM 中的值（作為我們的縮放值）和 NumReads（作為我們的“原始”計數）列，並將其摺疊到基因級別。
• scaledTPM：這是將 TPM 放大到文庫大小作為“原始”計數。
• lengthScaledTPM：這用於從 TPM 生成“原始”計數表。 “原始”計數值是通過使用 TPM 值 x featureLength x 文庫大小生成的。這些代表與原始計數在同一尺度上的數量。

TPM 計算過程：

1. reads per kilobase (RPK)：將讀取計數除以每個基因的長度（以千鹼基為單位）

2. “per million” scaling factor：計算樣本中的所有 RPK 值並將此數字除以 1,000,000。

3. 將 RPK 值除以“per million” scaling factor。

# 執行 `tximport`
txi <- tximport(files, type="salmon", tx2gene=tx2gene[,c("tx_id", "ensgene")], countsFromAbundance="lengthScaledTPM")

7. 資料檢視

txi 物件是一個簡單的列表，其中包含丰度、計數和長度的矩陣。另一個列表元素 countsFromAbundance 攜帶 tximport 中使用的字元引數。長度矩陣包含每個基因的平均轉錄本長度，可用作基因水平分析的偏移量。

attributes(txi)

$names
[1] "abundance"           "counts"              "length"        [4] "countsFromAbundance"

我們將按原樣使用 txi 物件作為 DESeq2 的輸入，但會將其儲存到下一課。現在讓我們看一下計數矩陣。你會注意到有十進位制值，所以讓我們四捨五入到最接近的整數並將其轉換為 dataframe。

# 檢視 count
txi$counts %>% View()

# 將 count 寫入到物件
data <- txi$counts %>% 
  round() %>% 
  data.frame()

歡迎Star -> 學習目錄

更多教程 -> 學習目錄

本文由mdnice多平臺釋出