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

1 目標(biāo)

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

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

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

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

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

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

image.png

將所有的樣本名字放在一個(gè)文件（samples.ids）中, 這樣方便后續(xù)進(jìn)行批量處理透乾。
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í)，有兩種通用策略可以在計(jì)數(shù)之前映射reads。
(i) Reads可以使用splice aware比對(duì)工具映射到參考基因組上：比對(duì)器將允許拆分讀數(shù)以跨越內(nèi)含子续徽。當(dāng)感興趣的物種有參考基因組可用時(shí)蚓曼，這種策略通常受到青睞，因?yàn)樗试S檢測(cè)以前沒(méi)有描述過(guò)的轉(zhuǎn)錄本或轉(zhuǎn)錄本變異钦扭。然而纫版，這種方法僅限于（相對(duì)）少數(shù)有高質(zhì)量參考基因組可用的物種。
(ii) Reads也可以映射到參考轉(zhuǎn)錄組上客情。比對(duì)更容易其弊，并且當(dāng)目標(biāo)物種不存在參考轉(zhuǎn)錄組時(shí)，可以使用 RNA-Seq 讀取組裝自定義參考轉(zhuǎn)錄組膀斋，使其幾乎適用于所有物種梭伐。然而，這種方法僅限于檢測(cè)已知轉(zhuǎn)錄本和轉(zhuǎn)錄本變體的表達(dá)仰担。

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

3.3.1 工具介紹

1. tophat 軟件

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

TopHat 如何找到連接點(diǎn)的原理：

TopHat 可以在沒(méi)有參考注釋的情況下找到拼接點(diǎn)贮尉。通過(guò)首先將 RNA-Seq的reads映射到基因組拌滋，TopHat可以識(shí)別出潛在的外顯子，這是因?yàn)樵S多 RNA-Seq reads將與基因組連續(xù)對(duì)齊猜谚。使用這個(gè)初始映射信息败砂，TopHat 建立一個(gè)可能的剪接連接的數(shù)據(jù)庫(kù)，然后將讀取映射到這些連接以確認(rèn)它們魏铅。

短read測(cè)序機(jī)目前可以產(chǎn)生 100 bp 或更長(zhǎng)的reads昌犹，但許多外顯子比這短，因此它們會(huì)在初始映射中被遺漏览芳。TopHat 主要通過(guò)將所有輸入reads分成更小的段來(lái)解決這個(gè)問(wèn)題祭隔，然后獨(dú)立映射。這些小段對(duì)齊在程序的最后一步重新組合在一起路操，以產(chǎn)生端到端的read對(duì)齊疾渴。

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

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

想要使用這個(gè)軟件右钾，首先需要使用一下命令：

image.png

Bowtie 2是一種超快且內(nèi)存高效的工具旱爆，用于將測(cè)序讀數(shù)與長(zhǎng)參考序列對(duì)齊舀射。它特別擅長(zhǎng)將大約 50 到 100 個(gè)字符的讀數(shù)與相對(duì)較長(zhǎng)的（例如哺乳動(dòng)物）基因組進(jìn)行比對(duì)。Bowtie 2 使用FM 索引（基于Burrows-Wheeler 變換或BWT）對(duì)基因組進(jìn)行索引怀伦，以保持其較小的內(nèi)存占用：對(duì)于人類基因組脆烟，其內(nèi)存占用通常約為 3.2 GB 的 RAM。Bowtie 2 支持間隙房待、局部和雙端對(duì)齊模式邢羔。可以同時(shí)使用多個(gè)處理器以實(shí)現(xiàn)更高的對(duì)齊速度桑孩。

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

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

要找到與 Tophat 的連接點(diǎn)猜煮，您首先需要為 RNA-Seq 實(shí)驗(yàn)中的organism安裝bowtie index次员。 bowtie 站點(diǎn)為人類、小鼠王带、果蠅等提供了預(yù)先構(gòu)建的bowtie index淑蔚。 如果你的organism沒(méi)有bowtie index，使用 bowtie2-build 很容易自己構(gòu)建一個(gè).

image.png

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

2. GFF/GTF 格式文件

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

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

image.png

GTF(general transfer format)是GFF第二個(gè)版本嗅绰，

3 htseq-count軟件

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

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

• 模式并集， 取所有集合 S(i) 的并集橡羞。 對(duì)于大多數(shù)用例眯停，建議使用此模式。
• 模式交集卿泽，嚴(yán)格的所有集合 S(i) 的交集莺债。

• 模式交集， 非空的所有非空集 S(i) 的交集签夭。
  如果 S 恰好包含一個(gè)特征齐邦，則針對(duì)該特征計(jì)算讀取（或讀取對(duì)）第租。 如果它包含多個(gè)特征侄旬，則讀取（或讀取對(duì)）計(jì)為不明確的（不計(jì)入任何特征）煌妈，如果 S 為空儡羔，則讀刃颉（或讀取對(duì)）計(jì)為 no_feature。
  看圖更清晰理解：

  
  
  image.png

3.3.2 下載擬南芥參考基因

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

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,花費(fèi)2分鐘

bowtie2-build Arabidopsis.fasta At_ref

檢查index汰蜘，花費(fèi)幾秒鐘

bowtie2-inspect -n At-ref

3.3.4 reads映射

reads是存在以逗號(hào)隔開(kāi)的FASTQ或FASTA格式文件中

使用tophat完成

一般使用命令：

image.png

更多的選擇閱讀文檔

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

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

開(kāi)始操作

軟鏈接參考基因組的FASTA：

ln -s Arabidopsis.fasta At_ref.fa

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

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

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

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

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

每個(gè)樣本創(chuàng)一個(gè)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

此步驟花費(fèi)大約1個(gè)小時(shí)
查看任務(wù)

qstat -f

對(duì)所有的樣本進(jìn)行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ù)是否來(lái)自特定鏈的檢測(cè)（默認(rèn)值：yes）.對(duì)于stranded=no妖谴，無(wú)論是映射到與特征相同還是相反的鏈窿锉，都認(rèn)為讀取與特征重疊。 對(duì)于 stranded=yes 和單端讀取膝舅，讀取必須映射到與特征相同的鏈嗡载。 對(duì)于雙端讀取，第一次讀取必須在同一條鏈上仍稀，第二次讀取必須在相反的鏈上洼滚。 對(duì)于stranded=reverse，這些規(guī)則是相反的琳轿。

reads計(jì)數(shù)模板

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

運(yùn)行花費(fèi)半個(gè)小時(shí)判沟。

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

shell 命令

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

. 組裝讀取計(jì)數(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

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

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

3.4.1 工具介紹

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

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

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

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

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

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

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

• 進(jìn)行差異表達(dá)分析残制。Trinity 直接支持多種 DE 分析方法，包括 edgeR掖疮、DESeq2闰蚕、Limma/Voom 和 ROTS颗管。

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

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

分析使用 perl 腳本嵌入了一組流行工具進(jìn)行分析排拷。所以我們將使用一些代碼來(lái)映射猎唁、對(duì)齊和計(jì)數(shù)reads.
首先使用：align_and_estimate_abundance.pl误证， 功能是index參考轉(zhuǎn)錄組，對(duì)齊和計(jì)數(shù)reads盖腿。其中index和對(duì)其依靠的是Bowtie2軟件爽待，估計(jì)轉(zhuǎn)錄的豐富度使用RSEM。從 RNA-Seq 數(shù)據(jù)進(jìn)行轉(zhuǎn)錄量化的一個(gè)關(guān)鍵挑戰(zhàn)是處理映射到多個(gè)基因或同種型的讀數(shù)翩腐。 這個(gè)問(wèn)題對(duì)于在沒(méi)有測(cè)序基因組的情況下用從頭轉(zhuǎn)錄組組裝進(jìn)行量化特別重要鸟款，因?yàn)楹茈y確定哪些轉(zhuǎn)錄本是同一基因的同種型。 RSEM 允許從單端或雙端 RNA-Seq 數(shù)據(jù)中量化基因和同種型豐度茂卦，無(wú)需測(cè)序基因組何什。

請(qǐng)注意，Trinity提供了依賴對(duì)齊和無(wú)對(duì)齊的讀取計(jì)數(shù)兩種方案等龙。 它提供基因和轉(zhuǎn)錄水平的reads計(jì)數(shù)估計(jì)处渣。

3 數(shù)據(jù)標(biāo)準(zhǔn)化
原始讀取計(jì)數(shù)必須針對(duì)許多偏差（例如基因長(zhǎng)度或每個(gè)樣本的讀取數(shù)）進(jìn)行校正。 Trinity 提供的任何豐度估計(jì)方法都將提供源自每個(gè)轉(zhuǎn)錄本的 RNA-Seq 片段計(jì)數(shù)的轉(zhuǎn)錄本水平估計(jì)而咆，此外還提供了轉(zhuǎn)錄本表達(dá)的標(biāo)準(zhǔn)化測(cè)量霍比，該測(cè)量值考慮了轉(zhuǎn)錄本長(zhǎng)度幕袱、映射的讀取數(shù) 到轉(zhuǎn)錄本暴备，以及映射到任何轉(zhuǎn)錄本的讀取總數(shù)。 標(biāo)準(zhǔn)化表達(dá)指標(biāo)可以報(bào)告為每千堿基轉(zhuǎn)錄本長(zhǎng)度的片段映射 (FPKM) 或每百萬(wàn)轉(zhuǎn)錄本 (TPM) 的轉(zhuǎn)錄本们豌。

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

來(lái)自Araport, 需要登錄進(jìn)行免費(fèi)的注冊(cè)涯捻。再使用以下代碼獲取浅妆。

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，可以對(duì)所有樣本一次完成障癌。

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

此過(guò)程花費(fèi)大約5分鐘凌外，會(huì)生成14個(gè)文件，具有.bowtie2.和.RSEM文件

3.4.4 對(duì)read進(jìn)行映射和計(jì)數(shù)

使用ltrinity的perl命令：align_and_estimate_abundance.pl涛浙，并且使用RSEM估計(jì)方法

image.png

image.png

2建立gene_trans_map
我們需要準(zhǔn)備一個(gè)gene和轉(zhuǎn)錄名字對(duì)應(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進(jìn)行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)建多個(gè)樣本的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

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

image.png

image.png

3.4.5 生成表達(dá)矩陣

使用：trinity下的abundance_estimates_to_matrix.pl命令
該腳本將非常簡(jiǎn)單地創(chuàng)建一個(gè)矩陣轿亮，將所有樣本的讀取計(jì)數(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分鐘

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

第三部分： 3.5 鑒定差異表達(dá)的基因

使用R 包DESeq2但骨。

3.5.1 包介紹

詳細(xì)文檔介紹：https://bioconductor.org/packages/release/bioc/html/DESeq2.html.
DESeq2 允許估計(jì)來(lái)自高通量測(cè)序分析的計(jì)數(shù)數(shù)據(jù)（依賴于方差-均值）励七，并基于使用負(fù)二項(xiàng)分布的廣義線性模型 (GLM) 測(cè)試差異表達(dá)。 一個(gè)非潮疾基本的 GLM 模型如下：

image.png

DESeq2 將首先對(duì)數(shù)據(jù)進(jìn)行建模以估計(jì)模型的系數(shù)掠抬。 一旦設(shè)置了系數(shù)，就可以為每個(gè) x 確定 y添坊。
這里只是簡(jiǎn)要介紹使用的函數(shù)剿另，其他的請(qǐng)看上文的連接文檔
1 DESeqDataSetFromMatrix
DESeqDataSet 是 RangedSummarizedExperiment 的子函數(shù)，用于存儲(chǔ)輸入值贬蛙、中間計(jì)算和差異表達(dá)式分析的結(jié)果雨女。 DESeqDataSet 在“計(jì)數(shù)”矩陣中強(qiáng)制執(zhí)行非負(fù)整數(shù)值，作為分析列表中的第一個(gè)元素存儲(chǔ)阳准。 此外氛堕，必須提供指定實(shí)驗(yàn)設(shè)計(jì)的公式。
使用：DESeqDataSetFromMatrix(countData, colData, design)
其中design:設(shè)計(jì)一個(gè)公式來(lái)表達(dá)每個(gè)基因的計(jì)數(shù)如何取決于 colData 中的變量野蝇。 許多 R 公式是有效的讼稚，包括具有多個(gè)變量的設(shè)計(jì)，例如 ~ 組 + 條件绕沈，以及具有交互作用的設(shè)計(jì)锐想，例如 ~ 基因型 + 治療 + 基因型：治療。 查看各種設(shè)計(jì)的結(jié)果以及如何提取結(jié)果表乍狐。
countData,為輸入的非負(fù)整數(shù)矩陣
colData: 輸入DataFrame或者data.frame格式赠摇，其只是有一列，行是對(duì)應(yīng)著countData。
2 DESeq
DESeq 基于負(fù)二項(xiàng)分布進(jìn)行差異表達(dá)分析藕帜。 它通過(guò)以下步驟執(zhí)行默認(rèn)分析：
? 估計(jì)大小因素：estimateSizeFactors
? 估計(jì)色散：estimateDispersions
? 負(fù)二項(xiàng)式 GLM 擬合和 Wald 統(tǒng)計(jì)：nbinomWaldTest

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

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

3 results
結(jié)果從 DESeq 分析中提取結(jié)果表，給出樣本的基本均值荒适、log2 倍數(shù)變化长捧、標(biāo)準(zhǔn)誤差、檢驗(yàn)統(tǒng)計(jì)量吻贿、p 值和調(diào)整后的 p 值串结。

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

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

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

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

3.5.2 下載DESeq2

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

3.5.3 鑒定差異表達(dá)基因(成對(duì)比較)

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

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ù)挖掘

我們將對(duì) DGE 數(shù)據(jù)進(jìn)行非常簡(jiǎn)單和基本的數(shù)據(jù)挖掘。 我們將使用 Thalemine 門(mén)戶珍剑，它實(shí)現(xiàn)了一個(gè) Intermine 接口掸宛，允許非常快速地獲取有關(guān)基因列表的功能信息招拙，包括 GO 富集分析唧瘾。連接：
https://bar.utoronto.ca/thalemine/

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

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注釋基因