PlantNGSTools使用说明

1 软件包安装

1.1 安装PlantNGSTools

devtools::install_github('biomarble/PlantNGSTools',upgrade=F)

如果GitHub访问有问题，可以使用下面的命令安装

devtools::install_url('https://archive.fastgit.org/biomarble/PlantNGSTools/archive/main.zip',upgrade = F)

1.2 修复EBSeq bug

由于EBSeq的依赖包 blockmodeling 作者名字包含特殊字符，在中文的windows下加载会报错 因此，使用如下方式修复此bug:

file.copy(system.file('replacement','CITATION',package='PlantNGSTools'),find.package('blockmodeling'),overwrite=T)

2 加载软件包

library(PlantNGSTools)

本软件使用到的示例文件，可用如下命令拷贝到当前目录：

file.copy(system.file('demo','',package = 'PlantNGSTools'),'.',overwrite = T,recursive = T)

3 差异分析

3.1 输入文件读取

3.1.1 基因表达量Count文件

差异分析软件（DESeq2和EBSeq）的输入文件并非标准化后的FPKM或TPM矩阵文件，而是未标准化前的read count矩阵文件。

其中，每行为一个基因，每列为一个样本：

count = read.delim("demo/count.txt", header = T, sep = "\t", row.names = 1)
S1 S2 S3 S4 S5 S6 S7 S8 S9
Os01g0100200 1 0 0 0 1 0 2 1 0
Os01g0100400 344 330 415 345 306 477 435 488 284
Os01g0100466 0 0 0 0 0 0 0 0 0
Os01g0100500 1281 1213 1406 912 1025 1329 1111 1336 769
Os01g0100600 101 108 118 84 80 106 93 129 78
Os01g0100900 1157 1025 1066 787 1051 1051 951 1182 613

3.1.2 样本分组文件

对于有生物学重复样本的情况，必须输入样本分组信息文件， 其中包含了每个样本的分组所属，同一个分组内的样本互为生物学重复。 差异表达分析时，会直接使用组名指定差异分组。

group = read.delim("demo/group.txt", header = T, sep = "\t")
Sample Group
S1 A
S2 A
S3 A
S4 B
S5 B
S6 B
S7 C
S8 C
S9 C

3.2 差异表达分析

3.2.1 有生物学重复-DESeq2

对于有生物学重复的情况，使用DESeq2软件分析。 参数中control和treat用于定义对照组和处理组的组名，如果不输入此两个参数则进入交互模式。

resultFilePath = DEGAnalysis_DESeq2(count, group, outdir = "output", useFDR = F, cut = 0.05, FCcut = 2, control = "A", treat = "B")
差异基因分析结果见以下目录:
output/A.vs.B.AllGene.tsv 
output/A.vs.B.DEG.tsv 
output/A.vs.B.DEGlist.tsv

差异分析结果，保存在outdir定义的文件夹下

文件名 含义
A.vs.B.AllGene.tsv 所有基因的log2FC，Pvalue等信息文件
A.vs.B.DEG.tsv 差异基因文件，包含基因的log2FC，Pvalue等信息
A.vs.B.DEGlist.tsv 差异基因列表，仅有差异基因名称

其中文件名中.vs.前面的A为control组，vs后面的B为case组，\(log_2FoldChange\)的计算方法为： \[ log_2(Fold\ Change) = log_2\frac{Expression(Case)}{Expression(Control)} \]

DEGAnalysis_DESeq2差异分析结果resultFilePath中包含了所有基因文件resultFilePath$allgene，以及差异基因文件resultFilePath$deg。

读取差异分析结果如下：

allgenes = read.delim(resultFilePath$allgene, header = T, row.names = 1)
degs = read.delim(resultFilePath$deg, header = T, row.names = 1)
S1 S2 S3 S4 S5 S6 MeanExpression log2FoldChange Pvalue FDR
Os01g0111700 78.796912 12.244902 19.766978 6.259099 16.000699 5.3400162 23.06810 -2.018871 0.0142885 0.1371295
Os01g0112400 259.296814 13.186818 79.067911 594.614457 3360.146851 180.6705496 747.83057 3.556148 0.0117534 0.1230981
Os01g0112600 92.540559 60.282594 21.485845 11.266379 1.230823 0.8900027 31.28270 -3.755901 0.0000148 0.0010328
Os01g0117500 5.497459 1.883831 7.734904 17.525479 24.616460 15.1300460 12.06470 1.906998 0.0107076 0.1168868
Os01g0121600 1401.852035 815.698851 184.778270 386.812352 158.776170 131.7204007 513.27301 -1.827102 0.0363427 0.2373220
Os01g0123900 50.393374 128.100512 31.799051 27.540038 13.539053 31.1500948 47.08702 -1.532813 0.0091684 0.1070213

3.2.2 无生物学重复-EBSeq

对于无生物学重复的情况，可以使用EBSeq进行差异分析。

指定control的样本名和treat的样本名，如果不输入此两个参数则进入交互模式。

resultFilePath = DEGAnalysis_EBSeq(count, outdir = "output", cut = 0.05, FCcut = 2, control = "S1", treat = "S2")
Removing transcripts with 100 th quantile < = 0 
14516 transcripts will be tested
差异基因分析结果见以下目录:
output/S1.vs.S2.AllGene.tsv 
output/S1.vs.S2.DEG.tsv 
output/S1.vs.S2.DEGlist.tsv 
allgenes = read.delim(resultFilePath$allgene, header = T, row.names = 1)
degs = read.delim(resultFilePath$deg, header = T, row.names = 1)
S1 S2 MeanExpression log2FoldChange FDR
Os01g0111700 84.69305 13.200611 48.94683 -2.632023 0.0000000
Os01g0112400 278.69922 14.216043 146.45763 -4.241132 0.0000000
Os01g0123900 54.16416 138.098700 96.13143 1.341436 0.0001719
Os01g0127600 1952.86414 7937.628966 4945.24656 2.022810 0.0000000
Os01g0614300 315.13693 136.067837 225.60238 -1.208344 0.0000004
Os01g0677400 24.62007 2.030863 13.32547 -3.285424 0.0056083

3.3 差异基因分析绘图

3.3.1 火山图

VolcanoPlot(allgenes, useFDR = T, cut = 0.05, FCcut = 2, MainTitle = "Title", showlabel = "topPvalue", showlabel.num = 10, xlim = 10)

选项 含义 可选范围
showlabel 是否显示显著基因的名称 “allSig”：标注所有显著基因;
‘topPvalue’：标注Pvalue最小的基因;
‘topFC’：标注FC最大的基因;
‘no’:不标注
showlabel.num 指定显示的基因数目
当showlabel为topPvalue或topFC时有效		 整数数字
xlim x轴显示范围 整数数字
useFDR 是否使用FDR作为显著性
具体显著性阈值由cut和FCcut来定义		 T或F
cut 显著性阈值，当useFDR=T时表示FDR阈值
当useFDR=F时表示Pvalue阈值		 0~1的小数
FCcut FoldChange阈值 大于1的任意数字
MainTitle 图片的标题 字符串

3.3.2 MA图

MAPlot(allgenes, useFDR = T, cut = 0.05, FCcut = 2, MainTitle = "Title", showlabel = "topFC", showlabel.num = 10, ylim = 10)

选项 含义 可选范围
showlabel 是否显示显著基因的名称 “allSig”：标注所有显著基因;
‘topPvalue’：标注Pvalue最小的基因;
‘topFC’：标注FC最大的基因;
‘no’:不标注
showlabel.num 指定显示的基因数目
当showlabel为topPvalue或topFC时有效		 整数数字
xlim x轴显示范围 整数数字
useFDR 是否使用FDR作为显著性
具体显著性阈值由cut和FCcut来定义		 T或F
cut 显著性阈值，当useFDR=T时表示FDR阈值
当useFDR=F时表示Pvalue阈值		 0~1的小数
FCcut FoldChange阈值 大于1的任意数字
MainTitle 图片的标题 字符串

4 富集分析

4.1 常见物种富集分析

4.1.1 KEGG Pathway富集分析

4.1.1.1 查看KEGG数据库中的物种信息

本软件包中已包含了常用物种的KEGGpathway信息，可直接通过参数指定调用，可以通过如下命令查看目前包含了哪些物种：

KEGGdbInfo()
  taxon      ReferenceGenomeVersion   update           IDformat
    osa                   IRGSP-1.0 20210518       Os05g0375100
    ath                      TAIR10 20210518          AT1G01090
    mdm                 ASM211411v1 20210518       MD17G0149700
    pop                  Pop_tri_v3 20210518 POPTR_001G122400v3
    sbi      Sorghum_bicolor_NCBIv3 20210518  SORBI_3008G059900
    sly                       SL3.0 20210518   Solyc09g064370.3
 zma.V4               B73 RefGen_v4 20210518     Zm00001d042869
 zma.V5    Zm-B73-REFERENCE-NAM-5.0 20210518    Zm00001eb148000
    gmx            Glycine_max_v2.1 20210613    GLYMA_19G000700
    pca Pogostemon_Cablin_sdata2018 20210615         Pcab093315

4.1.1.2 KEGG pathway富集分析

deg = read.delim("demo/demo.rice.deg.txt", header = F)$V1
KEGGresult = KEGGenrich(deglist = deg, useFDR = F, taxon = "osa", outdir = "output", outprefix = "demo.osa")
Enrichment Done! Result files saved:
output/demo.osa.tsv
output/demo.osa.html
选项 含义 可选范围
deglist 差异基因列表变量 必填
taxon 物种简称 KEGGdbInfo()结果中taxon列
useFDR 是否使用FDR作为显著性
具体显著性阈值由cut来定义		 T或F
outdir 输出结果目录 必填
outprefix 输出结果文件前缀 必填

4.1.1.3 KEGG富集结果

head(KEGGresult)
ID Pathway DEGsInPathway GenesInPathway Pvalue
21025 osa01110 Biosynthesis of secondary metabolites 33 950 0.0015104
14827 osa01100 Metabolic pathways 47 1658 0.0076116
13605 osa00940 Phenylpropanoid biosynthesis 9 162 0.0077934
3002 osa00073 Cutin, suberine and wax biosynthesis 3 20 0.0086401
13301 osa00906 Carotenoid biosynthesis 3 23 0.0128104
8632 osa00520 Amino sugar and nucleotide sugar metabolism 7 121 0.0151567

网页输出文件除了包含以上信息之外，具有如下额外内容：
1. pathway中高亮标出差异基因
2. DEGsInPathway和GenesInPathway对应的详细基因列表

示例网页文件链接

4.1.1.4 KEGG富集气泡图作图

KEGGbubble(dataset = KEGGresult, top = 10, MainTitle = "KEGG pathway Enrichment plot", xangle = 45, TrimName = F, NameWidth = 50, ColorScheme = "rw", ColorReverse = F)

选项 含义 可选范围
dataset KEGGenrich获得的注释结果变量 必填
top 输出BP、MF、CC分别的前n个功能节点
最终输出的功能同时满足top与cut参数		 >1即可
MainTitle 图片标题 任意字符串
xangle X轴文字旋转角度 30,45,90
ColorScheme 颜色方案 ryb, rwb, ryg, rwg, rw, bw, gw
其中r为红色，y为黄色，w为白色，b为蓝色，g为绿色
ColorReverse 是否反转ColorScheme的颜色的顺序 T或F
TrimName 是否对名称特别长的Pathway名称进行剪切缩短 T或F
NameWidth 如果 TrimName为T，则直接剪切前NameWidth个字符，否则按照NameWidth长度自动换行 默认为50

4.1.2 GO富集分析

4.1.2.1 查看GO数据库中的物种信息

本软件包中已包含了常用物种的GO注释信息，通过如下命令查看：

GOdbInfo()
  taxon      ReferenceGenomeVersion   update           IDformat
    osa                   IRGSP-1.0 20210518       Os05g0375100
    ath                      TAIR10 20210518          AT1G01090
    mdm                 ASM211411v1 20210518       MD17G0149700
    pop                  Pop_tri_v3 20210518 POPTR_001G122400v3
    sbi      Sorghum_bicolor_NCBIv3 20210518  SORBI_3008G059900
    sly                       SL3.0 20210518   Solyc09g064370.3
 zma.V4               B73 RefGen_v4 20210518     Zm00001d042869
 zma.V5    Zm-B73-REFERENCE-NAM-5.0 20210518    Zm00001eb148000
    pca Pogostemon_Cablin_sdata2018 20210613    GLYMA_19G000700
    gmx            Glycine_max_v2.1 20210613         Pcab093315

4.1.2.2 GO富集分析

deg = read.delim("demo/demo.rice.deg.txt", header = F)$V1
GOresult = GOEnrich(deglist = deg, taxon = "osa", useFDR = F, cut = 0.05, outdir = "output", outprefix = "demo.osa")
genome version: IRGSP-1.0 
Processing BP class.........done!
Processing MF class.........done!
Processing CC class.........done!

GO enrichment Done!
Result files saved at:
output/demo.osa.GO.significant.tsv
output/demo.osa.GO.all.tsv
选项 含义 可选范围
deglist 差异基因列表变量 必填
taxon 物种简称 GOdbInfo()结果中taxon列
useFDR 是否使用FDR作为显著性
具体显著性阈值由cut来定义		 T或F
cut 显著性阈值，当useFDR=T时表示FDR阈值
当useFDR=F时表示Pvalue阈值		 0~1的小数
outdir 输出结果目录 必填
outprefix 输出结果文件前缀 必填

4.1.2.3 GO富集气泡图作图

GObubble(dataset = GOresult, onlySig = T, useFDR = F, cut = 0.05, top = 10, TrimTerm = F, TermWidth = 50, MainTitle = "GO plot", ColorScheme = "rw", ColorReverse = F)

选项 含义 可选范围
dataset GOEnrich获得的注释结果变量 必填
onlySig 是否只显示显著的功能节点
useFDR 是否使用FDR作为显著性
具体显著性阈值由cut来定义		 T或F
cut 显著性阈值，当useFDR=T时表示FDR阈值
当useFDR=F时表示Pvalue阈值		 0~1的小数
top 输出BP、MF、CC分别的前n个功能节点
最终输出的功能同时满足top与cut参数		 >1即可
MainTitle 图片标题 任意字符串
ColorScheme 颜色方案 ryb, rwb, ryg, rwg, rw, bw, gw
其中r为红色，y为黄色，w为白色，b为蓝色，g为绿色
ColorReverse 是否反转ColorScheme的颜色的顺序 T或F
TrimTerm 是否对长名称的Term进行剪切缩短 T或F
TermWidth 如果 TrimTerm为T，则直接剪切前TermWidth个字符，否则按照TermWidth长度自动换行 默认为60

4.1.2.4 GO二级节点注释条形图

绘制GO数据库的所有二级节点的所有基因与差异基因条形图：

GOBar(deg, taxonid = "osa")

4.2 自定义注释的富集分析（无参物种，或自注释物种）

4.2.1 基于eggNOG-Mapper的GO富集

eggNOG-Mapper 可以对提交的蛋白/CDS序列进行GO注释。

获取tsv格式注释结果（下载页面中写为CSV格式），可直接输入GOEnrich_eggnog函数使用。

deg = read.delim("demo/demo.rice.deg.txt", header = F)$V1
GOresult = GOEnrich_eggnog(deglist = deg, useFDR = F, cut = 0.05, eggnogFile = "demo/demo.eggnog.tsv", outdir = "output", outprefix = "eggNOG")
Processing BP class.........done!
Processing MF class.........done!
Processing CC class.........done!

GO enrichment Done!
Result files saved at:
output/eggNOG.GO.significant.tsv
output/eggNOG.GO.all.tsv
GObubble(dataset = GOresult, onlySig = T, useFDR = F, cut = 0.05, top = 10, MainTitle = "GO plot", TrimTerm = F, TermWidth = 50, ColorScheme = "rw", ColorReverse = F)

选项 含义 可选范围
deglist 差异基因列表变量 必填
eggnogFile eggNOG-Mapper注释输出的tsv文件路径 必填
useFDR 是否使用FDR作为显著性
具体显著性阈值由cut来定义		 T或F
cut 显著性阈值，当useFDR=T时表示FDR阈值
当useFDR=F时表示Pvalue阈值		 0~1的小数
outdir 输出结果目录 必填
outprefix 输出结果文件前缀 必填

4.2.1.1 GO二级节点注释条形图

绘制GO数据库的所有二级节点的所有基因与差异基因条形图：

GOBar(deg, eggnog = "demo/demo.eggnog.tsv")

4.2.2 基于PANNZER2的GO富集

PANNZER2数据库可以对提交的序列文件进行GO注释，速度快，且注释全面。

下载GO prediction的文件（单击鼠标右键，另存为） ，示例文件链接

PANNZER2注释结果页面

deg = read.delim("demo/demo.rice.deg.txt", header = F)$V1
GOresult = GOEnrich_pannzer2(deglist = deg, useFDR = F, cut = 0.05, pannzerfile = "demo/demo.pannzer.GO.txt", outdir = "output", outprefix = "Pannzer2.go")
Processing BP class.........done!
Processing MF class.........done!
Processing CC class.........done!

GO enrichment Done!
Result files saved at:
output/Pannzer2.go.GO.significant.tsv
output/Pannzer2.go.GO.all.tsv
GObubble(dataset = GOresult, onlySig = T, useFDR = F, cut = 0.05, top = 10, TrimTerm = F, TermWidth = 50, MainTitle = "GO plot", ColorScheme = "rw", ColorReverse = F)

选项 含义 可选范围
deglist 差异基因列表变量 必填
pannzerfile PANNZER2注释输出的GO.out文件路径 必填
useFDR 是否使用FDR作为显著性
具体显著性阈值由cut来定义		 T或F
cut 显著性阈值，当useFDR=T时表示FDR阈值
当useFDR=F时表示Pvalue阈值		 0~1的小数
outdir 输出结果目录 必填
outprefix 输出结果文件前缀 必填

4.2.2.1 GO二级节点注释条形图

绘制GO数据库的所有二级节点的所有基因与差异基因条形图：

GOBar(deg, pannzer2 = "demo/demo.pannzer.GO.txt")

4.2.3 自主注释的GO富集分析

4.2.3.1 GeneID与GOID对应表: 一对一

如果已有的GO注释格式为：GeneID和GOID一对一，如下所示：

第一列为GeneID，第二列为GOID。

以txt格式存储，列与列之间以tab（制表符）分隔。

GeneID GOID
Pcab132295 GO:0005247
Pcab132295 GO:0006821
Pcab132295 GO:0016021
Pcab132295 GO:0034220
Pcab132296 GO:0005247
Pcab132296 GO:0006821
Pcab132296 GO:0016021
Pcab132296 GO:0034220
Pcab142741 GO:0004601
Pcab142741 GO:0098869

使用GOEnrich_customTable函数，tableFile定义文件路径：

deg = read.delim("demo/demo.rice.deg.txt", header = F)$V1
GOresult = GOEnrich_customTable(deglist = deg, useFDR = F, cut = 0.05, tableFile = "demo/demo.GOtable.txt", outdir = "output", outprefix = "table.GO")
Processing BP class.........done!
Processing MF class.........done!
Processing CC class.........done!

GO enrichment Done!
Result files saved at:
output/table.GO.GO.significant.tsv
output/table.GO.GO.all.tsv
GObubble(dataset = GOresult, onlySig = T, useFDR = F, cut = 0.05, top = 10, TrimTerm = F, TermWidth = 50, MainTitle = "GO plot", ColorScheme = "rw", ColorReverse = F)

选项 含义 可选范围
deglist 差异基因列表变量 必填
tableFile 一对一注释文件路径 必填
useFDR 是否使用FDR作为显著性
具体显著性阈值由cut来定义		 T或F
cut 显著性阈值，当useFDR=T时表示FDR阈值
当useFDR=F时表示Pvalue阈值		 0~1的小数
outdir 输出结果目录 必填
outprefix 输出结果文件前缀 必填

4.2.3.2 GO二级节点注释条形图

绘制GO数据库的所有二级节点的所有基因与差异基因条形图：

GOBar(deg, customTable = "demo/demo.GOtable.txt")

4.2.3.3 GeneID与GOID对应表：一对多

如果已有的GO注释格式为：GeneID和GOID一对多，如下所示：

第一列为GeneID，第二列为GOID，其中多个GOID之间以逗号(,)分隔。

以txt格式存储，列与列之间以tab（制表符）分隔。

GeneID GOID
Os01g0100100 GO:0090630,GO:0005096,GO:0016021,GO:0005507,GO:0006886,GO:0016491
Os01g0100400 GO:0005507,GO:0046658,GO:0006508,GO:0016491,GO:0009506,GO:0004185,GO:0016021
Os01g0100500 GO:0005507,GO:0016021,GO:0016491
Os01g0100600 GO:0003676,GO:0016310,GO:0016301
Os01g0100650 GO:0016021
Os01g0100700 GO:0015935,GO:0003735,GO:0006412,GO:0003723,GO:0022626,GO:0000028

示例文件链接

使用GOEnrich_customMapping函数，mappingfile定义文件路径：

deg = read.delim("demo/demo.rice.deg.txt", header = F)$V1
GOresult = GOEnrich_customMapping(deglist = deg, useFDR = F, cut = 0.05, mappingfile = "demo/demo.GOMapping.txt", outdir = "output", outprefix = "map.GO")
Processing BP class.........done!
Processing MF class.........done!
Processing CC class.........done!

GO enrichment Done!
Result files saved at:
output/map.GO.GO.significant.tsv
output/map.GO.GO.all.tsv
GObubble(dataset = GOresult, onlySig = T, useFDR = F, cut = 0.05, top = 10, TrimTerm = F, TermWidth = 50, MainTitle = "GO plot", ColorScheme = "rw", ColorReverse = F)

选项 含义 可选范围
deglist 差异基因列表变量 必填
mappingfile 一对多的注释文件路径 必填
useFDR 是否使用FDR作为显著性
具体显著性阈值由cut来定义		 T或F
cut 显著性阈值，当useFDR=T时表示FDR阈值
当useFDR=F时表示Pvalue阈值		 0~1的小数
outdir 输出结果目录 必填
outprefix 输出结果文件前缀 必填

4.2.3.4 GO二级节点注释条形图

绘制GO数据库的所有二级节点的所有基因与差异基因条形图：

GOBar(deg, customMapping = "demo/demo.GOMapping.txt")

4.2.4 基于blastkoala的KEGG富集

blastkoala工具可以将上传的序列与KEGG数据库中的序列进行blast比对，得到GeneID的对应关系（GeneID->K number），进而得到GeneID与KEGG Pathway的对应关系。

第一列为GeneID，第二列为K number，文件以txt格式存储，列与列之间以tab（制表符）分隔，无表头。 如下表所示：

Pcab047734
Pcab047735
Pcab047736
Pcab047737
Pcab047738
Pcab047748 K11649

示例文件链接

使用KEGGenrich_blastkoala函数，以blastkoalafile定义下载的文件路径：

deg = read.delim("demo/demo.rice.deg.txt", header = F)$V1
KEGGresult = KEGGenrich_blastkoala(deglist = deg, useFDR = F, blastkoalafile = "demo/demo.blastkoala.txt", outdir = "output", outprefix = "blastkoala")
Enrichment Done! Result files saved:
output/blastkoala.tsv
output/blastkoala.html
KEGGbubble(dataset = KEGGresult, top = 10, MainTitle = "KEGG pathway Enrichment plot", xangle = 45, ColorScheme = "rw", ColorReverse = F)

选项 含义 可选范围
deglist 差异基因列表变量 必填
blastkoalafile blastkoala注释输出文件路径 必填
useFDR 是否输出FDR值 T或F
outdir 输出结果目录 必填
outprefix 输出结果文件前缀 必填