RNA-seq下游数据分析-ballgown到报告。 以Rscript为主,对接PGx RNA-seq choppy现有pipeline,到生成RNA-seq分析报告所需的rds和csv文件。
r
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
yingyu 1de1106973 上传文件至 '' 2 weeks ago
FAQ.md 更新 'FAQ.md' 9 months ago
ID_convert_table.rds 上传文件至 '.' 9 months ago
LICENSE Initial commit 10 months ago
News.md 更新 'News.md' 9 months ago
README.md 更新 'README.md' 7 months ago
RNAseq_1_expr_ballgown.R 上传文件至 '.' 9 months ago
RNAseq_1_expr_stringtie.R 上传文件至 '.' 9 months ago
RNAseq_2_pca.R 上传文件至 '.' 9 months ago
RNAseq_3_cor.R 上传文件至 '.' 9 months ago
RNAseq_4_pwDEG.R 上传文件至 '.' 9 months ago
RNAseq_5_pwGSEA.R 上传文件至 '.' 9 months ago
RNAseq_6_enrichfunc.R 上传文件至 '' 2 weeks ago
human_c2_v5p2.rdata 上传文件至 '.' 9 months ago
human_c5_v5p2.rdata 上传文件至 '.' 9 months ago
index.md 更新 'index.md' 9 months ago

README.md

RNAseqDownstream2report

[TOC]

RNA-seq下游数据分析-ballgown到报告。 以Rscript为主,对接PGx RNA-seq choppy现有pipeline,到生成RNA-seq分析报告所需的rds和csv文件。

整体流程图

包括一下几个文件:

  1. RNAseq_1_expr_stringtie.R:从stringtie的文件夹到基因表达水平表格(每列为样本,每行为基因)
  2. RNAseq_1_expr_ballgown.R:从ballgown的文件夹到基因表达水平表格(每列为样本,每行为基因)
  3. RNAseq_2_PCA.R : 计算PCA。
  4. RNAseq_3_cor.R:计算correlation。
  5. RNAseq_4_pwDEG.R:根据分组信息,计算两两差异信息。
  6. RNAseq_5_pwGSEA.R:根据基因表达水平基于GSEA进行通路分析。
  7. 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

  1. 准备文件:

    1. ballgown/stringtie 文件夹
    2. summary_group 样本group信息
  2. 确认事项:

    1. 所在机器的R/bioconductot安装包已完成(请参考library节查看)

    服务器:10.157.72.53已完成包的安装。

    1. 所在机器联网(RNAseq_6_enrichFunc.R需联网计算)

    PGx服务器每天需重新联网。

    1. 代码、rdata、rds数据均已在运行目录下。

      1. 代码:RNAseq开头的1-6*.R
      2. rds: ID_convert_table.rds
      3. rdata: human_c2_v5p2.rdata、human_c5_v5p2.rdata
  3. 运行以下命令:

   # 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
  1. 将结果上传至报告系统所在服务器,并编辑{报告模板.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

输出结果

该代码将输出一系列输出文件:

  1. 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”

  1. 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

  1. 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

  1. 火山图
   @scatter-plot(dataFile='/yourdir/data/rnaseq_AvsB_degs.rds', dataType='rds', xAxis='logfc', xTitle='log2FC',yAxis='log10p',yTitle='-log10 (p)', colorAttr='sigene')
  1. 差异基因清单表
   @data-table-js(dataUrl='/yourdir/data/rnaseq_degs_acrossgroups.csv')
  1. 整体差异基因数量
   @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 

输出结果

  1. rnaseq_gsea_curatedgenesets.csv 基于curated gene sets的GSEA富集结果。

  2. 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")