下一代測序 (NGS)/RNA

本指南旨在為沒有分析下一代資料的經驗的人提供一個易於理解的 RNA-seq 資料分析指南。但是，假設對 R、下一代測序程式和使用 UNIX shell 有基本的瞭解。這裡描述的大多數步驟都在評論文章中概述。^[1]

• 它主要由 Matthew Young (myoung@wehi.edu.au) 編寫，目前仍在開發中。
• 病原體示例由 B. Usadel 提供，並使用了一組不同的工具。

RNA

轉錄組學簡介。

對於執行 RNA-seq 的每個樣本，您通常會收到一個包含數百萬個短 (25-300bp) DNA 序列（稱為讀）和指示每個鹼基呼叫置信度的質量分數的檔案。但是，有一些重要的常見變化，這些變化取決於使用的平臺和協議。這些包括但不限於

• 鹼基空間或顏色空間
• 雙端/配對或單端/非配對
• 鏈特異性或非鏈特異性

以下部分將詳細介紹每種型別。

鹼基空間與顏色空間

用於第二代 RNA 測序的兩個主要平臺由 Illumina 和 ABI Solid 生產。雖然兩者都產生數百萬個短讀，但它們是以略微不同的方式測序和報告的。Solid 平臺使用一種測序技術，該技術透過一次連線兩個鹼基對來生成讀資訊。讀中的每個鹼基對都透過兩個（可能不同的）二核苷酸測序兩次。為了利用這種測序化學，Solid 報告其讀不是作為核苷酸序列，而是作為四種顏色的序列，其中每種顏色代表鹼基之間的轉換，稱為“顏色空間編碼”。^[2] 另一方面，Illumina 平臺一次讀取片段上的一個鹼基，直到達到所需的讀長度。然後，觀察到的鹼基被報告為輸出。

結果是，分析 RNA-seq 資料的工具取決於用於生成短讀資料的平臺。雖然可以將顏色空間讀轉換為鹼基空間（反之亦然），但這樣做會導致資料中出現嚴重的偏差，應不惜一切代價避免。因此，讀應保留其本機格式，應使用適當的工具分析它們。

雙端/配對讀

標準 RNA-seq 協議涉及逆轉錄 mRNA (cDNA) 的隨機剪下，然後從片段的一端測序一個短“讀”。這意味著只有片段的前 25-300bp 是已知的（取決於讀的長度），而其餘片段仍然未測序。事實上，由於片段化是隨機的，每個片段的長度也是未知的，儘管通常會應用大小選擇步驟。

雖然化學物質的精確度不足以對整個片段進行測序，但可以應用一個巧妙的技巧，即從片段的兩端獲取一個短讀，從而產生一對短讀，一個來自片段的一端，一個來自另一端。執行此操作的讀稱為雙端讀或配對讀。雙端讀允許推斷有關間隔序列的更多資訊，並且對於從頭轉錄組構建和檢測結構變異（如插入缺失）特別有用。

鏈特異性讀

RNA-seq 資料可以分為鏈特異性和非鏈特異性兩種。如果資料是非鏈特異性的，則無法直接從序列中識別片段轉錄的鏈。此外，由於 RNA-seq 協議通常涉及形成雙鏈 cDNA 以便於測序，因此返回的序列與來源 DNA 序列的反向互補序列的可能性與原始 DNA 序列的可能性相同。在實際應用中，這意味著雖然有一半的 reads 對映到正向鏈，一半對映到反向鏈，但這種對映不包含任何關於 RNA 從哪條鏈轉錄的資訊。

另一方面，鏈特異性 RNA-seq 資料保留了鏈資訊，從而可以識別 RNA 從哪條鏈轉錄。

典型工作流程

TopHat

Tophat 是一款工具，可以作為 RNA-Seq reads 的快速拼接連線對映器。Tophat 使用 Bowtie 將 RNA-Seq reads 對映到哺乳動物基因組，Bowtie 是一款超高通量短讀比對器，利用 Burrows-Wheeler 索引。在 Bowtie 比對 reads 後，Tophat 會分析對映結果，以識別外顯子之間的拼接連線。Tophat 是為數不多的幾種能夠有效處理內含子區域的工具之一，內含子區域在長度大於 50 bp 的 reads 中更為普遍。我們可以使用 HISAT 代替 tophat

參考資料集

人類

基因組參考聯盟提供了一個人類參考基因組，該基因組由來自不同人群的幾個個體構建而成。並且，透過修正現有的組裝錯誤（包括但不限於在所使用個體之間錯誤混合單倍型結構）不斷提高參考的質量。

CCDS

一致性 CDS (CCDS) 專案是一項合作努力，旨在識別一組核心的人類和小鼠蛋白質編碼區域，這些區域是一致註釋的且質量很高。註釋是 NCBI 的 RefSeq、EBI 的 Ensembl 和桑格研究所的 Havana 之間的共識。

GENCODE

GENCODE 是 ENCODE 擴充套件專案的一個子專案，旨在註釋整個人類基因組中所有基於證據的基因特徵。

HAVANA

桑格研究所的 HAVANA 小組手動註釋人類基因組，為基因位點的全部複雜性和可能不被自動化註釋系統很好地滿足的特徵提供全面的註釋。

RefSeq

RefSeq 是 NCBI 提供的一種參考序列註釋服務。

UCSC 基因

UCSC 已知基因是透過完全自動化的過程構建的，該過程基於來自 Swiss-Prot/TrEMBL (UniProt) 的蛋白質資料以及來自 Genbank 的相關 mRNA 資料。

SEQwiki 內容轉儲

RNASeq 資料可以做什麼

為了分析你的 RNASeq 資料，不同的目標需要不同的預處理步驟。

尋找差異表達基因

如果你想尋找差異表達基因，許多 BioConductor 軟體包要求你不要將你的讀取計數轉換為例如 rpkm，而要使用原始讀取計數。原因是統計模型使用原始讀取計數資料。（讀取次數越少，短噪聲越多，而較低的 rpkm 值仍可能與許多讀取相關聯，例如在非常長的基因情況下）。請注意，即使你使用 RNASeq 資料，你仍然需要生物學重複。

此外，你將希望使用為 RNAseq（即計數資料）設計的軟體包，例如DESeq，bayseq，NBPSeq 或 edgeR，它們都使用負二項式模型。

一個典型的流程可能如下所示

1. 對你的資料進行質量控制
2. 讀取比對（參見說明）
3. 計算每個特徵的讀取次數
4. 統計分析

讀取計數/彙總

在將讀取比對到你的基因組後，你需要彙總這些讀取。這可能是一個關鍵步驟，因為有許多不同的彙總方法。一種方法可以是十分嚴格，只統計完全明確的讀取，但是在下游分析中需要考慮到這一點。可以使用例如 HTSeq count 或 Bioconductor Iranges 模組來彙總資料

統計分析

在測試差異表達時，你應該使用一個可以對計數資料進行建模的軟體包。在 Bioconductor 中，可以選擇 DESeq, edgeR 和 bayseq。

差異外顯子使用

除了基因表達的總體變化外，還可能存在等位基因或特定外顯子的差異丰度。可以使用 DEXSeq 來測試差異丰度外顯子，另請參見 Genome Research 論文。

類別富集（例如 GO 富集）

與所有富集分析一樣，你發現的類別可能會受到檢測能力的影響。使用 RNASeq，檢測差異表達的能力（在相同型別的錯誤下，例如相同的調整後的 p 值）對於計數更多的基因來說更高。例如，如果所有光合作用基因都高度表達，你更有可能發現其中任何一個基因發生差異表達，而不是其他基因。現在，如果你發現光合作用基因在你的差異表達基因列表中富集，這可能僅僅是因為它們更容易被發現。

Bioconductor 中的Goseq 軟體包的作者指出，更長的基因往往會產生更多的計數，並且他們提供了進行 GO 類別富集分析的軟體，該軟體可以根據基因長度或直接根據表達強度進行調整。

估計表達高度

一個典型的流程可能如下所示

1. 讀取比對
2. 計算每個特徵的讀取次數
3. 根據基因長度進行歸一化

質量控制 (QC)

有很多方法可以對測序資料進行質量檢查。一些流行的選擇是：FastQC， DNAA 或者你可以在 R 中進行操作或使用 HTSeq

儘管可以使用更復雜的質量指標，但一個基本的檢查是，序列組成在讀取的長度上沒有太大變化，並且質量得分沒有過低，這是一個良好的起點。這些檢查（以及其他一些檢查）可以透過將 fastq 檔案載入到 fastQC 程式中來執行。

比對讀取

除非你的目標是進行從頭轉錄組組裝，否則任何分析的第一步都是將你的數百萬條短讀取比對到某種參考。這通常是提取 RNA 的物種的參考基因組。有許多不同的比對器可以執行短讀取比對到參考（列表可在 http://en.wikipedia.org/wiki/List_of_sequence_alignment_software 和 http://seqanswers.com/wiki/Special:BrowseData/Bioinformatics%20application?Bioinformatics_method=Mapping 中找到）。每一個都有特定的優點和缺點，通常涉及速度和敏感性之間的權衡。在本指南中，我們將使用比對器 BOWTIE，因為它正在積極開發、速度快、使用廣泛，並支援基本空間和顏色空間讀取。

作為起點，我假設你擁有以下檔案

對於單端讀取

sample.fa 或 sample.cfa（cfa 用於顏色空間，fa 用於基本空間）

對於雙端讀取

sample_pair1.fa 和 sample_pair2.fa 或 sample_pair1.cfa 和 sample_pair2.cfa（cfa 用於顏色空間，fa 用於基本空間）

獲取參考

為了將短序列比對到參考基因組，該基因組必須被轉換為BOWTIE可以使用的索引。一些常見的基因組的預構建索引可以從BOWTIE網站下載。如果您的基因組不可用，您將需要從包含參考基因組的fasta檔案中自己構建索引，這些檔案可以從UCSC 透過FTP獲取。可以使用以下命令完成此操作：

   bowtie-build reference.fa reference_name

引數“reference_name”是用於從現在開始引用此參考基因組的唯一識別符號。如果您的序列為顏色空間，則透過新增-C命令構建顏色空間索引：

   bowtie-build -C reference.fa reference_name

此命令將輸出6個檔案，分別命名為reference_name.1.ebwt、reference_name.2.ebwt、reference_name.3.ebwt、reference_name.4.ebwt、reference_name.rev.1.ebwt和reference_name.rev.2.ebwt。

無論您如何獲取索引，為了使BOWTIE能夠使用它，上述六個檔案都需要放在BOWTIE索引目錄中。如果您不確定您的索引目錄是什麼，它由環境變數BOWTIE_INDEXES指向，因此：

   echo $BOWTIE_INDEXES

將顯示您應該放置索引檔案的路徑。

如果您希望使用除了BOWTIE之外的其他比對工具，您也需要從參考基因組構建索引。有關更多資訊，請參閱您首選比對工具的文件。

將序列比對到參考基因組

將參考基因組構建成BOWTIE可以使用的索引後，我們現在想將我們的資料比對到這個索引。BOWTIE提供了豐富的命令列選項，可用於調整比對演算法及其處理輸入/輸出的方式。這些命令列選項在BOWTIE手冊中進行了詳細描述，但有一些標誌在RNA-seq分析中經常使用，值得在這裡提及。

--sam或-S告訴BOWTIE以SAM格式而不是BOWTIE格式輸出比對結果。SAM正迅速成為報告短序列比對的標準，並得到許多下游分析工具的支援。除非有充分的理由，否則應始終指定此標誌。

--best標誌告訴BOWTIE保證它報告的比對在所有找到的匹配中與參考基因組的錯配數量最少。它還執行一些其他理想的功能，例如消除鏈偏向性。 這些好處的代價是--best稍微慢一些，但在幾乎所有情況下，這種差異都可以忽略不計。請注意，--best不適用於雙端序列。除非速度是一個主要考慮因素，否則應啟用此標誌。

每個鹼基都有一個與之相關的質量得分，該質量得分以PHRED量表 (http://en.wikipedia.org/wiki/Phred_quality_score) 報告，其中較低的得分意味著對鹼基呼叫準確性的信心較低。對於每個候選比對，BOWTIE會將與參考基因組不匹配的鹼基的質量得分加起來。任何錯配質量得分總和大於-e的匹配位置都被認為是無效的，不會被報告。-e的預設值為70，是在序列長度較短（~30bp）時最佳化的，但不適用於當今常用的較長的序列長度。此外，任何來自參考基因組的真正生物學變異（例如SNP）理論上都應該具有很高的質量得分。因此，具有SNP的序列在錯配質量得分總和指標上的得分將非常低。出於所有這些原因，建議使用者將-e提高到預設值以上，除非已知參考基因組是RNA生物學來源的極佳代表，並且序列中的錯誤數量很少。

-p標誌設定要使用的同時執行緒數。由於短序列比對實際上是無限可並行的，因此將其設定為可用的CPU核心數。

考慮到這些因素，我們使用以下命令將短序列資料對映到參考基因組：

   bowtie -p 8 --sam --best reference_name sample.fa aligned_reads.sam

如果使用顏色空間，則需要新增-C標誌。

   bowtie -p 8 --sam -C --best reference_name sample.cfa aligned_reads.sam

如果您的資料是雙端的，則需要使用-1和-2標誌指定兩個配對檔案。

   bowtie -p 8 --sam reference_name -1 sample_pair1.fa -2 sample_pair2.fa aligned_reads.sam

更復雜的比對

正在測序的cDNA片段不是來自基因組，而是來自轉錄組。轉錄組是透過組合來自基因組的外顯子形成的，這就是為什麼將序列對映到基因組是一個很好的近似。但是，這樣做會丟失任何跨越外顯子-外顯子邊界的序列的對映能力。序列越長，這個問題就越嚴重，因為序列更有可能跨越邊界而無法對映。根據所需的後續分析，這種在外顯子邊界處的覆蓋丟失可能是或可能不是一個問題。

外顯子連線庫

可以從已知的註釋外顯子構建外顯子連線庫，然後嘗試將那些無法對映到該序列的序列對映到庫。這種方法只能捕獲跨越已知註釋的序列，限制了其用途，並使結果偏向於註釋良好的基因/基因組。為此，需要構建一個新的參考基因組，其中包含參考基因組以及新的外顯子-外顯子連線序列。即使您只將無法比對到基因組的序列對映到庫，也要確保原始基因組參考仍然包含在參考基因組中，以便在確定比對時，序列可以在所有可能的對映位置之間競爭。您還需要決定從外顯子-外顯子邊界兩側取多少序列。從每個外顯子取的序列量必須小於序列長度，否則任何靠近外顯子邊界但沒有超過邊界的序列都將對映到兩個位置，即外顯子連線庫和基因組本身。這意味著，如果您打算修剪您的序列，那麼從外顯子連線兩側取的序列量必須小於修剪後的序列的長度。最後，您需要決定要考慮哪些外顯子-外顯子連線組合。您只考慮基因內、染色體內的剪接，而不考慮外顯子重新定位？所有這些情況在不同的樣本中以不同的頻率發生，您考慮的可能性越多，計算複雜度就越大。

要構建連線庫，您首先需要一個基因組註釋。這些可以從UCSC表格瀏覽器下載。然後，您可以編寫自己的工具來建立一個包含所有外顯子連線的fasta檔案，或者使用現有的工具，例如USeq軟體包中包含的“Make Splice Junction Fasta”應用程式。 一旦您擁有包含外顯子-外顯子連線庫的fasta檔案，您就需要將其與參考基因組的fasta檔案合併。

   cat reference.fa junctions.fa >reference_and_junctions.fa

然後，如上所述構建新的bowtie索引：

   bowtie-build reference_and_junctions.fa junction_lib

最後，您以與參考基因組相同的方式將序列對映到庫。例如：

   bowtie -p 8 --sam --best junction_lib sample.fa aligned_reads.sam

其他選項

如果這種方法仍然無法對映合理的序列數量，那麼還有其他替代方法可以探索，即使計算成本更高。例如，您可以嘗試使用“從頭”剪接連線查詢器從資料本身估計剪接連線。有很多工具嘗試這樣做，一些例子包括TopHat, SplitSeek, PerM, 和SOAPsplice（以前稱為SOAPals）。 轉錄組的完全從頭組裝也是一種可能性，儘管為此工作良好需要非常高的覆蓋率。雙端資料對此任務非常有益，因為片段的任一端對映到不同的外顯子，是外顯子連線的非常強烈的證據。在綜述文章中，您可以找到此類工具的列表（見表格1）。

最近發表在基因組生物學上的論文進一步提出了一種自動框架，可以根據資料集的關鍵特徵（例如參考基因組、SNP 率、讀長等）輕鬆地對比和選擇比對器。Teaser 可作為線上工具使用，或者您可以下載並自定義自己的版本。

差異表達

表達分析的一個常見用途是在兩個實驗條件（例如野生型與敲除）之間尋找基因或其他感興趣物件的表達水平差異。為此，我們需要將資料從基因組座標的讀序列列表轉換為計數表。我們在這裡採用的策略是使用 Rsamtools 包將短讀序列載入到 R 中，然後計算與某些註釋物件重疊的讀序列數量，該物件通常類似於從 UCSC 下載的基因集合。轉換後，可以執行測試以尋找表達水平的統計學顯著差異。

讀序列的彙總

使用 SAMtools 壓縮比對的讀序列

為了能夠將數百萬個比對的讀序列載入到 R 中的記憶體中，我們需要建立人類可讀的 SAM 輸出的二進位制壓縮版本。該檔案將包含所有相同的資訊，但記憶體佔用更小，並且可以快速搜尋。要建立此類檔案，我們首先需要安裝 SAMtools。 接下來，我們需要使用 fasta 檔案構建參考索引。執行以下操作：

   samtools faidx reference.fa

這將建立一個檔案 reference.fa.fai。接下來，我們將 SAM 轉換為 BAM。BAM 檔案包含與 SAM 檔案相同的所有資訊，但經過壓縮，使其更加節省空間且可搜尋。

   samtools import reference.fa.fai aligned_reads.sam aligned_reads.bam

最後，我們需要建立讀序列索引，以便可以快速搜尋它們。為此，我們首先需要對 BAM 檔案進行排序。

   samtools sort aligned_reads.bam aligned_reads_sorted

這將建立 aligned_reads_sorted.bam，我們現在對其進行索引。

   samtools index aligned_reads_sorted.bam

這將建立索引檔案 aligned_reads_sorted.bam.bai。

如果您沒有參考的原始 fasta 檔案，因為您從 BOWTIE 網站下載了預構建的索引（或者因為您在製作自己的索引後丟失了 fasta 檔案），您可以透過執行以下命令重建源 fasta 檔案。

   bowtie-inspect reference_name>reference.fa

在 R 中使用 BAM 檔案

我們的目標是使用 R 根據每個樣本的基因對讀序列進行彙總。為此，我們將使用我們新建立的短讀序列的壓縮表示（已排序、已索引的 BAM 檔案）。

獲取基因資訊

我們要做的第一件事是定義基因在染色體座標中的位置。為此，我們使用 GenomicFeatures 包。 該包允許我們使用以下命令從 UCSC 基因組瀏覽器下載基因資訊

   library(GenomicFeatures)
   txdb=makeTranscriptDbFromUCSC(genome='hg19',tablename='ensGene')

可以使用各種基因組和基因 ID，但作為示例，我們將使用最新的 human 基因組和 ENSEMBL 基因 ID。變數“txdb”現在包含了我們所需的所有資訊，但為了使用它，我們需要進行一些處理。

   tx_by_gene=transcriptsBy(txdb,'gene')

這將生成一個 GRangesList 物件，它是一個 GRanges 物件列表，其中每個 GRanges 物件都是一個基因，條目是其轉錄本的基因組座標。我們將使用 countOverlaps 函式將此物件與包含短讀序列的物件重疊，從而計算出哪些讀序列與哪些基因重疊。

雖然我們選擇透過包含落入基因內的所有讀序列來進行彙總，但該過程對於任何 GRanges 或 GRangesList 物件都將相同。例如，如果您希望只包含與外顯子重疊的讀序列，您可以建立一個不同的 GRangesList 物件。

   ex_by_gene=exonsBy(txdb,'gene')

在 R 中彙總讀序列

首先，我們將 aligned_reads_sorted.bam 檔案載入到 R 中。

   library(Rsamtools)
   reads=readBamGappedAlignments("aligned_reads_sorted.bam")

檢查註釋和讀序列的相容性

在我們嘗試將讀序列與註釋進行比較之前，我們需要進行一些檢查，以確保一切正常。如果您的 RNA-seq 資料沒有附帶鏈資訊，那麼我們無法知道讀序列是從哪個鏈轉錄的。但是，比對過程會將其比對到一個鏈或另一個鏈（理論上，兩者可能性相同），因此讀序列將被分配一個人工鏈。當我們計算與基因（或其他特徵）重疊的讀序列數量時，只有比對到基因所在的鏈的讀序列才會被計數，大約一半的讀序列會丟失。為了避免這種情況，我們需要將讀序列的鏈值設定為“*”（未知）。

此外，比對軟體使用的染色體名稱（最終由參考基因組 fasta 檔案中的染色體名稱決定）通常與註釋中給出的染色體名稱不同。為了使比較函式正常工作，這些名稱需要轉換為相同的命名約定。透過列出讀序列和註釋的全部染色體，可以輕鬆檢查名稱是否匹配。請記住，我們的註釋資料儲存在“tx_by_gene”中，我們的短讀序列儲存在“reads”中。

   #The annotations have chromosomes called
   names(seqlengths(tx_by_gene))
   #The reads have chromosomes called
   as.character(unique(rname(reads)))

如果染色體名稱相同（或對於您關心的染色體名稱相同），則不需要進行名稱轉換。另一方面，如果它們不同，我們需要更改讀序列或註釋的命名約定。事實證明，更改讀序列的名稱通常更容易，但無論哪種情況，過程都是相同的。以下示例可以很好地說明這一點。假設您有以下情況

   #The annotations have chromosomes called
   > names(seqlengths(tx_by_gene))
   [1] "chr1" "chr10" "chr11" "chr12" "chr13"
   [6] "chr13_random" "chr14" "chr15" "chr16" "chr17"
   [11] "chr17_random" "chr18" "chr19" "chr1_random" "chr2"
   [16] "chr3" "chr3_random" "chr4" "chr4_random" "chr5"
   [21] "chr5_random" "chr6" "chr7" "chr7_random" "chr8"
   [26] "chr8_random" "chr9" "chr9_random" "chrM" "chrUn_random"
   [31] "chrX" "chrX_random" "chrY" "chrY_random"

   #The reads have chromosomes called
   >as.character(unique(rname(reads)))
   [1] "10.1-129993255" "11.1-121843856" "1.1-197195432" "12.1-121257530"
   [5] "13.1-120284312" "14.1-125194864" "15.1-103494974" "16.1-98319150"
   [9] "17.1-95272651" "18.1-90772031" "19.1-61342430" "2.1-181748087"
   [13] "3.1-159599783" "4.1-155630120" "5.1-152537259" "6.1-149517037"
   [17] "7.1-152524553" "8.1-131738871" "9.1-124076172" "MT.1-16299"
   [21] "X.1-166650296" "Y.1-15902555"

因此，我們需要將讀序列染色體名稱“NO.1-Length”轉換為註釋名稱“chrNO”。

   new_read_chr_names=gsub("(.*)[T]*\\..*","chr\\1",rname(reads))

如果您不熟悉正則表示式，請參閱 R 中 gsub 的幫助檔案。new_read_chr_names 現在將包含讀序列染色體名稱，已轉換為與註釋物件 tx_by_genes 相同的格式。

現在，我們可以透過從每個讀序列物件構建 GenomicRanges 物件來同時修復染色體名稱和鏈問題。如果您有無鏈 RNA-seq 資料，並且需要轉換染色體，我們會執行以下命令

   reads=GRanges(seqnames=new_read_chr_names,ranges=IRanges(start=start(reads),end=end(reads)), 
     strand=rep("*",length(reads)))

如果我們只是想轉換染色體名稱。

       reads=GRanges(seqnames=new_read_chr_names,ranges=IRanges(start=start(reads),end=end(reads)), 
     strand=strand(reads))

如果我們只是想使每個讀序列在鏈方面都模稜兩可。

   reads=GRanges(seqnames=rname(reads),ranges=IRanges(start=start(reads),end=end(reads)),
   strand=rep("*",length(reads)))

請注意，如果您將讀序列比對到外顯子連線庫，則每個外顯子連線將具有自己的“染色體”。如果您希望這些連線讀序列包含在彙總中，則必須將它們中的每一個都轉換為基因組座標。由於每個讀序列來自兩個不同的基因組位置（外顯子連線的任一側），因此您必須決定如何為每個讀序列分配基因組座標。

計算讀序列數量

最後，我們獲得了與每個基因（或您感興趣的其他任何內容）重疊的讀序列數量。

   counts=countOverlaps(tx_by_gene,reads)

“counts”現在將包含一個數值向量，其中第 i 個條目是與“tx_by_gene”中第 i 個基因重疊的讀序列數量。

差異表達測試

由於我們的目標是比較條件，因此我們將擁有多個讀序列泳道，每個條件可能有多個泳道。使用上一節中概述的程式，我們可以計算與每個實驗條件下每個重複的感興趣特徵（例如基因）重疊的讀序列數量。接下來，我們將它們組合成一個計數表。

   toc=data.frame(condition1_rep1=counts1.1,condition1_rep2=counts1.2,
     condition2_rep1=counts2.1,condition2_rep2=counts2.2,stringsAsFactors=FALSE)
   rownames(toc)=names(tx_by_gene)

等等，只要有條件和重複，就重複多少次。這裡的約定是 countsn.m 是一個向量，包含條件 n 中重複 m 中與 tx_by_gene 給出的基因重疊的讀取次數。

歸一化

已經證明，少數高表達基因會消耗大量的總序列。由於這種情況會在泳道和實驗條件之間發生變化，以及文庫大小，因此在測試差異表達時，需要進行某種樣本間歸一化。歸一化的選擇與用於確定任何基因在條件之間是否存在顯著差異表達 (DE) 的檢驗無關。例如，分位數歸一化會產生非整數計數，使得基於計數資料的檢驗（例如廣泛使用的泊松模型或負二項式模型）變得不適用。我們選擇使用比例因子歸一化方法，因為它保留了資料的計數性質，並且已被證明是提高 DE 檢測的有效方法。

為了進行歸一化和差異表達檢驗，我們將使用 R 包 edgeR， 儘管還有其他選擇。我們現在可以使用 TMM 方法計算歸一化因子。

   library(edgeR)
   norm_factors=calcNormFactors(as.matrix(toc))

統計檢驗

接下來，我們需要建立一個 edgeR 使用的 DGE 物件。使用 TMM 方法計算的比例因子透過按歸一化因子加權文庫大小（然後在統計模型中用作偏移量）來合併到統計檢驗中。

   DGE=DGEList(toc,lib.size=norm_factors*colSums(toc),group=rep(c("Condition1","Condition2"),c(2,2)))

組變數標識內容表中的哪些列來自哪個實驗條件或“組”。為了進行顯著性統計檢驗，我們首先估計公共離散引數

   disp=estimateCommonDisp(DGE)

最後，我們計算基因是 DE 的 p 值

   tested=exactTest(disp)

基因集檢驗 (GO)

為了使用 RNA-seq 資料準確地檢驗 DE 基因中基因集的過度表達，我們需要使用一種方法來考慮這種技術的特定偏差。GOseq 包 是一種用於在執行 GO（和其他基於基因集的檢驗）分析時解釋某些 RNA-seq 特定偏差的方法。

首先，我們必須格式化 edgeR 的輸出以供 goseq 讀取。我們將任何 Benjamini-Hochberg FDR 小於 .05 的基因稱為 DE。

   library(goseq)
   genes = as.integer(p.adjust(tested$table$p.value, method = "BH") < 0.05)
   names(genes) = row.names(tested$table)

接下來，我們計算機率加權函式，校正長度偏差，這是所有形式的 RNA-seq 資料中存在的技術偏差。

   pwf=nullp(genes,'hg19','ensGene')

如果我們想校正總讀取計數偏差，我們會使用每個基因的計數數量來計算 pwf，如下所示

   pwf=nullp(genes,bias.data=rowsum(counts[match(names(genes),rownames(counts))]))

最後，我們計算每個 GO 類別在 DE 基因中過度表達的 p 值。

   GO.pvals=goseq(pwf,'hg19','ensGene')

示例 1：差異表達和 GO 術語分析：李前列腺癌資料集

本節提供了一個易於理解的示例，用於說明上述分析流程。

描述

該資料集比較了經過和未經過類固醇激素雄激素處理的前列腺癌 LNcap 細胞系。 測序使用 Illumina GA I 完成，產生 36 bp 的單端非鏈讀取。機器的輸出是 7 個檔案（每個來自不同的測序泳道）

untreated1.fa

untreated2.fa

untreated3.fa

untreated4.fa

treated1.fa

treated2.fa

treated3.fa

請注意，此資料集有點不尋常，因為讀取中的質量得分缺失。因此，我們在進行分析時需要牢記這一點。

質量控制

這是已釋出的資料，所以這是一個相當不錯的質量控制，人們希望……

序列比對

構建參考

流程中的第一步是將所有讀取比對到參考序列。由於該資料取自人類 LNcap 細胞，因此人類基因組的最新版本是一個顯而易見的選擇。我們已使用所有標準選項安裝了 BOWTIE（版本 0.10.0）。從 BOWTIE 網站可以獲得人類索引的預構建副本，但是，為了說明從頭開始構建基因組，我們改為從 UCSC 下載基因組的 .fa 檔案。我們建立一個包含 7 個 RNA-seq 資料檔案和從 加州大學聖克魯茲分校 下載的 chromFA.tar.gz 檔案的工作目錄。檔案 chromFA.tar.gz 包含所有人類染色體的序列，包括未分配的重疊群。為了製作 BOWTIE 索引，我們需要將它們連線到一個檔案中。我們將從我們的 fasta 檔案中排除所有重疊群。

   tar -zxvf chromFA.tar.gz

我們只想要 chr1-22.fa、chrX.fa、chrY.fa 和 chrM.fa，所以刪除其他所有內容

   rm chr*_*.fa

現在，我們將所需的檔案連線起來（最好將參考序列分別儲存在每個染色體一個檔案和每個基因組一個檔案格式中，儘管 BOWTIE 索引可以從任何一種格式中生成）。

   cat chr*.fa>hg19.fa

我們構建 BOWTIE 索引並將其命名為 hg19。

   bowtie-build hg19.fa hg19

並將 BOWTIE 索引檔案移動到 BOWTIE 可以找到它們的適當位置。

   mv *.ebwt $BOWTIE_INDEXES

比對讀取

構建完人類基因組的 BOWTIE 索引後，我們現在繼續將每個泳道的讀取對映到索引。我們想使用 --best 和 --sam 標記。此時，我們回想起我們的資料缺少讀取的質量資訊。因此，我們使用 -v 3 選項，該選項在將讀取比對到基因組時忽略質量得分。

   for src_fastafile in *treated*.fa
   do
     bowtie -v 3 -p 8 --best --sam hg19 ${src_fastafile} ${src_fastafile%%.fa}.sam
   done

for 迴圈僅用於方便，等同於編寫

   bowtie -v 3 -p 8 --best --sam hg19 untreated1.fa untreated1.sam
   bowtie -v 3 -p 8 --best --sam hg19 untreated2.fa untreated2.sam
   bowtie -v 3 -p 8 --best --sam hg19 untreated3.fa untreated3.sam
   bowtie -v 3 -p 8 --best --sam hg19 untreated4.fa untreated4.sam
   bowtie -v 3 -p 8 --best --sam hg19 treated1.fa treated1.sam
   bowtie -v 3 -p 8 --best --sam hg19 treated2.fa treated2.sam
   bowtie -v 3 -p 8 --best --sam hg19 treated3.fa treated3.sam

BOWTIE 將輸出 7 個 SAM 檔案，其中包含比對的讀取。為了統計比對讀取的比例，我們在 shell 中執行以下命令（此資訊也由 BOWTIE 直接報告，但能夠自己計算它很有用）。

   awk '$3!="*"' untreated1.fa|wc -l
   wc -l untreated1.fa

第一個命令列印 BOWTIE 輸出檔案中對映的讀取數量，第二個命令輸出輸入檔案中讀取數量的 4 倍（因為 fasta 格式是每個讀取 4 行）。

我們將忽略那些跨越外顯子-外顯子邊界的讀取，並繼續進行分析。

讀序列的彙總

轉換為 BAM

為了將比對後的reads彙總成基因，我們首先需要將BOWTIE輸出的SAM檔案轉換為壓縮的BAM索引格式。首先，我們需要為人類基因組建立samtools格式的索引。為此，我們需要參考序列的fasta檔案（應該與用於建立比對reads索引的fasta檔案相同），我們已經有了，但我們會從BOWTIE索引重建它。

   bowtie-inspect hg19>hg19.fa

現在我們構建samtools索引。以下命令會生成一個.fai檔案，我們可以用它將SAM檔案轉換為BAM檔案。

   samtools faidx hg19.fa

接下來，我們將所有7個SAM檔案轉換為BAM檔案。

   for SAM in *.sam
   do
     #Convert to BAM format
     samtools import hg19.fa.fai ${SAM} ${SAM%%.sam}.bam
     #Sort everything
     samtools sort ${SAM%%.sam}.bam ${SAM%%.sam}_sorted
     #Create an index for fast searching
     samtools index ${SAM%%.sam}_sorted.bam
     #Delete temporary files
     rm ${SAM%%.sam}.bam
   done

同樣，bash for迴圈只是為了方便，上面等同於對所有7個檔案執行以下命令。

   samtools import hg19.fa.fai untreated1.sam untreated1.bam
   samtools sort untreated1.bam untreated1_sorted
   samtools index untreated1_sorted.bam
   rm untreated1.bam

在R中進行處理

接下來我們需要將排序後的、索引的BAM檔案載入到R中。由於我們使用的是非鏈特異性RNA-seq，我們需要將每個read的鏈指示器設定為模糊（透過將其設定為"*"）。在啟動R後，我們執行：

   library(Rsamtools)
   #Create a list of bam file object containing the short reads
   bamlist=list()
   src_files=list.files(pattern="*_sorted.bam$")
   for(filename in src_files){
     #Since we do not know which strand the reads were originally transcribed,
     #so set the strand to be ambiguous
     tmp=readBamGappedAlignments(filename)
     bamlist[[length(bamlist)+1]]=GRanges(seqnames=rname(tmp),
       ranges=IRanges(start=start(tmp),end=end(tmp)),
       strand=rep("*",length(tmp)))
   }
   names(bamlist)=src_files

將檔案載入到R中後，我們需要建立一個註釋物件。由於我們使用的是hg19，我們可以使用GenomicFeatures包從UCSC輕鬆下載一個。我們選擇使用ENSEMBL基因註釋。

   library(GenomicFeatures)
   txdb=makeTranscriptDbFromUCSC(genome="hg19",tablename="ensGene")

我們想比較基因的差異表達，因此我們將按基因進行彙總，並選擇將所有落在基因體內（包括內含子）的reads計入基因的計數。

   tx_by_gene=transcriptsBy(txdb,"gene")

最後，我們計算每個lane中落在每個基因的reads數量，並將結果記錄在一個計數表中。

   #Initialize table of counts
   toc=data.frame(rep(NA,length(tx_by_gene)))
   for(i in 1:length(bamlist)){
     toc[,i]=countOverlaps(tx_by_gene,bamlist[[i]])
   }
   #Fix up labels
   rownames(toc)=names(tx_by_gene)
   colnames(toc)=names(bamlist)

差異表達測試

最終獲得一個計數表後，我們現在想比較處理組和未處理組，並尋找每個基因的計數數量是否存在任何統計學上的顯著差異。我們將使用edgeR使用的負二項式模型來做到這一點。

歸一化

我們使用TMM方法（以第一個lane為參考）計算合適的縮放因子進行歸一化。

   library(edgeR)
   norm_factors=calcNormFactors(as.matrix(toc))

計數本身並沒有改變，而是將這些縮放因子用作負二項式模型中的偏移量。這被包含在edgeR需要的DGE列表物件中。

   DGE=DGEList(toc,lib.size=norm_factors*colSums(toc),group=gsub("[0-9].*","",colnames(toc)))

統計測試

接下來我們計算一個共同的離散引數，它表示資料中額外的泊松可變性。

   disp=estimateCommonDisp(DGE)

這使我們能夠計算基因差異表達的p值。

   tested=exactTest(disp)

基因本體測試

為了測試DE基因中過度表達的GO類別，我們首先需要選擇一個閾值，用於在應用多重假設校正後將基因稱為差異表達。我們選擇流行的顯著性閾值為0.05。

   library(goseq)
   #Apply benjamini hochberg correction for multiple testing
   #choose genes with a p-value less than .05
   genes=as.integer(p.adjust(tested$table$p.value,method="BH") <.05)
   names(genes)=row.names(tested$table)

現在我們計算機率加權函式，它量化了長度偏差效應。

   pwf=nullp(genes,"hg19","ensGene")

最後，我們計算每個GO類別在DE基因中過度表達的p值。

   GO.pvals=goseq(pwf,"hg19","ensGene")

示例2 差異表達：擬南芥病原體資料

 Donwnload the six fastq.gz files from here: http://www.ebi.ac.uk/ena/data/view/SRP004047&display=html

這些檔案來自一項擬南芥研究，使用受感染植物和模擬感染植物的三個重複樣本。這是NBPSeq R包的基礎資料集。

質量控制

 Download FastQC and open the files in FastQC one by one.
 You can open the files by using File->Open. 

截至版本0.94，如果您使用的是Windows，請使用Linux版本並雙擊run_fastqc.bat。

 You will see a very wiggly line for the first library. If you just look at the peaks and note the sequence you will see the pattern
 AAGAGCTCGTATGC starting at the green plateau towards the right. This is an illumina adapter sequence, which you will also see in the
 overrepresented counts tab.

檔案：FastQC.png

讀取比對

從這裡下載擬南芥基因組版本。

 ftp://ftp.arabidopsis.org/home/tair/Sequences/blast_datasets/other_datasets/CURRENT/arabidopsis.seq.gz
 It doesn't include mitochondrial or chloroplast sequences but is good enough for the tutorial purpose

 unzip all read files 

使用bowtie需要構建索引。

 Unix systems
 bowtie-build arabidopsis.seq arabidopsis

 Windows
 C:\Prog\bowtie-0.12.7\bowtie-build.exe arabidopsis.seq arabidopsis

或者，如果您使用的是支援的生物體之一，您可以從bowtie網站下載提供的索引檔案。

現在您可以使用bowtie比對reads。 Bowtie 具有許多選項，您最好檢查它們。這裡我們告訴它使用兩個處理器（-p 2）來報告基於SAM的比對結果 -S，所有（-a）比對最多隻有一個錯配（-v 1）。我們進一步限制了這一點，因為只有每個reads提供多個有效比對結果才應該被丟棄（-m 1）。

 bowtie -p 2 -S -a -m 1 -v 1 arabidopsis SRR074262.fastq aligned_074262
 bowtie -p 2 -S -a -m 1 -v 1 arabidopsis SRR074263.fastq aligned_074263
 bowtie -p 2 -S -a -m 1 -v 1 arabidopsis SRR074264.fastq aligned_074264
 bowtie -p 2 -S -a -m 1 -v 1 arabidopsis SRR074284.fastq aligned_074284
 bowtie -p 2 -S -a -m 1 -v 1 arabidopsis SRR074285.fastq aligned_074285
 bowtie -p 2 -S -a -m 1 -v 1 arabidopsis SRR074286.fastq aligned_074286

 C:\Prog\bowtie-0.12.7\bowtie.exe -p 2 -S -a -m 1 -v 1 arabidopsis SRR074262.fastq aligned_074262
 

這需要對所有6個庫重複。

輸出SRR074262.fastq

1. 已處理的reads：9619406
2. 至少有一個報告的比對結果的reads：5034111（52.33%）
3. 無法比對的reads：4132070（42.96%）
4. 由於-m被抑制的比對結果的reads：453225（4.71%）

報告了5034111個比對結果到1個輸出流中。

輸出SRR074263.fastq

1. 已處理的reads：4199495
2. 至少有一個報告的比對結果的reads：2298490（54.73%）
3. 無法比對的reads：1701822（40.52%）
4. 由於-m被抑制的比對結果的reads：199183（4.74%）

報告了2298490個比對結果到1個輸出流中。

輸出SRR074264.fastq

1. 已處理的reads：4810892
2. 至少有一個報告的比對結果的reads：2730242（56.75%）
3. 無法比對的reads：1852040（38.50%）
4. 由於-m被抑制的比對結果的reads：228610（4.75%）

報告了2730242個比對結果到1個輸出流中。

輸出SRR074284.fastq

1. 已處理的reads：4763394
2. 至少有一個報告的比對結果的reads：2622118（55.05%）
3. 無法比對的reads：1908441（40.06%）
4. 由於-m被抑制的比對結果的reads：232835（4.89%）

報告了2622118個比對結果到1個輸出流中。

輸出SRR074285.fastq

1. 已處理的reads：8217239
2. 至少有一個報告的比對結果的reads：4515790（54.96%）
3. 無法比對的reads：3206578（39.02%）
4. 由於-m被抑制的比對結果的reads：494871（6.02%）

報告了4515790個比對結果到1個輸出流中。

輸出SRR074286.fastq

1. 已處理的reads：3776979
2. 至少有一個報告的比對結果的reads：2014498（53.34%）
3. 無法比對的reads：1506142（39.88%）
4. 由於-m被抑制的比對結果的reads：256339（6.79%）

報告了2014498個比對結果到1個輸出流中。

彙總Reads

如果您想使用HTSeq-count 並且使用的是Fedora或CentOS，您將需要付出一些額外的努力，因為HTSeq使用的是python 2.6，而Fedora和CentOS只安裝了python 2.4。（在執行以下任何操作之前，請諮詢您的系統管理員是否允許這樣做）

這裡有一些建議 http://stackoverflow.com/questions/1465036/install-python-2-6-in-centos 我最喜歡的是使用EPEL（見下文），但YMMV。

 rpm -Uvh http://download.fedora.redhat.com/pub/epel/5/i386/epel-release-5-4.noarch.rpm
 yum install python26
 yum install python26-numpy
 yum install python26-numpy-devel

然後按照網站上的詳細說明繼續，但鍵入

 python26 setup.py build
 python26 setup.py install

而不是

 python setup.py build
 python setup.py install

---

現在我們準備好下載gff檔案了： ftp://ftp.arabidopsis.org/home/tair/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.gff

該檔案不包含RNA基因和轉座因子。

現在，父ID可能適用於HTSeq的一些用例。如果您需要一個真正的GTF檔案，您可以使用提到的perl指令碼 這裡 底部。如果我們將其命名為celeste_script.pl，這裡將是執行的內容。（我們不能使用 -i Parent 的原因是，否則外顯子可能屬於兩個剪接變體，例如 AT1G01040.1 和 AT1G01040.2，而不會被計入，因為我們想使用交集嚴格模式）

 cat TAIR10_GFF3_genes.gff |perl celeste_script.pl > TAIR10_GTF.gtf 

現在終於

 htseq-count -m intersection-strict -s no aligned_074262 TAIR10_GTF.gff > counts_074262
 htseq-count -m intersection-strict -s no aligned_074263 TAIR10_GTF.gtf > counts_074263
 htseq-count -m intersection-strict -s no aligned_074264 TAIR10_GTF.gtf > counts_074264
 htseq-count -m intersection-strict -s no aligned_074284 TAIR10_GTF.gtf > counts_074284
 htseq-count -m intersection-strict -s no aligned_074285 TAIR10_GTF.gtf > counts_074285
 htseq-count -m intersection-strict -s no aligned_074286 TAIR10_GTF.gtf > counts_074286

現在我們開始得到一些我們可以使用的東西。

 head counts_074262

 AT1G01010       75
 AT1G01020       73
 AT1G01030       33
 AT1G01040       109
 AT1G01046       0
 AT1G01050       41
 AT1G01060       10
 AT1G01070       2
 AT1G01073       0
 AT1G01080       281

 head counts_074263

 AT1G01010       35
 AT1G01020       34
 AT1G01030       18
 AT1G01040       52
 AT1G01046       0
 AT1G01050       51
 AT1G01060       5
 AT1G01070       19
 AT1G01073       0
 AT1G01080       212

 head counts_074264

 AT1G01010       80
 AT1G01020       46
 AT1G01030       45
 AT1G01040       50
 AT1G01046       0
 AT1G01050       69
 AT1G01060       13
 AT1G01070       35
 AT1G01073       0
 AT1G01080       226

 head counts_074284

 AT1G01010       45
 AT1G01020       50
 AT1G01030       33
 AT1G01040       65
 AT1G01046       0
 AT1G01050       50
 AT1G01060       0
 AT1G01070       7
 AT1G01073       0
 AT1G01080       178

 head counts_074285

 AT1G01010       35
 AT1G01020       43
 AT1G01030       26
 AT1G01040       72
 AT1G01046       0
 AT1G01050       51
 AT1G01060       3
 AT1G01070       6
 AT1G01073       0
 AT1G01080       531

 head counts_074286

 AT1G01010       61
 AT1G01020       27
 AT1G01030       36
 AT1G01040       25
 AT1G01046       0
 AT1G01050       25
 AT1G01060       16
 AT1G01070       20
 AT1G01073       0
 AT1G01080       128

差異表達基因的分析

啟動R

使用DESeq（現在參考 vignette）

 R
 

在R中

 library(DESeq)
 setwd("where the data is")
 c1<-read.table("counts_074284",row.names=1) 
 c2<-read.table("counts_074286",row.names=1)
 c3<-read.table("counts_074262",row.names=1)
 c4<-read.table("counts_074263",row.names=1)
 c5<-read.table("counts_074264",row.names=1)
 c6<-read.table("counts_074285",row.names=1)
 counts<-cbind(c1,c2,c3,c4,c5,c6)
 counts<-counts[-c(32679:32683),]  #remove the more general lines
 colnames(counts)<-c("P1","P2","P3","M1","M2","M3")
 design <- rep (c("P","Mo"),each=3)
 de  <-  newCountDataSet(counts, design)
 de  <-  estimateSizeFactors( de)
 # de  <-  estimateVarianceFunctions(  de  ) # This function has been removed. Use 'estimateDispersions' instead.  
 de <- estimateDispersions( de )
 res  <-  nbinomTest(  de,  "P",  "Mo")

有多少基因是顯著的？

 sum(na.omit(res$padj<0.05))