RNAseqDownstream2report/README.md

# RNAseqDownstream2report



[TOC]

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

## 整体流程图

包括一下几个文件：

1. RNAseq_1_ballgown.R

2. RNAseq_2_PCA.R

3. RNAseq_3_cor.R

4. RNAseq_4_pwDEG.R

5. RNAseq_5_pwGSEA.R

6. RNAseq_6_enrichFunc.R

7. RNAseq_7_sexcheck.R


```mermaid
graph LR;
    A(input: ballgown dir)-->B[RNAseq_1_ballgown.R];
    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];
    B-->J[RNAseq_7_sexcheck.R];
```



## RNAseq_1_ballgown.R

### 功能简介

从ballgown的文件夹到基因表达水平表格（每列为样本，每行为基因）



### 代码与参数

```shell
Usage: Rscript RNAseq_1_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_fpkm_c4r58395_2019-04-30.txt

文件内容例如：

> 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



### 运行示例

```shell
#最少输入
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文件。



### 代码与参数

```shell
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"



### 运行示例

```shell
#最少输入
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文件。



### 代码与参数

```shell
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



### 运行示例

```shell
#最少输入
 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)



### 代码参数

```shell
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"

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

3. PROJECT_CODE_GROUPversus_deg_stats.csv(rds)

   差异基因统计，用于choppy报告系统中的barplot图。

   > "number","type","versus"
   > 65,"upregulated","A vs B"
   > 138,"downregulated","A vs B"



### 运行示例

```shell
#最少输入
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")

2. 差异基因清单表

   @data-table-js(dataUrl=/*yourdir*/data/rnaseq_degs_acrossgroups.csv')

3. 整体差异基因数量

   @stack-barplot-r(dataFile='/*yourdir*/data/rnaseq_degs_stats.rds', dataType='rds', xAxis='versus', yAxis='number',labelAttr='type',barPos='stack') 


## RNAseq_5_pwGSEA.R



## RNAseq_6_enrichfunc.R

### 功能简介

利用 clusterProfiler包对不同比较中的差异基因进行GO和KEGG通路分析。

输入：RNAseq_4_pwDEG.R输出的差异基因清单表（PROJECT_CODE_deg_acrossgroups.csv）。

#注意联网

### 代码参数

```shell
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通路结果。



### 运行示例

```shell
#最少输入
Rscript RNAseq_6_enrichfunc.R -i example_3_degs_acrossgroups.csv
```



### choppy report

@data-table-js(dataUrl=/*yourdir*/data/rnaseq_GOenrich.csv')

@data-table-js(dataUrl=/*yourdir*/data/rnaseq_KEGGenrich.csv')



## RNAseq_7_sexcheck.R