|
|
@@ -10,23 +10,22 @@ RNA-seq下游数据分析-ballgown到报告。 以Rscript为主,对接PGx RNA- |
|
|
|
|
|
|
|
包括一下几个文件:
|
|
|
|
|
|
|
|
1. **RNAseq_1_ballgown.R**:从ballgown的文件夹到基因表达水平表格(每列为样本,每行为基因)
|
|
|
|
|
|
|
|
2. **RNAseq_2_PCA.R** : 计算PCA。
|
|
|
|
|
|
|
|
3. **RNAseq_3_cor.R**:计算correlation。
|
|
|
|
|
|
|
|
4. **RNAseq_4_pwDEG.R**:根据分组信息,计算两两差异信息。
|
|
|
|
|
|
|
|
5. **RNAseq_5_pwGSEA.R**:根据基因表达水平基于GSEA进行通路分析。
|
|
|
|
|
|
|
|
6. **RNAseq_6_enrichFunc.R**:根据差异基因进行GO和KEGG通路分析。
|
|
|
|
1. **RNAseq_1_expr_stringtie.R**:从ballgown的文件夹到基因表达水平表格(每列为样本,每行为基因)
|
|
|
|
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通路分析。
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
```mermaid
|
|
|
|
graph LR;
|
|
|
|
A(input: ballgown dir)-->B[RNAseq_1_ballgown.R];
|
|
|
|
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];
|
|
|
@@ -62,11 +61,14 @@ graph LR; |
|
|
|
3. 运行以下命令:
|
|
|
|
|
|
|
|
```shell
|
|
|
|
Rscript RNAseq_1_ballgown.R -i ./ballgown/
|
|
|
|
Rscript RNAseq_2_pca.R -i rnaseq_geneexp_log2fpkm_floor0p01_c13r58395_2019-04-30.txt -g group.txt
|
|
|
|
Rscript RNAseq_3_cor.R -o -i rnaseq_geneexp_log2fpkm_floor0p01_c13r58395_2019-04-30.txt -g group.txt
|
|
|
|
Rscript RNAseq_4_pwDEG.R -i rnaseq_geneexp_log2fpkm_floor0p01_c13r58395_2019-04-30.txt -g group.txt -a TRUE -b TRUE -f -4.05
|
|
|
|
Rscript RNAseq_5_pwGSEA.R -i rnaseq_geneexp_log2fpkm_floor0p01_c13r58395_2019-04-30.txt -g group.txt
|
|
|
|
# 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
|
|
|
|
```
|
|
|
|
|
|
|
@@ -75,7 +77,90 @@ graph LR; |
|
|
|
4. 将结果上传至报告系统所在服务器,并编辑{报告模板.md}的第一部分(实验目的、主要结论),完成并运行报告。
|
|
|
|
|
|
|
|
|
|
|
|
## RNAseq_1_ballgown.R
|
|
|
|
## RNAseq_1_expr_stringtie.R
|
|
|
|
|
|
|
|
### 功能简介
|
|
|
|
|
|
|
|
从stringtie的输出结果到基因表达水平表格(每列为样本,每行为基因)。
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
### 代码参数
|
|
|
|
|
|
|
|
```shell
|
|
|
|
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
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
### 运行示例
|
|
|
|
|
|
|
|
```shell
|
|
|
|
#最少输入
|
|
|
|
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
|
|
|
|
|
|
|
|
### 功能简介
|
|
|
|
|
|
|
@@ -86,7 +171,7 @@ graph LR; |
|
|
|
### 代码参数
|
|
|
|
|
|
|
|
```shell
|
|
|
|
Usage: Rscript RNAseq_1_ballgown.R [options]
|
|
|
|
Usage: Rscript RNAseq_1_expr_ballgown.R [options]
|
|
|
|
|
|
|
|
Options:
|
|
|
|
-o OUT_DIR, --out_dir=OUT_DIR
|
|
|
@@ -107,7 +192,7 @@ Options: |
|
|
|
-h, --help
|
|
|
|
Show this help message and exit
|
|
|
|
```
|
|
|
|
参数解释
|
|
|
|
**参数解释**
|
|
|
|
|
|
|
|
| 参数 | 取值类型 | 解释 | 例如 |
|
|
|
|
| -------------------------------------------- | ---------- | ------------------------------------------------------------ | ----------- |
|
|
|
@@ -122,7 +207,7 @@ Options: |
|
|
|
|
|
|
|
### 输出结果
|
|
|
|
|
|
|
|
tab分隔的基因表达谱,文件名为:rnaseq_geneexp_fpkm_c4r58395_2019-04-30.txt
|
|
|
|
tab分隔的基因表达谱,文件名为:rnaseq_geneexp_log2FPKM.txt (-l TRUE /默认)或rnaseq_geneexp_FPKM.txt (-l FALSE) 。
|
|
|
|
|
|
|
|
文件内容例如:
|
|
|
|
|
|
|
@@ -221,11 +306,11 @@ Rscript RNAseq_2_pca.R -o -i ballgown_geneexp_log2fpkm_floor0p01_c3r58395_2019-0 |
|
|
|
|
|
|
|
有group信息:
|
|
|
|
|
|
|
|
@scatter-plot(dataFile='/*yourdir*/data/rnaseq_pca.rds', dataType='rds', xAxis='PC1', xTitle="",yAxis='PC2',yTitle="",colorAttr="group1")
|
|
|
|
``` @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="")
|
|
|
|
``` @scatter-plot(dataFile='/yourdir/data/rnaseq_pca.rds', dataType='rds', xAxis='PC1', xTitle='',yAxis='PC2',yTitle='')```
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
@@ -301,7 +386,7 @@ Rscript RNAseq_3_cor.R -o -i ballgown_geneexp_log2fpkm_floor0p01_c3r58395_2019-0 |
|
|
|
|
|
|
|
### choppy report
|
|
|
|
|
|
|
|
@heatmap-d3(dataFile='/*yourdir*/data/rnaseq_cor.rds', dataType='rds', labCol='TRUE')
|
|
|
|
```@heatmap-d3(dataFile='/yourdir/data/rnaseq_cor.rds', dataType='rds', labCol='TRUE')```
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
@@ -421,15 +506,16 @@ Rscript RNAseq_4_pwDEG.R -i example_geneexp_log2fpkm_floor0p01_c13r58395_2019-0 |
|
|
|
|
|
|
|
1. 火山图
|
|
|
|
|
|
|
|
@scatter-plot(dataFile='/*yourdir*/data/rnaseq_AvsB_degs.rds', dataType='rds', xAxis='logfc', xTitle="log2FC",yAxis='log10p',yTitle="-log10 (p)", colorAttr="sigene")
|
|
|
|
```@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')
|
|
|
|
```@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')
|
|
|
|
``` @stack-barplot-r(dataFile='/yourdir/data/rnaseq_degs_stats.rds', dataType='rds', xAxis='versus', yAxis='number',labelAttr='type',barPos='stack') ```
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
## RNAseq_5_pwGSEA.R
|
|
|
@@ -438,6 +524,8 @@ Rscript RNAseq_4_pwDEG.R -i example_geneexp_log2fpkm_floor0p01_c13r58395_2019-0 |
|
|
|
|
|
|
|
利用fgsea包对不同比较的基因进行GSEA通路分析。GSEA原理请参看:<https://www.pnas.org/content/102/43/15545>
|
|
|
|
|
|
|
|
注意准备参考文件。
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
### 代码参数
|
|
|
@@ -464,22 +552,26 @@ Options: |
|
|
|
|
|
|
|
-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 |
|
|
|
|
| -h, --help | | 查看帮助文档并退出 | -h |
|
|
|
|
| 参数 | 取值类型 | 解释 | 例如 |
|
|
|
|
| ----------------------------------------------- | --------- | ------------------------------------------------------------ | ----------- |
|
|
|
|
| -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 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
@@ -493,8 +585,6 @@ Rscript RNAseq_5_pwGSEA.R -o /home/yuying/rnaseqreport_test -i example_geneexp_l |
|
|
|
|
|
|
|
### 输出结果
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1. rnaseq_gsea_curatedgenesets.csv 基于curated gene sets的GSEA富集结果。
|
|
|
|
|
|
|
|
2. rnaseq_gsea_go.csv 基于GO 功能的GSEA富集结果。
|
|
|
@@ -522,9 +612,9 @@ A table with GSEA results. Each row corresponds to a tested pathway. The columns |
|
|
|
|
|
|
|
### choppy report
|
|
|
|
|
|
|
|
@data-table-js(dataUrl='/*yourdir*/data/rnaseq_gsea_curatedgenesets.csv')
|
|
|
|
``` @data-table-js(dataUrl='/yourdir/data/rnaseq_gsea_curatedgenesets.csv')```
|
|
|
|
|
|
|
|
@data-table-js(dataUrl='/*yourdir*/data/rnaseq_gsea_go.csv')
|
|
|
|
```@data-table-js(dataUrl='/yourdir/data/rnaseq_gsea_go.csv')```
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
@@ -611,9 +701,9 @@ 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_GOenrich.csv')```
|
|
|
|
|
|
|
|
@data-table-js(dataUrl='/*yourdir*/data/rnaseq_KEGGenrich.csv')
|
|
|
|
```@data-table-js(dataUrl='/yourdir/data/rnaseq_KEGGenrich.csv')```
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
@@ -650,5 +740,4 @@ BiocManager::install("clusterProfiler") |
|
|
|
install.packages("optparse")
|
|
|
|
```
|
|
|
|
|
|
|
|
|
|
|
|
|