本文是由比利時列日大學(xué)Marc HANIKENNE課程整理。陸陸續(xù)續(xù)花費1個月完成，根據(jù)需要做的主要內(nèi)容分為四個部分。
第一部分：將RNA-seq數(shù)據(jù)映射到參考基因組上
第二部分：將RNA-seq數(shù)據(jù)映射到參考轉(zhuǎn)錄組上，并且生成基因表達矩陣，用于第三部分分析
第三部分：使用DESeq包鑒定差異表達基因（成對），
第四部分：對差異基因進行后續(xù)的GO和KEGG注釋

1 目標

RNA-Seq 模塊的目標是說明如何處理和分析 RNA-Seq 數(shù)據(jù)以識別差異表達基因 (DGE)。
練習(xí)中使用真實數(shù)據(jù)集，來自暴露于兩種生長條件的擬南芥的兩種基因型的 Illumina RNA 測序。
需要做：
1). 在參考基因組（工具：tophat 和 htseq-count）或參考轉(zhuǎn)錄組（工具trinity）上映射reads，作為reads映射和計數(shù)的兩種相互替代策略；
2). 使用 DESeq2包（ R 語言）鑒定不同表達的基因（DGEs）
3). 簡單進行數(shù)據(jù)挖掘(GO和KEGG注釋)。

2 數(shù)據(jù)介紹

模式植物擬南芥的兩種基因型（wt 和突變體）在對照 (c) 和處理 (t) 條件下生長。 在實驗結(jié)束時，分別收獲根 (r) 和芽 (s) 組織（每個樣品 4 株植物）。 本實驗獨立重復(fù) 3 次。 使用 NextSeq Illumina 測序儀以高通量模式（~400 Mio. clusters）在兩次運行中制備和測序總 RNA，以獲得 75 bp 的單端讀數(shù)。 注意，已經(jīng)使用 fastqc 和trimmomatic 分析和質(zhì)控了原始reads（這兩個軟件我會在基因組組裝分析中介紹）。

3 分析數(shù)據(jù)

3.1 查看數(shù)據(jù)

head <your_sample>.fastq
查看每個文件數(shù)據(jù)的reads數(shù)：

image.png

將所有的樣本名字放在一個文件（samples.ids）中, 這樣方便后續(xù)進行批量處理。
shell code

for f in *.fastq; do echo `basename $f .fastq`; done > samples.ids

3.2 RNA-Seq數(shù)據(jù)分析中reads映射的一般考慮

在分析 RNA-Seq 數(shù)據(jù)時，有兩種通用策略可以在計數(shù)之前映射reads。
(i) Reads可以使用splice aware比對工具映射到參考基因組上：比對器將允許拆分讀數(shù)以跨越內(nèi)含子。當感興趣的物種有參考基因組可用時，這種策略通常受到青睞，因為它允許檢測以前沒有描述過的轉(zhuǎn)錄本或轉(zhuǎn)錄本變異。然而，這種方法僅限于（相對）少數(shù)有高質(zhì)量參考基因組可用的物種。
(ii) Reads也可以映射到參考轉(zhuǎn)錄組上。比對更容易，并且當目標物種不存在參考轉(zhuǎn)錄組時，可以使用 RNA-Seq 讀取組裝自定義參考轉(zhuǎn)錄組，使其幾乎適用于所有物種。然而，這種方法僅限于檢測已知轉(zhuǎn)錄本和轉(zhuǎn)錄本變體的表達。

第一部分：3.3 Read映射到參考基因組

3.3.1 工具介紹

1. tophat 軟件

我們將使用流行的 tophat，這是一個將 RNA-Seq 讀數(shù)與基因組對齊以識別外顯子-外顯子剪接連接的程序。 它建立在超快的短讀映射程序 bowtie2 之上。 該軟件針對 75bp 或更長的讀取進行了優(yōu)化。
更多介紹點擊查看：Tophat 鏈接。

TopHat 如何找到連接點的原理：

TopHat 可以在沒有參考注釋的情況下找到拼接點。通過首先將 RNA-Seq的reads映射到基因組，TopHat可以識別出潛在的外顯子，這是因為許多 RNA-Seq reads將與基因組連續(xù)對齊。使用這個初始映射信息，TopHat 建立一個可能的剪接連接的數(shù)據(jù)庫，然后將讀取映射到這些連接以確認它們。

短read測序機目前可以產(chǎn)生 100 bp 或更長的reads，但許多外顯子比這短，因此它們會在初始映射中被遺漏。TopHat 主要通過將所有輸入reads分成更小的段來解決這個問題，然后獨立映射。這些小段對齊在程序的最后一步重新組合在一起，以產(chǎn)生端到端的read對齊。

TopHat 從兩個證據(jù)來源生成可能的剪接點數(shù)據(jù)庫。第一個也是最有力的證據(jù): 剪接點來源是當來自同一reads的兩個片段（對于至少 45bp 的reads）在同一基因組序列上的某個距離處作圖時，或者當一個內(nèi)部片段未能作圖時——再次表明這種讀取跨越多個外顯子。使用這種方法，“GT-AG”、“GC-AG”和“AT-AC”內(nèi)含子將從頭開始找到。第二個來源是“coverage islands,”的配對，它們是初始映射中堆積讀取的不同區(qū)域。轉(zhuǎn)錄組中的相鄰島通常拼接在一起，因此 TopHat 尋找將這些與內(nèi)含子連接起來的方法。我們只建議用戶將第二個選項 (--coverage-search) 用于短讀取 (< 45bp) 和少量讀取 (<= 1000 萬)。后一個選項只會報告“GT-AG”內(nèi)含子之間的比對.

Tophat可以使用FASTA,FASTQ(推薦)格式的reads。

想要使用這個軟件，首先需要使用一下命令：

image.png

Bowtie 2是一種超快且內(nèi)存高效的工具，用于將測序讀數(shù)與長參考序列對齊。它特別擅長將大約 50 到 100 個字符的讀數(shù)與相對較長的（例如哺乳動物）基因組進行比對。Bowtie 2 使用FM 索引（基于Burrows-Wheeler 變換或BWT）對基因組進行索引，以保持其較小的內(nèi)存占用：對于人類基因組，其內(nèi)存占用通常約為 3.2 GB 的 RAM。Bowtie 2 支持間隙、局部和雙端對齊模式?？梢酝瑫r使用多個處理器以實現(xiàn)更高的對齊速度。

Bowtie 2 以SAM格式輸出對齊方式，從而能夠與大量使用同樣使用SAM文件格式的其他工具（例如SAMtools、GATK）進行互操作。Bowtie 2 在GPLv3 許可下分發(fā)，它在 Windows、Mac OS X 和 Linux 和 BSD 下的命令行上運行。

Bowtie 2通常是比較基因組學(xué)管道的第一步，包括變異調(diào)用、ChIP-seq、RNA-seq、BS-seq。Bowtie 2和Bowtie（此處也稱為“ Bowtie 1 ”）也緊密集成到許多其他工具中，此處列出了其中一些。

要找到與 Tophat 的連接點，您首先需要為 RNA-Seq 實驗中的organism安裝bowtie index。 bowtie 站點為人類、小鼠、果蠅等提供了預(yù)先構(gòu)建的bowtie index。 如果你的organism沒有bowtie index，使用 bowtie2-build 很容易自己構(gòu)建一個.

image.png

Bowtie2-inspect 從 bowtie index 中提取信息，允許確定它是什么類型的索引以及使用什么參考序列來構(gòu)建它。

2. GFF/GTF 格式文件

通過提供描述基因組特征（例如外顯子/內(nèi)含子）染色體位置的基因組注釋文件，可以幫助通過 tophat 在參考基因組上進行reads映射。 注釋文件以 GFF/GTF 格式提供。

Tophat使用的的genome annotation文件就是GFF/GTF格式。
GFF一般就9列，是以Tab間隔，對應(yīng)名字：

image.png

GTF(general transfer format)是GFF第二個版本，

3 htseq-count軟件

給定一個具有對齊測序讀數(shù)的文件和基因組特征列表，htseq-count會計算有多少reads映射到每個特征。這里的特征是染色體上的區(qū)間（即，位置范圍）或這些區(qū)間的聯(lián)合。在 RNA-Seq 的情況下，特征通常是基因，其中每個基因在這里被視為其所有外顯子的聯(lián)合。人們也可以將每個外顯子視為一個特征，例如，為了檢查可變剪接。對于比較 ChIP-Seq，特征可能是預(yù)定列表中的結(jié)合區(qū)域。

必須特別注意決定如何處理重疊多個特征的reads。 htseq-count 腳本允許在三種模式之間進行選擇。 htseq-count 的三種重疊解析模式的工作原理如下：對于 read 中的每個位置 i，定義一個集合 S(i) 為所有與位置 i 重疊的特征的集合。 然后，考慮集合 S，它是（i 遍歷讀取或讀取對中的所有位置）

• 模式并集， 取所有集合 S(i) 的并集。 對于大多數(shù)用例，建議使用此模式。
• 模式交集，嚴格的所有集合 S(i) 的交集。

• 模式交集， 非空的所有非空集 S(i) 的交集。
  如果 S 恰好包含一個特征，則針對該特征計算讀取（或讀取對）。 如果它包含多個特征，則讀?。ɑ蜃x取對）計為不明確的（不計入任何特征），如果 S 為空，則讀?。ɑ蜃x取對）計為 no_feature。
  看圖更清晰理解：

  
  
  image.png

3.3.2 下載擬南芥參考基因

網(wǎng)站：https://www.araport.org/（需要注冊）
也可以使用以下命令：

curl -sO -H 'Authorization: Bearer <your_id_key>' \
https://api.araport.org/files/v2/media/system/araport-public-files// \
TAIR10_genome_release/assembly/TAIR10_Chr.all.fasta.gz
curl -sO -H 'Authorization: Bearer <your_id_key>' \
https://api.araport.org/files/v2/media/system/araport-public-files// \
Araport11_latest/annotation/Araport11_GFF3_genes_transposons.201606.gff.gz
curl -sO -H 'Authorization: Bearer <your_id_key>' \
https://api.araport.org/files/v2/media/system/araport-public-files// \
Araport11_latest/annotation/Araport11_GFF3_genes_transposons.201606.gtf.gz

3.3.3 給參考基因建index

使用bowtie2-build.

To index the Arabidopsis,花費2分鐘

bowtie2-build Arabidopsis.fasta At_ref

檢查index，花費幾秒鐘

bowtie2-inspect -n At-ref

3.3.4 reads映射

reads是存在以逗號隔開的FASTQ或FASTA格式文件中

使用tophat完成

一般使用命令：

image.png

更多的選擇閱讀文檔

其中： --num-threads 4 ## 可以多線程運行
-o --output-dir <string> ## tophat輸出結(jié)果的文目錄
--min-intron-length <int> ##內(nèi)含子的最短長度：默認70
--max-intron-length <int>##內(nèi)含子的最長長度：默認500000
-G --GTF <GTF/GFF3 file> # 為 TopHat 提供一組基因模型注釋和/或已知轉(zhuǎn)錄本，作為 GTF 2.2 或 GFF3 格式的文件。 如果提供此選項，TopHat 將首先提取轉(zhuǎn)錄本序列，并首先使用 Bowtie 將讀數(shù)與該虛擬轉(zhuǎn)錄組對齊。 只有沒有完全映射到轉(zhuǎn)錄組的讀數(shù)才會被映射到基因組上。 在轉(zhuǎn)錄組上進行映射的讀數(shù)將被轉(zhuǎn)換為基因組映射（根據(jù)需要拼接）并與最終 tophat 輸出中的新映射和連接點合并。

請注意，所提供的 GTF/GFF 文件的第一列（指示特征所在的染色體或重疊群的列）中的值必須與用于 TopHat 的 Bowtie 索引中的參考序列名稱相匹配。 你可以使用 bowtie-inspect 進行檢查。
--transcriptome-index <dir/prefix> ##當為 TopHat 提供已知的轉(zhuǎn)錄文件（上面的 -G/--GTF 選項）時，會構(gòu)建轉(zhuǎn)錄組序列文件，并且必須為其創(chuàng)建 Bowtie index，以便將讀取與已知轉(zhuǎn)錄本對齊。創(chuàng)建此 Bowtie 索引可能非常耗時，并且在許多情況下，相同的轉(zhuǎn)錄組數(shù)據(jù)被用于將多個樣本與 TopHat 對齊。因此，轉(zhuǎn)錄組索引和相關(guān)數(shù)據(jù)文件（原始 GFF 文件）可以在使用此選項的多次 TopHat 運行中重復(fù)使用，因此這些文件僅針對具有給定轉(zhuǎn)錄本數(shù)據(jù)的第一次運行創(chuàng)建。如果計劃使用相同的轉(zhuǎn)錄組數(shù)據(jù)進行多次 TopHat 運行，則應(yīng)首先使用 -G/--GTF 選項以及指向目錄和名稱前綴的 --transcriptome-index 選項運行 TopHat，該選項將指示轉(zhuǎn)錄組數(shù)據(jù)的位置文件將被存儲。然后使用相同的 --transcriptome-index 選項值的后續(xù) TopHat 運行將直接使用在第一次運行中創(chuàng)建的轉(zhuǎn)錄組數(shù)據(jù)（第一次運行后不需要 -G 選項）。

開始操作

軟鏈接參考基因組的FASTA：

ln -s Arabidopsis.fasta At_ref.fa

創(chuàng)建轉(zhuǎn)錄組的index.只需要做一次，即可用于全部樣本，耗時5分鐘

tophat -G Arabidopsis.gtf --transcriptome-index=transcriptome_data/At_ref At_ref

會在transcriptome_data/ 下產(chǎn)生10個文件

映射reads， 先創(chuàng)建一個模板

tophat -o output_[% basename %] --read-mismatches 2 --min-intron-length 40 \
--max-intron-length 2000 --num-threads 2 --report-secondary-alignments \
--no-novel-juncs --transcriptome-index=transcriptome_data/At_ref At_ref \
[% basename %].fastq

每個樣本創(chuàng)一個sh

for f in `cat samples.ids`
do tpage --define queue=smallnodes --define basename=$f tophat.tt \
> tophat_$f.sh
done

提交任務(wù)：

for f in `cat samples.ids`
do qsub -pe snode 2 tophat_$f.sh
done

此步驟花費大約1個小時
查看任務(wù)

qstat -f

對所有的樣本進行summary查看

for f in `cat samples.ids`
do head output_$f/align_summary.txt
done

3.3.5 Reads counting

使用htseq-count

image.png

image.png

-s <yes/no/reverse>
數(shù)據(jù)是否來自特定鏈的檢測（默認值：yes）.對于stranded=no，無論是映射到與特征相同還是相反的鏈，都認為讀取與特征重疊。 對于 stranded=yes 和單端讀取，讀取必須映射到與特征相同的鏈。 對于雙端讀取，第一次讀取必須在同一條鏈上，第二次讀取必須在相反的鏈上。 對于stranded=reverse，這些規(guī)則是相反的。

reads計數(shù)模板

htseq-count -f bam -s reverse output_[% basename %]/accepted_hits.bam Arabidopsis.gtf

運行花費半個小時。

檢索對齊統(tǒng)計信息

shell 命令

for f in <your_name>_htseqcount_*.o*; do tail -n 5 $f; done

. 組裝讀取計數(shù)矩陣

得到基因名字

cut -f1 <your_name>_htseqcount_<your_sample>.o<job_number> > gene_lists

得到 read count

for f in `cat samples.ids`
do cut -f2 <your_name>_htseqcount_$f.o* > $f.count
done

組裝基因 list和count

paste gene_lists *.count > <your_name>_htseqcount.matrix

得到的這個結(jié)果文件，將用于后續(xù)的統(tǒng)計分析，發(fā)現(xiàn)DGEs

第二部分： 3.4 read映射到參考轉(zhuǎn)錄組

3.4.1 工具介紹

1. trinity 軟件
   由Broad 研究所和耶路撒冷希伯來大學(xué)開發(fā)的 Trinity代表了一種從 RNA-seq 數(shù)據(jù)高效、穩(wěn)健地從頭重建轉(zhuǎn)錄組的新方法。Trinity 結(jié)合了三個獨立的軟件模塊：Inchworm、Chrysalis 和 Butterfly，依次應(yīng)用于處理大量 RNA-seq 讀數(shù)。Trinity 將序列數(shù)據(jù)劃分為許多單獨的 de Bruijn 圖，每個圖代表給定基因或基因座的轉(zhuǎn)錄復(fù)雜性，然后獨立處理每個圖以提取全長剪接異構(gòu)體并梳理源自旁系同源基因的轉(zhuǎn)錄本。簡而言之，過程是這樣工作的：

• Inchworm將 RNA-seq 數(shù)據(jù)組裝成獨特的轉(zhuǎn)錄本序列，通常為顯性同種型生成全長轉(zhuǎn)錄本，但隨后僅報告選擇性剪接轉(zhuǎn)錄本的獨特部分。

• Chrysalis將 Inchworm contigs 聚類成簇，并為每個簇構(gòu)建完整的 de Bruijn 圖。每個簇代表給定基因（或共享序列的基因組）的完整轉(zhuǎn)錄復(fù)雜性。然后 Chrysalis 在這些不相交的圖之間劃分完整的讀取集。

• Butterfly然后并行處理各個圖，跟蹤圖中讀取和讀取對的路徑，最終報告可變剪接同種型的全長轉(zhuǎn)錄本，并梳理出對應(yīng)于旁系同源基因的轉(zhuǎn)錄本。

2 轉(zhuǎn)錄組組裝下游分析
組裝完成后，您可能需要進行多項分析，以根據(jù)組裝的轉(zhuǎn)錄本和輸入的 RNA-Seq 數(shù)據(jù)探索生物體生物學(xué)的各個方面。此類研究可能包括：

• 量化轉(zhuǎn)錄本和基因豐度。這是許多其他分析的先決條件，例如檢查樣本中差異表達的轉(zhuǎn)錄本。

• 質(zhì)量檢查您的樣品和生物復(fù)制品。確保您的樣本和生物學(xué)重復(fù)之間的關(guān)系有意義。如果數(shù)據(jù)存在任何混雜因素，例如異常值或批處理效應(yīng)，您將希望盡早發(fā)現(xiàn)它們并在進一步的數(shù)據(jù)探索中考慮到它們。

• 進行差異表達分析。Trinity 直接支持多種 DE 分析方法，包括 edgeR、DESeq2、Limma/Voom 和 ROTS。

• 使用提取的編碼區(qū)TransDecoder和功能注釋使用的成績單Trinotate。

• 如果您的生物體具有組裝的基因組，請考慮使用 Trinity 轉(zhuǎn)錄組組裝使用PASA進行基因結(jié)構(gòu)注釋。

分析使用 perl 腳本嵌入了一組流行工具進行分析。所以我們將使用一些代碼來映射、對齊和計數(shù)reads.
首先使用：align_and_estimate_abundance.pl， 功能是index參考轉(zhuǎn)錄組，對齊和計數(shù)reads。其中index和對其依靠的是Bowtie2軟件，估計轉(zhuǎn)錄的豐富度使用RSEM。從 RNA-Seq 數(shù)據(jù)進行轉(zhuǎn)錄量化的一個關(guān)鍵挑戰(zhàn)是處理映射到多個基因或同種型的讀數(shù)。 這個問題對于在沒有測序基因組的情況下用從頭轉(zhuǎn)錄組組裝進行量化特別重要，因為很難確定哪些轉(zhuǎn)錄本是同一基因的同種型。 RSEM 允許從單端或雙端 RNA-Seq 數(shù)據(jù)中量化基因和同種型豐度，無需測序基因組。

請注意，Trinity提供了依賴對齊和無對齊的讀取計數(shù)兩種方案。 它提供基因和轉(zhuǎn)錄水平的reads計數(shù)估計。

3 數(shù)據(jù)標準化
原始讀取計數(shù)必須針對許多偏差（例如基因長度或每個樣本的讀取數(shù)）進行校正。 Trinity 提供的任何豐度估計方法都將提供源自每個轉(zhuǎn)錄本的 RNA-Seq 片段計數(shù)的轉(zhuǎn)錄本水平估計，此外還提供了轉(zhuǎn)錄本表達的標準化測量，該測量值考慮了轉(zhuǎn)錄本長度、映射的讀取數(shù) 到轉(zhuǎn)錄本，以及映射到任何轉(zhuǎn)錄本的讀取總數(shù)。 標準化表達指標可以報告為每千堿基轉(zhuǎn)錄本長度的片段映射 (FPKM) 或每百萬轉(zhuǎn)錄本 (TPM) 的轉(zhuǎn)錄本。

3.4.2 擬南芥參考轉(zhuǎn)錄組

來自Araport, 需要登錄進行免費的注冊。再使用以下代碼獲取。

curl -sO -H 'Authorization: Bearer <your_id_key>' \
https://api.araport.org/files/v2/media/system/araport-public-files// \
Araport11_Release_201606/annotation/Araport11_genes.201606.cds.fasta.gz

3.4.3 indexing擬南芥參考轉(zhuǎn)錄組

使用ltrinity的perl命令：align_and_estimate_abundance.pl，可以對所有樣本一次完成。

image.png

更多options:https://github.com/trinityrnaseq/trinityrnaseq/wiki/Trinity-Transcript-Quantification#estimating-transcript-abundance

index操作的命令

perl /media/vol1/apps/trinityrnaseq-2.2.0/util/align_and_estimate_abundance.pl \
--transcripts Arabidopsis_transcripts.fasta --est_method RSEM --aln_method bowtie2 \
--prep_reference --output_dir ref_transcriptome_index

此過程花費大約5分鐘，會生成14個文件，具有.bowtie2.和.RSEM文件

3.4.4 對read進行映射和計數(shù)

使用ltrinity的perl命令：align_and_estimate_abundance.pl，并且使用RSEM估計方法

image.png

image.png

2建立gene_trans_map
我們需要準備一個gene和轉(zhuǎn)錄名字對應(yīng)的文件，并且需要將轉(zhuǎn)錄組的fasta文件排序成參考轉(zhuǎn)錄組的方式
從fasta文件中提取轉(zhuǎn)錄的ID， shell命令

grep \> Arabidopsis_transcripts.fasta | cut -f2 -d '>' | cut -f1 -d '|' > transcripts.ids
# Let's paste twice this list in the same file
$ paste transcripts.ids transcripts.ids > double_transcripts.ids
$ head double_transcripts.ids
# And apply the following perl one liner to remove the transcript number
# from 1st column
$ perl -nle 's/^(AT\w+)\.\d+/$1/g; print' double_transcripts.ids > gene_trans_map.txt

3進行read的map和count

align_and_estimate_abundance.pl命令

使用模板：

perl /media/vol1/apps/trinityrnaseq-2.2.0/util/align_and_estimate_abundance.pl \
--transcripts Arabidopsis_transcripts.fasta --seqType fq \
--single [% basename %].fastq --est_method RSEM --aln_method bowtie2 \
--SS_lib_type R --thread_count [% thread %] \
--gene_trans_map gene_trans_map.txt --output_prefix [% basename %] \
--output_dir trinity_[% basename %]

創(chuàng)建多個樣本的sh文件：

for f in `cat samples.ids`
do tpage --define queue=smallnodes --define basename=$f \
--define thread=2 trinity_align_estimate.tt > align_estimate_$f.sh
done

批量提交任務(wù)：

for f in `cat samples.ids`
do qsub -pe snode 2 align_estimate_$f.sh
done

這一步大概花費90分鐘
再查看你的結(jié)果：

image.png

image.png

3.4.5 生成表達矩陣

使用：trinity下的abundance_estimates_to_matrix.pl命令
該腳本將非常簡單地創(chuàng)建一個矩陣，將所有樣本的讀取計數(shù)數(shù)據(jù)組合在一起。
使用：

perl /media/vol1/apps/trinityrnaseq-2.2.0/util/abundance_estimates_to_matrix.pl \
--est_method RSEM trinity_*/*.genes.results --out_prefix <your_name>

大概需要2分鐘

該腳本輸出多個文件。 <your_name>.counts.matrix 文件提供未經(jīng)標準化的原始讀取計數(shù)。 這是必須用于使用 DESeq2 進行進一步 DGE 分析的文件。
額外的 .matrix 文件對應(yīng)于 TPM 表達值矩陣（未跨樣本歸一化）和 TMM 歸一化表達值矩陣（應(yīng)用了跨樣本歸一化）。 有關(guān)此查看的更多詳細信息：https://github.com/trinityrnaseq/trinityrnaseq/wiki/Trinity-Transcript-Quantification

第三部分： 3.5 鑒定差異表達的基因

使用R 包DESeq2。

3.5.1 包介紹

詳細文檔介紹：https://bioconductor.org/packages/release/bioc/html/DESeq2.html.
DESeq2 允許估計來自高通量測序分析的計數(shù)數(shù)據(jù)（依賴于方差-均值），并基于使用負二項分布的廣義線性模型 (GLM) 測試差異表達。 一個非?；镜?GLM 模型如下：

image.png

DESeq2 將首先對數(shù)據(jù)進行建模以估計模型的系數(shù)。 一旦設(shè)置了系數(shù)，就可以為每個 x 確定 y。
這里只是簡要介紹使用的函數(shù)，其他的請看上文的連接文檔
1 DESeqDataSetFromMatrix
DESeqDataSet 是 RangedSummarizedExperiment 的子函數(shù)，用于存儲輸入值、中間計算和差異表達式分析的結(jié)果。 DESeqDataSet 在“計數(shù)”矩陣中強制執(zhí)行非負整數(shù)值，作為分析列表中的第一個元素存儲。 此外，必須提供指定實驗設(shè)計的公式。
使用：DESeqDataSetFromMatrix(countData, colData, design)
其中design:設(shè)計一個公式來表達每個基因的計數(shù)如何取決于 colData 中的變量。 許多 R 公式是有效的，包括具有多個變量的設(shè)計，例如 ~ 組 + 條件，以及具有交互作用的設(shè)計，例如 ~ 基因型 + 治療 + 基因型：治療。 查看各種設(shè)計的結(jié)果以及如何提取結(jié)果表。
countData,為輸入的非負整數(shù)矩陣
colData: 輸入DataFrame或者data.frame格式，其只是有一列，行是對應(yīng)著countData。
2 DESeq
DESeq 基于負二項分布進行差異表達分析。 它通過以下步驟執(zhí)行默認分析：
? 估計大小因素：estimateSizeFactors
? 估計色散：estimateDispersions
? 負二項式 GLM 擬合和 Wald 統(tǒng)計：nbinomWaldTest

有關(guān)每個步驟的完整詳細信息，請參閱相應(yīng)功能的手冊頁。 DESeq 函數(shù)返回 DESeqDataSet 對象后，可以使用結(jié)果函數(shù)生成結(jié)果表（log2 倍數(shù)變化和 p 值）。 有關(guān)針對多重測試校正的獨立過濾和 p 值調(diào)整的信息，請參閱結(jié)果手冊頁。

使用DESeq(object)， 是一個DESeqDataSet的object. 如： DESeqDataSetFromMatrix.

3 results
結(jié)果從 DESeq 分析中提取結(jié)果表，給出樣本的基本均值、log2 倍數(shù)變化、標準誤差、檢驗統(tǒng)計量、p 值和調(diào)整后的 p 值。

resultsNames 返回模型的估計效應(yīng)（系數(shù)）的名稱。
使用

results(object, contrast, lfcThreshold = 0, alpha = 0.1)
resultsNames(object)

object是 DESeqDataSet，已在其上調(diào)用以下函數(shù)之一：DESeq、nbinocontrastmWaldTest 或 nbinomLRT。
contrast參數(shù)指定從對象中提取什么比較以構(gòu)建結(jié)果表。
在我們的例子中，一個字符向量正好包含三個元素：設(shè)計公式中一個因子的名稱、折疊變化的分子級別的名稱以及折疊變化的分母級別的名稱（最簡單的情況）。
lfcThreshold 是一個非負值，指定 log2 倍變化閾值。 默認值為 0，對應(yīng) log2 倍數(shù)變化為零的測試。

alpha 用于優(yōu)化獨立過濾的顯著性截止值（默認為 0.1）。 如果調(diào)整后的 p 值截止值 (FDR) 不是 0.1，則 alpha 應(yīng)設(shè)置為該值。
plotCounts
plotCounts 允許在對數(shù)尺度上為單個基因創(chuàng)建歸一化計數(shù)圖。
使用：plotCounts(dds, gene, intgroup = "condition")
dds是DESeqDataSet.， gene是一個特殊基因名稱，intgroup:在colData(x)中，用于進行分組的名稱。

3.5.2 下載DESeq2

library(BiocManager)
BiocManager::install("openssl")
BiocManager::install("RCurl")
BiocManager::install(c("DESeq2","limma","gplots"), force = T)

3.5.3 鑒定差異表達基因(成對比較)

我們將在下面發(fā)現(xiàn)如何編寫允許這樣做的 R 腳本。 您需要在 RStudio 內(nèi)的同一個 R 腳本中按順序添加每個新步驟。 我們將首先基于獨立于治療的基因型（WT vs 突變體）識別顯示 DGE 的基因，然后基于獨立于基因型的治療（Ctrl vs Treat）的 DGE，最后看看治療對 每個基因型。 這些分析中的每一個都對應(yīng)于統(tǒng)計測試的不同設(shè)計，在我們的腳本中必須考慮到這一點。

Step 1. Load the data and describe the dataset

#Load data
countData=read.table("tophat_root.matrix",header=TRUE,row.names=1,sep="\t")
head(countData)
#Describe the dataset for each variable
genot=rep(c("WT","mut"),each=6)
treat=(rep(rep(c("Ctrl","Treat"),each=3),2))
g_t=rep(c("WT-Ctrl", "WT-Treat", "mut-Ctrl", "mut-Treat"),each=3)
#Load dataset description in a data frame
colData=data.frame(g_t,genot,treat,row.names=names(countData))
colData

Step 2. Build model for genotype effect analysis

#Genotype effect
#####
#Load data using the DESeqDataSetFromMatrix command
genotDesign=DESeqDataSetFromMatrix(countData = countData,colData = colData,
                                   design = ~ genot)
#Build model using the DESeq command
genot_DESeq <- DESeq(genotDesign)
#Observe parameters of the model
resultsNames(genot_DESeq)

Step 3. Summary statistics of the data with PCA

rld<-rlog(genot_DESeq)
#tiff(filename = "PCA_genot.tiff", width = 1500, height = 1500, units = "px", res = 150)
plotPCA(rld, intgroup=c("g_t"))
dev.off()

Step 4. Build a heatmap of sample distances

#Build sample distance
sampleDist <- dist(t(assay(rld)))
#Build heatmap
sampleDistMatrix<-as.matrix(sampleDist)
rownames(sampleDistMatrix)<-paste(rld$g_t)
colnames(sampleDistMatrix)<-NULL
colours=colorRampPalette(rev(brewer.pal(9, "Blues")))(300)
tiff(filename = "heatmap_sampledist_Treat_root.tiff", width = 1500, 
     height = 1500, units = "px", res = 150)
heatmap.2(sampleDistMatrix, dendrogram = "both", trace = "none", col = colours,
           main = "Treat Root Sample Distance", margin=c(6, 8))
dev.off()

Step 5. Identify DGEs for genotype effect

#Extract results (contrast WT and mutant) with set lfc and pvalue
res_genot=results(genot_DESeq, contrast = c("genot", "mut", "WT"), 
                  lfcThreshold = 1, alpha = 0.05)
#Observe the summary of the analysis
summary(res_genot)
#Look at the results
head(res_genot,2)
#Export data into a table
write.table(res_genot,"pairwise_root_WT_vs_mut.txt",sep="\t")
#Filter data to extract up-regulated genes with a certain lfc and pvalue
fc_genotM<- res_genot[which(res_genot$log2FoldChange > 1 & res_genot$padj<0.05),]
#Filter data to extract down-regulated genes with a certain lfc and pvalue
fc_genotL<- res_genot[which(res_genot$log2FoldChange < -1 & res_genot$padj<0.05),]
#Export data into tables
write.table(fc_genotM,"root_higher_mut_vs_WT.txt",sep="\t")
write.table(fc_genotL,"root_lower_mut_vs_WT.txt",sep="\t")

Step 6. Build graph showing gene expression for genes of interest

plotCounts(genot_DESeq, "AT2G19110", intgroup = "genot")

第四部分： 3.6 數(shù)據(jù)挖掘

我們將對 DGE 數(shù)據(jù)進行非常簡單和基本的數(shù)據(jù)挖掘。 我們將使用 Thalemine 門戶，它實現(xiàn)了一個 Intermine 接口，允許非?？焖俚孬@取有關(guān)基因列表的功能信息，包括 GO 富集分析。連接：
https://bar.utoronto.ca/thalemine/

為了使用這個工具，我們首先需要從DESeq2生成的基因列表中提取基因名稱
DESeq2 分析生成了 24 個文件（8 個成對.txt，8 個高.txt 和 8 個低.txt）。 由于它們對應(yīng)于我們過濾后的數(shù)據(jù)，我們只對 high.txt 和 lower.*txt 文件感興趣。
使用shell對文件信息提取，并且進行合并：

mkdir full_DGE_data
mv pairwise*.txt full_DGE_data
ls
# have a look at one of the files
 head higher_root_Ctrl_mut_vs_WT.txt
cut -f2 -d '"' higher_root_Ctrl_mut_vs_WT.txt | head
cut -f2 -d '"' higher_root_Ctrl_mut_vs_WT.txt | sed '1d' | head
# Let's do that for all files
for f in *root*.txt; do cut -f2 -d '"' $f | sed '1d' > $f.gene.list; done
 ls

使用Thalemine注釋基因