yingyu
1de1106973
上传文件至 ''
|
5 年之前 | |
|---|---|---|
| FAQ.md | 6 年之前 | |
| ID_convert_table.rds | 6 年之前 | |
| LICENSE | 6 年之前 | |
| News.md | 6 年之前 | |
| README.md | 6 年之前 | |
| RNAseq_1_expr_ballgown.R | 6 年之前 | |
| RNAseq_1_expr_stringtie.R | 6 年之前 | |
| RNAseq_2_pca.R | 6 年之前 | |
| RNAseq_3_cor.R | 6 年之前 | |
| RNAseq_4_pwDEG.R | 6 年之前 | |
| RNAseq_5_pwGSEA.R | 6 年之前 | |
| RNAseq_6_enrichfunc.R | 5 年之前 | |
| human_c2_v5p2.rdata | 6 年之前 | |
| human_c5_v5p2.rdata | 6 年之前 | |
| index.md | 6 年之前 | |
README.md
RNAseqDownstream2report
[TOC]
RNA-seq下游数据分析-ballgown到报告。 以Rscript为主，对接PGx RNA-seq choppy现有pipeline，到生成RNA-seq分析报告所需的rds和csv文件。
整体流程图
包括一下几个文件：
- RNAseq_1_expr_stringtie.R：从stringtie的文件夹到基因表达水平表格（每列为样本，每行为基因）
- RNAseq_1_expr_ballgown.R：从ballgown的文件夹到基因表达水平表格（每列为样本，每行为基因）
- RNAseq_2_PCA.R ： 计算PCA。
- RNAseq_3_cor.R：计算correlation。
- RNAseq_4_pwDEG.R：根据分组信息，计算两两差异信息。
- RNAseq_5_pwGSEA.R：根据基因表达水平基于GSEA进行通路分析。
- RNAseq_6_enrichFunc.R：根据差异基因进行GO和KEGG通路分析。
graph LR;
A1[input: stringtie dir]-->B1[RNAseq_1_expr_stringtie.R];
A2[input: ballgown dir]-->B2[RNAseq_1_expr_ballgown.R];
B1-->B{Expression profile};
B2-->B;
B-->C[RNAseq_2_PCA.R];
B-->D[RNAseq_3_cor.R];
B-->E[RNAseq_4_pwDEG.R];
B-->F[RNAseq_5_pwGSEA.R];
E-->G[RNAseq_6_enrichFunc.R];
Quick start
准备文件：
- ballgown/stringtie 文件夹
- summary_group 样本group信息
确认事项：
- 所在机器的R/bioconductot安装包已完成（请参考library节查看）
服务器：10.157.72.53已完成包的安装。
- 所在机器联网（RNAseq_6_enrichFunc.R需联网计算）
PGx服务器每天需重新联网。
代码、rdata、rds数据均已在运行目录下。
- 代码：RNAseq开头的1-6*.R
- rds: ID_convert_table.rds
- rdata: human_c2_v5p2.rdata、human_c5_v5p2.rdata
运行以下命令：
# RNAseq_1_expr_stringtie.R与 RNAseq_1_expr_ballgown.R 2选1
Rscript RNAseq_1_expr_stringtie.R
Rscript RNAseq_1_expr_ballgown.R -i ./ballgown/
Rscript RNAseq_2_pca.R -i rnaseq_geneexp_log2FPKM.txt -g group.txt
Rscript RNAseq_3_cor.R -i rnaseq_geneexp_log2FPKM.txt -g group.txt
Rscript RNAseq_4_pwDEG.R -i rnaseq_geneexp_log2FPKM.txt -g group.txt -a TRUE -b TRUE -f -4.05
Rscript RNAseq_5_pwGSEA.R -i rnaseq_geneexp_log2FPKM.txt -g group.txt
Rscript RNAseq_6_enrichfunc.R -i rnaseq_degs_acrossgroups.csv
- 将结果上传至报告系统所在服务器，并编辑{报告模板.md}的第一部分(实验目的、主要结论)，完成并运行报告。
RNAseq_1_expr_stringtie.R
功能简介
从stringtie的输出结果到基因表达水平表格（每列为样本，每行为基因）。
代码参数
Usage: Rscript RNAseq_1_expr_stringtie.R [options]
Options:
-o OUT_DIR, --out_dir=OUT_DIR
The output directory [default ./]
-i INPUT_DIR, --input_dir=INPUT_DIR
The directory input of expression files. It is output from stringtie software named as ".gene.abundance.txt". [default ./]
-f NUMBER, --floor_value=NUMBER
A number to add to each value before log2 transformation to avoid infinite value.[default: 0.01]
-l TRUE, --log2_norm=TRUE
Perform log2 transformation on FPKM/TPM value. [default: TRUE]
-s SAMPLE_NAME, --sample_name=SAMPLE_NAME
File in tab-delimited format for sample name if usr want to rename sample name. The input file containing sample name as file name and sample name to be renamed. [default NULL]
-p PROJECT_CODE, --project_code=PROJECT_CODE
Project code, which is used as prefix of output file. [default: rnaseq]
-h, --help
Show this help message and exit
参数解释
| 参数 | 取值类型 | 解释 | 例如 |
|---|---|---|---|
| -o OUT_DIR, --out_dir=OUT_DIR | character | 输出路径，默认为./。可加“/”也可不加“/” | ./ |
| -i INPUT, --input=INPUT | character | 输入文件所在路径，stringtie输出文件所在的文件夹。默认会读取所有该路径下所有”.gene.abundance.txt”.的文件进行合并。 | ./ |
| -f NUMBER, --floor_value=NUMBER | number | 在log2转换时，需在0中加入一个底值。默认为0.01 | 0.01 |
| -l TRUE, --log2_norm=TRUE | TRUE/FALSE | 是否进行log2转换 | TRUE |
| -s SAMPLE_NAME, --sample_name=SAMPLE_NAME | character | 文件名与样本名对应文件。如需重命名样本，可以输入对应文件。文件要求：tab分隔的文件。第一列为原始文件名。第二列为修改的样本名（表达谱的列名）。第一行为标题。 | NULL |
| -p PROJECT_CODE, --project_code=PROJECT_CODE | character | project代号，输出文件的前缀，默认rnaseq | rnaseq |
| -h, --help | 查看帮助文档并退出 | -h |
输出结果
tab分隔的基因以FPKM和TPM标准化的表达谱。
文件名为：rnaseq_geneexp_log2FPKM(log2TPM).txt 和(-l TRUE /默认)或rnaseq_geneexp_FPKM(TPM).txt (-l FALSE) 。
文件内容例如：
Gene P1 P2 P3 ENSG00000000003 2.951 5.085 3.592 ENSG00000000005 -6.644 -4.248 -3.085 ENSG00000000419 4.966 6.197 5.332 ENSG00000000457 0.854 1.838 0.665 ENSG00000000460 -0.19 1.693 0.145 ENSG00000000938 -6.644 -6.148 -6.644 ENSG00000000971 -5.919 -6.644 -6.644 ENSG00000001036 5.134 5.47 4.998 ENSG00000001084 2.676 2.303 2.638
运行示例
#最少输入
Rscript RNAseq_1_expr_ballgown.R -i ./ballgown/
# 其他输入
Rscript RNAseq_1_expr_ballgown.R -o yourdir -i ./ballgown/ -l FALSE -p test
RNAseq_1_expr_ballgown.R
功能简介
从ballgown的文件夹到基因表达水平表格（每列为样本，每行为基因）
代码参数
Usage: Rscript RNAseq_1_expr_ballgown.R [options]
Options:
-o OUT_DIR, --out_dir=OUT_DIR
The output directory [default ~]
-i INPUT, --input=INPUT
The directory input of expression files. It is output from ballgown software.
-f NUMBER, --floor_value=NUMBER
A number to add to each value before log2 transformation to avoid infinite value.[default: 0.01]
-l TRUE, --log2_norm=TRUE
Perform log2 transformation on FPKM value. [default: TRUE]
-p PROJECT_CODE, --project_code=PROJECT_CODE
Project code, which is used as prefix of output file. [default: rnaseq]
-h, --help
Show this help message and exit
参数解释
| 参数 | 取值类型 | 解释 | 例如 |
|---|---|---|---|
| -o OUT_DIR, --out_dir=OUT_DIR | character | 输出路径，默认为./。可加“/”也可不加“/” | ./ |
| -i INPUT, --input=INPUT | character | 输入路径， ballgown的文件夹，必须输入。要求ballgown文件夹中文件如下所示<http://www.bioconductor.org/packages/release/bioc/vignettes/ballgown/inst/doc/ballgown.html | ./ballgown/ |
| -f NUMBER, --floor_value=NUMBER | number | 在log2转换时，需在0中加入一个底值。默认为0.01 | 0.01 |
| -l TRUE, --log2_norm=TRUE | TRUE/FALSE | 是否进行log2转换 | TRUE |
| -p PROJECT_CODE, --project_code=PROJECT_CODE | character | project代号，输出文件的前缀，默认rnaseq | rnaseq |
| -h, --help | 查看帮助文档并退出 | -h |
输出结果
tab分隔的基因表达谱，文件名为：rnaseq_geneexp_log2FPKM.txt (-l TRUE /默认)或rnaseq_geneexp_FPKM.txt (-l FALSE) 。
文件内容例如：
Gene P1 P2 P3 ENSG00000000003 2.951 5.085 3.592 ENSG00000000005 -6.644 -4.248 -3.085 ENSG00000000419 4.966 6.197 5.332 ENSG00000000457 0.854 1.838 0.665 ENSG00000000460 -0.19 1.693 0.145 ENSG00000000938 -6.644 -6.148 -6.644 ENSG00000000971 -5.919 -6.644 -6.644 ENSG00000001036 5.134 5.47 4.998 ENSG00000001084 2.676 2.303 2.638
运行示例
#最少输入
Rscript RNAseq_1_ballgown.R -i ./ballgown/
# 其他输入
Rscript RNAseq_1_ballgown.R -o /home/yuying/rnaseqreport_test -i ./ballgown/ -l FALSE -p test
RNAseq_2_pca.R
功能简介
计算PCA，输出choppy report所需的scatterplot图的rds和csv文件。
代码参数
Usage: Rscript RNAseq_2_pca.R [options]
Options:
-o OUT_DIR, --out_dir=OUT_DIR
The output directory [default ~]
-i INPUT, --input=INPUT
The input expression files. required!
-g SAMPLE_GROUP, --sample_group=SAMPLE_GROUP
File for sample group infomation.The input file containing sample name and group infomation. note colname must be like: sample group1 group2...
-p PROJECT_CODE, --project_code=PROJECT_CODE
Project code, which is used as prefix of output file. [default: rnaseq]
-h, --help
Show this help message and exit
| 参数 | 取值类型 | 解释 | 例如 |
|---|---|---|---|
| -o OUT_DIR, --out_dir=OUT_DIR | character | 输出路径，默认为./。可加“/”也可不加“/” | ./ |
| -i INPUT, --input=INPUT | character | 输入文件名，必须输入。输入表达谱必须是log scaled的tab分隔的表达谱，可以是RNAseq_1_ballgown.R的输出文件。 | example.txt |
| -g SAMPLE_GROUP, --sample_group=SAMPLE_GROUP | character | 有tab分隔的样本的分组信息，一行为一个样本，每列为分组信息。分组信息可以是多列。如果没有，可以不输入。 | group.txt |
| -p PROJECT_CODE, --project_code=PROJECT_CODE | character | project代号，输出文件的前缀，默认rnaseq | rnaseq |
| -h, --help | 查看帮助文档并退出 | -h |
输出结果
各样本的各PC值，choppy report所需的scatterplot图的rds和csv文件，（其中绘图时仅需rds文件，csv文件就看看）：
rnaseq_pca.csv：逗号分隔的文件 rnaseq_pca.rds：R对象
内容如下：
”“,“PC1”,“PC2”,“PC3”,“sample”,“group1”,“group2” “P1”,-135.940,-151.769,-4.017e-15,“P1”,“A”,“test1” “P2”,259.848,-4.758,-1.906e-13,“P2”,“B”,“test2” “P3”,-123.908,156.528,2.464e-13,“P3”,“B”,“test1”
运行示例
#最少输入
Rscript RNAseq_2_pca.R -i ballgown_geneexp_log2fpkm_floor0p01_c3r58395_2019-04-29.txt
#其他输入
Rscript RNAseq_2_pca.R -o -i ballgown_geneexp_log2fpkm_floor0p01_c3r58395_2019-04-29.txt -g group2.txt -p test
choppy report
有group信息：
@scatter-plot(dataFile='/yourdir/data/rnaseq_pca.rds', dataType='rds', xAxis='PC1', xTitle='',yAxis='PC2',yTitle="",colorAttr='group1')
无group信息：
@scatter-plot(dataFile='/yourdir/data/rnaseq_pca.rds', dataType='rds', xAxis='PC1', xTitle='',yAxis='PC2',yTitle='')
RNAseq_3_cor.R
功能简介
计算correlation，输出choppy report所需的scatterplot图的rds和csv文件。
代码参数
Usage: Rscript RNAseq_3_cor.R [options]
Options:
-o OUT_DIR, --out_dir=OUT_DIR
The output directory [default ./]
-i INPUT, --input=INPUT
The input expression files. required!
-g SAMPLE_GROUP, --sample_group=SAMPLE_GROUP
File for sample group infomation.The input file containing sample name and group infomation. note colname must be like: sample group1 group2...
-p PROJECT_CODE, --project_code=PROJECT_CODE
Project code, which is used as prefix of output file. [default: rnaseq]
-h, --help
Show this help message and exit
| 参数 | 取值类型 | 解释 | 例如 |
|---|---|---|---|
| -o OUT_DIR, --out_dir=OUT_DIR | character | 输出路径，默认为./。可加“/”也可不加“/” | ./ |
| -i INPUT, --input=INPUT | character | 输入文件名，必须输入。输入表达谱必须是log scaled的tab分隔的表达谱，可以是RNAseq_1_ballgown.R的输出文件。 | example.txt |
| -g SAMPLE_GROUP, --sample_group=SAMPLE_GROUP | character | 有tab分隔的样本的分组信息，一行为一个样本，每列为分组信息。分组信息可以是多列。这项输入在该代码中与输出无关。该参数设置仅为串行的逻辑完整性。 | group.txt |
| -p PROJECT_CODE, --project_code=PROJECT_CODE | character | project代号，输出文件的前缀，默认rnaseq | rnaseq |
| -h, --help | 查看帮助文档并退出 | -h |
输出结果
样本两两关系的peason correlation matrix，choppy report所需的heatmap图的rds和csv文件，（其中绘图时仅需rds文件，csv文件就看看）。
rnaseq_cor.csv：逗号分隔的文件 rnaseq_cor.rds：R对象
内容如下：
”“,“P1”,“P2”,“P3” “P1”,1,0.898,0.944 “P2”,0.898,1,0.901 “P3”,0.944,0.901,1
运行示例
#最少输入
Rscript RNAseq_3_cor.R -i ballgown_geneexp_log2fpkm_floor0p01_c3r58395_2019-04-29.txt
#其他输入
Rscript RNAseq_3_cor.R -o -i ballgown_geneexp_log2fpkm_floor0p01_c3r58395_2019-04-29.txt -g group2.txt -p test
choppy report
@heatmap-d3(dataFile='/yourdir/data/rnaseq_cor.rds', dataType='rds', labCol='TRUE')
RNAseq_4_pwDEG.R
功能简介
根据sample_group的分组信息，两两计算差异基因，输出差异基因列表和火山图所需文件，用于choppy报告系统。差异基因的cutoff：t test p<0.05 同时 log2 fold change >=1 或<= (-1)
代码参数
Usage: Rscript RNAseq_4_pwDEG.R [options]
Options:
-o OUT_DIR, --out_dir=OUT_DIR
The output directory [default ./]
-i INPUT, --input=INPUT
The input expression files. required!
-g SAMPLE_GROUP, --sample_group=SAMPLE_GROUP
File for sample group infomation.The input file containing sample name and group infomation. note colname must be like: sample group1 group2...
-p PROJECT_CODE, --project_code=PROJECT_CODE
Project code, which is used as prefix of output file. [default: rnaseq]
-a FALSE, --output_all_genes=FALSE
Output rds files for choppy contains all genes. By default, only DEGs are listed in the output rds and csv for report. NOTE choppy report may not be availble to display correctly if too many points exit. [default: FALSE]
-b FALSE, --low_expr_filter=FALSE
Conduct low expression filtering before DEG analysis. [default: FALSE]
-f NUMBER, --low_expr_filter_cutoff=NUMBER
Genes across all samples with expreesion lower than this value will be filtered out [default: 0]
-h, --help
Show this help message and exit
| 参数 | 取值类型 | 解释 | 例如 |
|---|---|---|---|
| -o OUT_DIR, --out_dir=OUT_DIR | character | 输出路径，默认为./。可加“/”也可不加“/” | ./ |
| -i INPUT, --input=INPUT | character | 输入文件名，必须输入。输入表达谱必须是log scaled的tab分隔的表达谱，可以是RNAseq_1_ballgown.R的输出文件。 | example.txt |
| -g SAMPLE_GROUP, --sample_group=SAMPLE_GROUP | character | 有tab分隔的样本的分组信息，必须输入。格式为：每行一个样本，每列为分组信息。分组信息可以是多列。 | group.txt |
| -p PROJECT_CODE, --project_code=PROJECT_CODE | character | project代号，输出文件的前缀，默认rnaseq | rnaseq |
| -a FALSE, --output_all_genes=FALSE | TRUE/FALSE | 是否输出所有基因 | FALSE |
| -b FALSE, --low_expr_filter=FALSE | TRUE/FALSE | 在差异其因寻找之前是否进行低表达基因过滤 | FALSE |
| -f NUMBER, --low_expr_filter_cutoff=NUMBER | number | 在所有样本中，低于-f的表达水平的基因会作为低表达基因过滤掉 | 0 |
| -h, --help | 查看帮助文档并退出 | -h |
输出结果
该代码将输出一系列输出文件：
- PROJECT_CODE_GROUPversus_deg.csv(rds)
根据分组情况进行两两比较，每次比较的结果输出1个PROJECT_CODE_GROUPversus_deg.csv，一个PROJECT_CODE_GROUPversus_deg.rds。因此根据比较次数，将输出比较次数*2个文件。
choppy report所需的火山图的rds和csv文件，（其中绘图时仅需rds文件，csv文件就看看）。
例如：
rnaseq_AvsB_degs.csv
rnaseq_AvsB_degs.rds
文件内容为：每次比较的差异基因、比较的组别信息、log2FC、-log10 P value、是否为DEG。默认情况下只输出DEG，若计算时设置参数“-a TRUE”，则输出所有基因，sigene标记为nonDEG和DEG两类。
“gene”,“versus”,“logfc”,“log10p”,“sigene” “ENSG00000004776”,“A vs B”,-2.82273809523809,1.50125858539637,“DEG” “ENSG00000007171”,“A vs B”,3.06833333333333,1.62715851425345,“DEG” “ENSG00000011347”,“A vs B”,-1.10492857142857,1.31435980565864,“DEG” “ENSG00000011677”,“A vs B”,-2.19445238095238,1.55184587515544,“DEG”
- PROJECT_CODE_deg_acrossgroups.csv
将1产生的所有差异基因（仅差异基因）集合到一起，用于差异基因清单表展示。
“gene”,“versus”,“pvalue”,“log2FC” “ENSG00000004776”,“A vs B”,0.03153,-2.823 “ENSG00000007171”,“A vs B”,0.0236,3.068 “ENSG00000011347”,“A vs B”,0.04849,-1.105 “ENSG00000011677”,“A vs B”,0.02806,-2.194 “ENSG00000012504”,“A vs B”,0.04607,3.357
- PROJECT_CODE_GROUPversus_deg_stats.csv(rds)
差异基因统计，用于choppy报告系统中的barplot图。
“number”,“type”,“versus” 65,“upregulated”,“A vs B” 138,“downregulated”,“A vs B”
运行示例
#最少输入
Rscript RNAseq_4_pwDEG.R -i example_geneexp_log2fpkm_floor0p01_c13r58395_2019-04-30.txt -g group13_1.txt
#其他输入
Rscript RNAseq_4_pwDEG.R -i example_geneexp_log2fpkm_floor0p01_c13r58395_2019-04-30.txt -g group13_1.txt -p example_2 --low_expr_filter=TRUE -f -1
choppy report
- 火山图
@scatter-plot(dataFile='/yourdir/data/rnaseq_AvsB_degs.rds', dataType='rds', xAxis='logfc', xTitle='log2FC',yAxis='log10p',yTitle='-log10 (p)', colorAttr='sigene')
- 差异基因清单表
@data-table-js(dataUrl='/yourdir/data/rnaseq_degs_acrossgroups.csv')
- 整体差异基因数量
@stack-barplot-r(dataFile='/yourdir/data/rnaseq_degs_stats.rds', dataType='rds', xAxis='versus', yAxis='number',labelAttr='type',barPos='stack')
RNAseq_5_pwGSEA.R
功能简介
利用fgsea包对不同比较的基因进行GSEA通路分析。GSEA原理请参看：https://www.pnas.org/content/102/43/15545
注意准备参考文件。
代码参数
Usage: Rscript RNAseq_5_pwGSEA.R [options]
Options:
-o OUT_DIR, --out_dir=OUT_DIR
The output directory [default ./]
-i INPUT, --input=INPUT
The input expression files. Required!
-e TYPE_GENE_ID, --type_gene_id=TYPE_GENE_ID
The type of gene symbol. Could be either of EnsemblGID/EntrezID/GeneSymbol [default: EnsemblGID]
-g SAMPLE_GROUP, --sample_group=SAMPLE_GROUP
File for sample group infomation.The input file containing sample name and group infomation. note colname must be like: sample group1 group2... Required!
-f NUMBER, --padjvalueCutoff=NUMBER
Cutoff value of adjusted p value. [default: 0.2]
-p PROJECT_CODE, --project_code=PROJECT_CODE
Project code, which is used as prefix of output file. [default: rnaseq]
-d REF_RDATA_DIR, --ref_rdata_dir=REF_RDATA_DIR
The directory of reference files: human_c2_v5p2.rdata, human_c5_v5p2.rdata and ID_convert_table.rds. [default: ./]
-h, --help
Show this help message and exit
| 参数 | 取值类型 | 解释 | 例如 |
|---|---|---|---|
| -o OUT_DIR, --out_dir=OUT_DIR | character | 输出路径，默认为./。可加“/”也可不加“/” | ./ |
| -i INPUT, --input=INPUT | character | 输入文件名，必须输入。输入表达谱必须是log scaled的tab分隔的表达谱，可以是RNAseq_1_ballgown.R的输出文件。 | example.txt |
| -e TYPE_GENE_ID, --type_gene_id=TYPE_GENE_ID | character | 基因ID类型，可以是：Ensembl gene ID (EnsemblGID)、Entrez Gene ID (EntrezID)或Gene Symbol (GeneSymbol)。[default: EnsemblGID] | EnsemblGID |
| -g SAMPLE_GROUP, --sample_group=SAMPLE_GROUP | character | 有tab分隔的样本的分组信息，必须输入。格式为：每行一个样本，每列为分组信息。分组信息可以是多列。 | group.txt |
| -q NUMBER, --padjvalueCutoff=NUMBER | number | 富集分析adjust p值 cutoff | 0.2 |
| -p PROJECT_CODE, --project_code=PROJECT_CODE | character | project代号，输出文件的前缀，默认rnaseq | rnaseq |
| -d REF_RDATA_DIR, --ref_rdata_dir=REF_RDATA_DIR | character | 参考文件所在路径，human_c2_v5p2.rdata, human_c5_v5p2.rdata and ID_convert_table.rds. [default: ./] | ./ |
| -h, --help | 查看帮助文档并退出 | -h |
运行示例
Rscript RNAseq_5_pwGSEA.R -o /home/yuying/rnaseqreport_test -i example_geneexp_log2fpkm_floor0p01_c13r58395_2019-04-30.txt -g group13_1.txt
输出结果
rnaseq_gsea_curatedgenesets.csv 基于curated gene sets的GSEA富集结果。
rnaseq_gsea_go.csv 基于GO 功能的GSEA富集结果。
”“,“versus”,“pathway”,“pval”,“padj”,“ES”,“NES”,“nMoreExtreme”,“size”,“leadingEdge” “1”,“test1 vs test2”,“GO_REGULATION_OF_CELL_ACTIVATION”,0.001007,0.1232,-0.4512,-1.458,0,458,“2302, 3127, 3123, 6352, 3119, 912, 972, 2625, 3122, 958, 8808, 53833, 284021, 10451, 154, 301, 55024, 3109, 923, 348, 8456, 84433, 51237, 4323, 7409, 919, 11005, 3606, 727897, 3929, 7056, 114548, 857, 10673, 695, 163747, 3120, 683, 124912, 29126, 114771, 2150, 3113, 3956, 22914, 4773, 83639, 634, 1029, 1236, 29108, 90865, 441478, 3600, 5896, 51083, 89780, 10148, 3965, 9173, 2323, 84106, 51744, 282618, 2852, 2056, 1269, 5592, 5724, 3623, 3567, 11314, 23529, 7474, 558, 10461, 283234, 11148, 5341, 3273, 9466, 22890, 22806, 917, 8832, 5588, 3952, 282616, 2207, 3111, 3574, 84807, 2268, 282617, 6869, 3659, 6441, 8772, 3575, 8546, 1948, 246778, 1178, 5585, 84959, 2064” “2”,“test1 vs test2”,“GO_IMMUNE_EFFECTOR_PROCESS”,0.001007,0.1232,-0.4688,-1.515,0,455,“5473, 6374, 3127, 1380, 3123, 8519, 3119, 972, 2625, 10581, 958, 722, 284021, 3437, 3627, 8284, 56892, 10451, 3075, 1755, 3428, 51191, 4353, 644150, 28984, 4939, 55061, 114836, 717, 117157, 9245, 8809, 566, 7409, 919, 154064, 3929, 23705, 340061, 55601, 114548, 4599, 9844, 730, 2633, 3433, 10410, 3764, 2150, 7098, 91607, 60489, 3956, 3439, 3426, 29108, 90865, 10584, 725, 10417, 3078, 715, 3383, 84871, 3434, 10964, 51744, 64218, 1621, 282618, 23586, 3665, 4600, 78989, 710, 733, 3815, 3936, 1636”
每列的含义：
A table with GSEA results. Each row corresponds to a tested pathway. The columns are the following:
- versus - compared group
- pathway – name of the pathway as in ‘names(pathway)’;
- pval – an enrichment p-value;
- padj – a BH-adjusted p-value;
- ES – enrichment score, same as in Broad GSEA implementation;
- NES – enrichment score normalized to mean enrichment of random samples of the same size;
- nMoreExtreme’ – a number of times a random gene set had a more extreme enrichment score value;
- size – size of the pathway after removing genes not present in ‘names(stats)’.
- leadingEdge – vector with indexes of leading edge genes that drive the enrichment, see http://software.broadinstitute.org/gsea/doc/GSEAUserGuideTEXT.htm#_Running_a_Leading.
choppy report
@data-table-js(dataUrl='/yourdir/data/rnaseq_gsea_curatedgenesets.csv')
@data-table-js(dataUrl='/yourdir/data/rnaseq_gsea_go.csv')
RNAseq_6_enrichfunc.R
功能简介
利用 clusterProfiler包对不同比较中的差异基因进行GO和KEGG通路分析。
输入：RNAseq_4_pwDEG.R输出的差异基因清单表（PROJECT_CODE_deg_acrossgroups.csv）。
注意联网
本功能较慢，每组的分析约需5分钟。
代码参数
Usage: Rscript RNAseq_6_enrichfunc.R [options]
Options:
-o OUT_DIR, --out_dir=OUT_DIR
The output directory [default ./]
-i INPUT, --input=INPUT
The input DEG list in csv format. The first column: gene; second column: group. Required!
-e TYPE_GENE_ID, --type_gene_id=TYPE_GENE_ID
The type of gene symbol. Could be either of EnsemblGID/EntrezID/GeneSymbol [default: EnsemblGID]
-f NUMBER, --pvalueCutoff=NUMBER
Cutoff value of p value. [default: 0.05]
-m PADJUSTMETHOD, --pAdjustMethod=PADJUSTMETHOD
Method of adjust p value. One of "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none". [default: BH]
-q NUMBER, --qvalueCutoff=NUMBER
Cutoff value of q value. [default: 0.2]
-p PROJECT_CODE, --project_code=PROJECT_CODE
Project code, which is used as prefix of output file. [default: rnaseq]
-h, --help
Show this help message and exit
| 参数 | 取值类型 | 解释 | 例如 |
|---|---|---|---|
| -o OUT_DIR, --out_dir=OUT_DIR | character | 输出路径，默认为./。可加“/”也可不加“/” | ./ |
| -i INPUT, --input=INPUT | character | 输入文件名，必须输入。输入表达谱必须是log scaled的tab分隔的表达谱，可以是RNAseq_1_ballgown.R的输出文件。 | example.txt |
| -e TYPE_GENE_ID, --type_gene_id=TYPE_GENE_ID | character | 基因ID类型，可以是：Ensembl gene ID (EnsemblGID)、Entrez Gene ID (EntrezID)或Gene Symbol (GeneSymbol)。[default: EnsemblGID] | EnsemblGID |
| -f NUMBER, --pvalueCutoff=NUMBER | number | 富集分析p值 cutoff | 0.05 |
| -m PADJUSTMETHOD, --pAdjustMethod=PADJUSTMETHOD | character | p adjust method | BH |
| -q NUMBER, --qvalueCutoff=NUMBER | number | 富集分析q值 cutoff | 0.2 |
| -p PROJECT_CODE, --project_code=PROJECT_CODE | character | project代号，输出文件的前缀，默认rnaseq | rnaseq |
| -h, --help | 查看帮助文档并退出 | -h |
输出结果
GO和KEGG通路结果。
内容如下：
”“,“versus”,“ID”,“Description”,“GeneRatio”,“BgRatio”,“pvalue”,“p.adjust”,“qvalue”,“geneID”,“Count” “1”,“A vs B”,“hsa05168”,“Herpes simplex virus 1 infection”,“39/185”,“492/7847”,9.625e-12,2.117e-09,2.057e-09,“256051/684/7568/10172/84765/80095/162963/84436/55762/55769/55786/90594/148268/342908/30832/348327/100129543/81931/390927/147837/57573/388566/91120/113835/163059/84671/65251/79973/126017/147949/374900/7594/3111/3135/728927/100129842/59348/26974/7772”,39
运行示例
#最少输入
Rscript RNAseq_6_enrichfunc.R -i rnaseq_degs_acrossgroups.csv
choppy report
@data-table-js(dataUrl='/yourdir/data/rnaseq_GOenrich.csv')
@data-table-js(dataUrl='/yourdir/data/rnaseq_KEGGenrich.csv')
Example data
oss://choppy-app-example-data/RNAseqDownstream2report_exampledata/ 13例RNA-seq数据。
R library
R library used:
library(optparse) #all
library(ballgown) #RNAseq_1_ballgown.R
library(fgsea) #RNAseq_5_pwGSEA.R
library(clusterProfiler) #RNAseq_6_enrichfunc.R
library(org.Hs.eg.db) #RNAseq_5_pwGSEA.R &RNAseq_6_enrichfunc.R
如果需要安装，
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ballgown")
BiocManager::install("org.Hs.eg.db")
BiocManager::install("fgsea")
BiocManager::install("clusterProfiler")
install.packages("optparse")