本文是由比利时列日大学Marc HANIKENNE课程整理。陆陆续续花费1个月完成，根据需要做的主要内容分为四个部分。
第一部分：将RNA-seq数据映射到参考基因组上
第二部分：将RNA-seq数据映射到参考转录组上，并且生成基因表达矩阵，用于第三部分分析
第三部分：使用DESeq包鉴定差异表达基因（成对），
第四部分：对差异基因进行后续的GO和KEGG注释

1 目标

RNA-Seq 模块的目标是说明如何处理和分析 RNA-Seq 数据以识别差异表达基因 (DGE)。
练习中使用真实数据集，来自暴露于两种生长条件的拟南芥的两种基因型的 Illumina RNA 测序。
需要做：
1). 在参考基因组（工具：tophat 和 htseq-count）或参考转录组（工具trinity）上映射reads，作为reads映射和计数的两种相互替代策略；
2). 使用 DESeq2包（ R 语言）鉴定不同表达的基因（DGEs）
3). 简单进行数据挖掘(GO和KEGG注释)。

2 数据介绍

模式植物拟南芥的两种基因型（wt 和突变体）在对照 (c) 和处理 (t) 条件下生长。 在实验结束时，分别收获根 (r) 和芽 (s) 组织（每个样品 4 株植物）。 本实验独立重复 3 次。 使用 NextSeq Illumina 测序仪以高通量模式（~400 Mio. clusters）在两次运行中制备和测序总 RNA，以获得 75 bp 的单端读数。 注意，已经使用 fastqc 和trimmomatic 分析和质控了原始reads（这两个软件我会在基因组组装分析中介绍）。

3 分析数据

3.1 查看数据

head <your_sample>.fastq
查看每个文件数据的reads数：

image.png

将所有的样本名字放在一个文件（samples.ids）中, 这样方便后续进行批量处理。
shell code

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

3.2 RNA-Seq数据分析中reads映射的一般考虑

在分析 RNA-Seq 数据时，有两种通用策略可以在计数之前映射reads。
(i) Reads可以使用splice aware比对工具映射到参考基因组上：比对器将允许拆分读数以跨越内含子。当感兴趣的物种有参考基因组可用时，这种策略通常受到青睐，因为它允许检测以前没有描述过的转录本或转录本变异。然而，这种方法仅限于（相对）少数有高质量参考基因组可用的物种。
(ii) Reads也可以映射到参考转录组上。比对更容易，并且当目标物种不存在参考转录组时，可以使用 RNA-Seq 读取组装自定义参考转录组，使其几乎适用于所有物种。然而，这种方法仅限于检测已知转录本和转录本变体的表达。

第一部分：3.3 Read映射到参考基因组

3.3.1 工具介绍

1. tophat 软件

我们将使用流行的 tophat，这是一个将 RNA-Seq 读数与基因组对齐以识别外显子-外显子剪接连接的程序。 它建立在超快的短读映射程序 bowtie2 之上。 该软件针对 75bp 或更长的读取进行了优化。
更多介绍点击查看：Tophat 链接。

TopHat 如何找到连接点的原理：

TopHat 可以在没有参考注释的情况下找到拼接点。通过首先将 RNA-Seq的reads映射到基因组，TopHat可以识别出潜在的外显子，这是因为许多 RNA-Seq reads将与基因组连续对齐。使用这个初始映射信息，TopHat 建立一个可能的剪接连接的数据库，然后将读取映射到这些连接以确认它们。

短read测序机目前可以产生 100 bp 或更长的reads，但许多外显子比这短，因此它们会在初始映射中被遗漏。TopHat 主要通过将所有输入reads分成更小的段来解决这个问题，然后独立映射。这些小段对齐在程序的最后一步重新组合在一起，以产生端到端的read对齐。

TopHat 从两个证据来源生成可能的剪接点数据库。第一个也是最有力的证据: 剪接点来源是当来自同一reads的两个片段（对于至少 45bp 的reads）在同一基因组序列上的某个距离处作图时，或者当一个内部片段未能作图时——再次表明这种读取跨越多个外显子。使用这种方法，“GT-AG”、“GC-AG”和“AT-AC”内含子将从头开始找到。第二个来源是“coverage islands,”的配对，它们是初始映射中堆积读取的不同区域。转录组中的相邻岛通常拼接在一起，因此 TopHat 寻找将这些与内含子连接起来的方法。我们只建议用户将第二个选项 (--coverage-search) 用于短读取 (< 45bp) 和少量读取 (<= 1000 万)。后一个选项只会报告“GT-AG”内含子之间的比对.

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

想要使用这个软件，首先需要使用一下命令：

image.png
Bowtie2 用于对齐基因组上的reads。

Bowtie 2是一种超快且内存高效的工具，用于将测序读数与长参考序列对齐。它特别擅长将大约 50 到 100 个字符的读数与相对较长的（例如哺乳动物）基因组进行比对。Bowtie 2 使用FM 索引（基于Burrows-Wheeler 变换或BWT）对基因组进行索引，以保持其较小的内存占用：对于人类基因组，其内存占用通常约为 3.2 GB 的 RAM。Bowtie 2 支持间隙、局部和双端对齐模式。可以同时使用多个处理器以实现更高的对齐速度。

Bowtie 2 以SAM格式输出对齐方式，从而能够与大量使用同样使用SAM文件格式的其他工具（例如SAMtools、GATK）进行互操作。Bowtie 2 在GPLv3 许可下分发，它在 Windows、Mac OS X 和 Linux 和 BSD 下的命令行上运行。

Bowtie 2通常是比较基因组学管道的第一步，包括变异调用、ChIP-seq、RNA-seq、BS-seq。Bowtie 2和Bowtie（此处也称为“ Bowtie 1 ”）也紧密集成到许多其他工具中，此处列出了其中一些。

要找到与 Tophat 的连接点，您首先需要为 RNA-Seq 实验中的organism安装bowtie index。 bowtie 站点为人类、小鼠、果蝇等提供了预先构建的bowtie index。 如果你的organism没有bowtie index，使用 bowtie2-build 很容易自己构建一个.

image.png

Bowtie2-inspect 从 bowtie index 中提取信息，允许确定它是什么类型的索引以及使用什么参考序列来构建它。

2. GFF/GTF 格式文件

通过提供描述基因组特征（例如外显子/内含子）染色体位置的基因组注释文件，可以帮助通过 tophat 在参考基因组上进行reads映射。 注释文件以 GFF/GTF 格式提供。

Tophat使用的的genome annotation文件就是GFF/GTF格式。
GFF一般就9列，是以Tab间隔，对应名字：

image.png

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

3 htseq-count软件

给定一个具有对齐测序读数的文件和基因组特征列表，htseq-count会计算有多少reads映射到每个特征。这里的特征是染色体上的区间（即，位置范围）或这些区间的联合。在 RNA-Seq 的情况下，特征通常是基因，其中每个基因在这里被视为其所有外显子的联合。人们也可以将每个外显子视为一个特征，例如，为了检查可变剪接。对于比较 ChIP-Seq，特征可能是预定列表中的结合区域。

必须特别注意决定如何处理重叠多个特征的reads。 htseq-count 脚本允许在三种模式之间进行选择。 htseq-count 的三种重叠解析模式的工作原理如下：对于 read 中的每个位置 i，定义一个集合 S(i) 为所有与位置 i 重叠的特征的集合。 然后，考虑集合 S，它是（i 遍历读取或读取对中的所有位置）

• 模式并集， 取所有集合 S(i) 的并集。 对于大多数用例，建议使用此模式。
• 模式交集，严格的所有集合 S(i) 的交集。

• 模式交集， 非空的所有非空集 S(i) 的交集。
  如果 S 恰好包含一个特征，则针对该特征计算读取（或读取对）。 如果它包含多个特征，则读取（或读取对）计为不明确的（不计入任何特征），如果 S 为空，则读取（或读取对）计为 no_feature。
  看图更清晰理解：

  
  image.png

3.3.2 下载拟南芥参考基因

网站：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输出结果的文目录
--min-intron-length <int> ##内含子的最短长度：默认70
--max-intron-length <int>##内含子的最长长度：默认500000
-G --GTF <GTF/GFF3 file> # 为 TopHat 提供一组基因模型注释和/或已知转录本，作为 GTF 2.2 或 GFF3 格式的文件。 如果提供此选项，TopHat 将首先提取转录本序列，并首先使用 Bowtie 将读数与该虚拟转录组对齐。 只有没有完全映射到转录组的读数才会被映射到基因组上。 在转录组上进行映射的读数将被转换为基因组映射（根据需要拼接）并与最终 tophat 输出中的新映射和连接点合并。

请注意，所提供的 GTF/GFF 文件的第一列（指示特征所在的染色体或重叠群的列）中的值必须与用于 TopHat 的 Bowtie 索引中的参考序列名称相匹配。 你可以使用 bowtie-inspect 进行检查。
--transcriptome-index <dir/prefix> ##当为 TopHat 提供已知的转录文件（上面的 -G/--GTF 选项）时，会构建转录组序列文件，并且必须为其创建 Bowtie index，以便将读取与已知转录本对齐。创建此 Bowtie 索引可能非常耗时，并且在许多情况下，相同的转录组数据被用于将多个样本与 TopHat 对齐。因此，转录组索引和相关数据文件（原始 GFF 文件）可以在使用此选项的多次 TopHat 运行中重复使用，因此这些文件仅针对具有给定转录本数据的第一次运行创建。如果计划使用相同的转录组数据进行多次 TopHat 运行，则应首先使用 -G/--GTF 选项以及指向目录和名称前缀的 --transcriptome-index 选项运行 TopHat，该选项将指示转录组数据的位置文件将被存储。然后使用相同的 --transcriptome-index 选项值的后续 TopHat 运行将直接使用在第一次运行中创建的转录组数据（第一次运行后不需要 -G 选项）。

开始操作

软链接参考基因组的FASTA：

ln -s Arabidopsis.fasta At_ref.fa

创建转录组的index.只需要做一次，即可用于全部样本，耗时5分钟

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

会在transcriptome_data/ 下产生10个文件

映射reads， 先创建一个模板

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

每个样本创一个sh

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

提交任务：

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

此步骤花费大约1个小时
查看任务

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
脚本输出一个表，其中包含每个feature(这里指基因)的计数，然后是特殊计数器，用于计算由于各种原因未针对任何feature进行计数的reads。 特殊计数器的名称都以双下划线开头，以便于过滤。特殊计数器是：
image.png
重要提示：strandedness的默认值为 yes。 如果你的 RNA-Seq 数据不是用特定链的协议制作的，这会导致一半的读数丢失。 因此，除非您有特定于链的数据，否则请确保设置选项 –stranded=no！
htseq-count具有很多opition，请查看链接文档
一些opition:
-f <sam or bam># 输入文件，sam或者bam格式

-s <yes/no/reverse>
数据是否来自特定链的检测（默认值：yes）.对于stranded=no，无论是映射到与特征相同还是相反的链，都认为读取与特征重叠。 对于 stranded=yes 和单端读取，读取必须映射到与特征相同的链。 对于双端读取，第一次读取必须在同一条链上，第二次读取必须在相反的链上。 对于stranded=reverse，这些规则是相反的。

reads计数模板

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

运行花费半个小时。

检索对齐统计信息

shell 命令

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

. 组装读取计数矩阵

得到基因名字

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

得到的这个结果文件，将用于后续的统计分析，发现DGEs

第二部分： 3.4 read映射到参考转录组

3.4.1 工具介绍

1. trinity 软件
   由Broad 研究所和耶路撒冷希伯来大学开发的 Trinity代表了一种从 RNA-seq 数据高效、稳健地从头重建转录组的新方法。Trinity 结合了三个独立的软件模块：Inchworm、Chrysalis 和 Butterfly，依次应用于处理大量 RNA-seq 读数。Trinity 将序列数据划分为许多单独的 de Bruijn 图，每个图代表给定基因或基因座的转录复杂性，然后独立处理每个图以提取全长剪接异构体并梳理源自旁系同源基因的转录本。简而言之，过程是这样工作的：

• Inchworm将 RNA-seq 数据组装成独特的转录本序列，通常为显性同种型生成全长转录本，但随后仅报告选择性剪接转录本的独特部分。

• Chrysalis将 Inchworm contigs 聚类成簇，并为每个簇构建完整的 de Bruijn 图。每个簇代表给定基因（或共享序列的基因组）的完整转录复杂性。然后 Chrysalis 在这些不相交的图之间划分完整的读取集。

• Butterfly然后并行处理各个图，跟踪图中读取和读取对的路径，最终报告可变剪接同种型的全长转录本，并梳理出对应于旁系同源基因的转录本。

2 转录组组装下游分析
组装完成后，您可能需要进行多项分析，以根据组装的转录本和输入的 RNA-Seq 数据探索生物体生物学的各个方面。此类研究可能包括：

• 量化转录本和基因丰度。这是许多其他分析的先决条件，例如检查样本中差异表达的转录本。

• 质量检查您的样品和生物复制品。确保您的样本和生物学重复之间的关系有意义。如果数据存在任何混杂因素，例如异常值或批处理效应，您将希望尽早发现它们并在进一步的数据探索中考虑到它们。

• 进行差异表达分析。Trinity 直接支持多种 DE 分析方法，包括 edgeR、DESeq2、Limma/Voom 和 ROTS。

• 使用提取的编码区TransDecoder和功能注释使用的成绩单Trinotate。

• 如果您的生物体具有组装的基因组，请考虑使用 Trinity 转录组组装使用PASA进行基因结构注释。

分析使用 perl 脚本嵌入了一组流行工具进行分析。所以我们将使用一些代码来映射、对齐和计数reads.
首先使用：align_and_estimate_abundance.pl， 功能是index参考转录组，对齐和计数reads。其中index和对其依靠的是Bowtie2软件，估计转录的丰富度使用RSEM。从 RNA-Seq 数据进行转录量化的一个关键挑战是处理映射到多个基因或同种型的读数。 这个问题对于在没有测序基因组的情况下用从头转录组组装进行量化特别重要，因为很难确定哪些转录本是同一基因的同种型。 RSEM 允许从单端或双端 RNA-Seq 数据中量化基因和同种型丰度，无需测序基因组。

请注意，Trinity提供了依赖对齐和无对齐的读取计数两种方案。 它提供基因和转录水平的reads计数估计。

3 数据标准化
原始读取计数必须针对许多偏差（例如基因长度或每个样本的读取数）进行校正。 Trinity 提供的任何丰度估计方法都将提供源自每个转录本的 RNA-Seq 片段计数的转录本水平估计，此外还提供了转录本表达的标准化测量，该测量值考虑了转录本长度、映射的读取数 到转录本，以及映射到任何转录本的读取总数。 标准化表达指标可以报告为每千碱基转录本长度的片段映射 (FPKM) 或每百万转录本 (TPM) 的转录本。

3.4.2 拟南芥参考转录组

来自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拟南芥参考转录组

使用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进行映射和计数

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

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

2建立gene_trans_map
我们需要准备一个gene和转录名字对应的文件，并且需要将转录组的fasta文件排序成参考转录组的方式
从fasta文件中提取转录的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 %]

创建多个样本的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

批量提交任务：

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

这一步大概花费90分钟
再查看你的结果：

image.png
image.png

3.4.5 生成表达矩阵

使用：trinity下的abundance_estimates_to_matrix.pl命令
该脚本将非常简单地创建一个矩阵，将所有样本的读取计数数据组合在一起。
使用：

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 文件提供未经标准化的原始读取计数。 这是必须用于使用 DESeq2 进行进一步 DGE 分析的文件。
额外的 .matrix 文件对应于 TPM 表达值矩阵（未跨样本归一化）和 TMM 归一化表达值矩阵（应用了跨样本归一化）。 有关此查看的更多详细信息：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 允许估计来自高通量测序分析的计数数据（依赖于方差-均值），并基于使用负二项分布的广义线性模型 (GLM) 测试差异表达。 一个非常基本的 GLM 模型如下：

image.png

DESeq2 将首先对数据进行建模以估计模型的系数。 一旦设置了系数，就可以为每个 x 确定 y。
这里只是简要介绍使用的函数，其他的请看上文的连接文档
1 DESeqDataSetFromMatrix
DESeqDataSet 是 RangedSummarizedExperiment 的子函数，用于存储输入值、中间计算和差异表达式分析的结果。 DESeqDataSet 在“计数”矩阵中强制执行非负整数值，作为分析列表中的第一个元素存储。 此外，必须提供指定实验设计的公式。
使用：DESeqDataSetFromMatrix(countData, colData, design)
其中design:设计一个公式来表达每个基因的计数如何取决于 colData 中的变量。 许多 R 公式是有效的，包括具有多个变量的设计，例如 ~ 组 + 条件，以及具有交互作用的设计，例如 ~ 基因型 + 治疗 + 基因型：治疗。 查看各种设计的结果以及如何提取结果表。
countData,为输入的非负整数矩阵
colData: 输入DataFrame或者data.frame格式，其只是有一列，行是对应着countData。
2 DESeq
DESeq 基于负二项分布进行差异表达分析。 它通过以下步骤执行默认分析：
• 估计大小因素：estimateSizeFactors
• 估计色散：estimateDispersions
• 负二项式 GLM 拟合和 Wald 统计：nbinomWaldTest

有关每个步骤的完整详细信息，请参阅相应功能的手册页。 DESeq 函数返回 DESeqDataSet 对象后，可以使用结果函数生成结果表（log2 倍数变化和 p 值）。 有关针对多重测试校正的独立过滤和 p 值调整的信息，请参阅结果手册页。

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

3 results
结果从 DESeq 分析中提取结果表，给出样本的基本均值、log2 倍数变化、标准误差、检验统计量、p 值和调整后的 p 值。

resultsNames 返回模型的估计效应（系数）的名称。
使用

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

object是 DESeqDataSet，已在其上调用以下函数之一：DESeq、nbinocontrastmWaldTest 或 nbinomLRT。
contrast参数指定从对象中提取什么比较以构建结果表。
在我们的例子中，一个字符向量正好包含三个元素：设计公式中一个因子的名称、折叠变化的分子级别的名称以及折叠变化的分母级别的名称（最简单的情况）。
lfcThreshold 是一个非负值，指定 log2 倍变化阈值。 默认值为 0，对应 log2 倍数变化为零的测试。

alpha 用于优化独立过滤的显着性截止值（默认为 0.1）。 如果调整后的 p 值截止值 (FDR) 不是 0.1，则 alpha 应设置为该值。
plotCounts
plotCounts 允许在对数尺度上为单个基因创建归一化计数图。
使用：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 鉴定差异表达基因(成对比较)

我们将在下面发现如何编写允许这样做的 R 脚本。 您需要在 RStudio 内的同一个 R 脚本中按顺序添加每个新步骤。 我们将首先基于独立于治疗的基因型（WT vs 突变体）识别显示 DGE 的基因，然后基于独立于基因型的治疗（Ctrl vs Treat）的 DGE，最后看看治疗对 每个基因型。 这些分析中的每一个都对应于统计测试的不同设计，在我们的脚本中必须考虑到这一点。

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 数据挖掘

我们将对 DGE 数据进行非常简单和基本的数据挖掘。 我们将使用 Thalemine 门户，它实现了一个 Intermine 接口，允许非常快速地获取有关基因列表的功能信息，包括 GO 富集分析。连接：
https://bar.utoronto.ca/thalemine/

为了使用这个工具，我们首先需要从DESeq2生成的基因列表中提取基因名称
DESeq2 分析生成了 24 个文件（8 个成对.txt，8 个高.txt 和 8 个低.txt）。 由于它们对应于我们过滤后的数据，我们只对 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注释基因