lbwang/cptac_rna_expression

CPTAC RNA表达 pipeline

本pipeline从GDC对齐的链特异性RNA-Seq BAM文件中生成每个基因的读计数。

• pipeline地址：https://github.com/ding-lab/cptac_rna_expression
• 作者：Liang-Bo Wang <***>

版本发布

• v1.0 (2019-10-18)：初始版本

处理说明

读计数通过featureCounts（subread v1.6.4）在链特异性模式下生成，参数为：-g gene_id -t exon -Q 10 -p -B -s 2。之后使用https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/#upper-quartile-fpkm将读计数转换为FPKM和FPKM-UQ。

使用的基因注释与GDC一致（gencode.v22.annotation.gtf.gz；md5：291330bdcff1094bc4d5645de35e0871），可在https://gdc.cancer.gov/about-data/data-harmonization-and-generation/gdc-reference-files页面获取。请注意，尽管基因注释与GDC相同，但GDC的读计数是通过HTSeq在非链特异性模式下完成的，因此本pipeline与GDC报告的读计数并不完全一致。

输出格式

每个样本会生成一个TSV文件，基因顺序完全相同。输出TSV文件包含以下列：

可使用https://gdc.cancer.gov/about-data/data-harmonization-and-generation/gdc-reference-files页面上的基因信息（https://api.gdc.cancer.gov/data/b011ee3e-14d8-4a97-aed4-e0b10f6bbe82；md5：0a3f1d9b0a679e2a426de36d8d74fbf9）获取每个基因的额外信息。

每批执行会生成一个TSV文件analysis_summary.dat，包含所有输出，列如下：

安装

使用Docker Hub上的镜像 https://hub.docker.com/r/lbwang/cptac_rna_expression。具体工具要求可查看Dockerfile。

输入准备

pipeline需要以下输入，通过snakemake_config.json传递：

json
{
    "sample_list": "样本名称列表",
    "bam_map": "（CPTAC）BAM map的路径",
    "gdc_gtf": "GDC注释GTF的路径（例如gencode.v22.annotation.gtf）",
    "gdc_gene_info": "GDC基因信息TSV的路径（例如gencode.gene.info.v22.tsv）",
    "workflow_root": "本仓库的本地文件路径（例如/path/to/cptac_rna_expression）"
}

Pipeline执行

克隆本仓库为pipeline_workflow：

git clone https://github.com/ding-lab/cptac_rna_expression pipeline_workflow

文件夹结构应如下：

<root>
- pipeline_workflow
- <batch1>
- <batch2>

在katmai上运行pipeline

若conda环境不存在，创建环境：

conda create -n cptac_expression python=3.7 \
    snakemake-minimal=5.10.0 \
    pandas=1.0.1 \
    samtools=1.9 htslib=1.9 \
    subread=1.6.4

创建新的批处理文件夹，例如cptac2_prospective：

mkdir cptac2_prospective

创建snakemake_config.json：

json
{
    "sample_list": "samples.list",
    "bam_map": "bammap.converted.tsv",
    "gdc_gtf": "/diskmnt/Datasets/Reference/GDC/gencode.v22.annotation.gtf",
    "gdc_gene_info": "/diskmnt/Datasets/Reference/GDC/gencode.gene.info.v22.tsv",
    "workflow_root": "/diskmnt/Projects/cptac_scratch/CPTAC_expression/pipeline_workflow"
}

运行完整pipeline：

snakemake \
    --configfile=snakemake_config.json \
    -s ../pipeline_workflow/Snakefile \
    -j 50 --resources io_heavy=4 -- \
    make_analysis_summary

在compute1上运行pipeline

参考示例路径：/storage1/fs1/lding/Active/CPTAC3/Analysis/rna_expression_pipeline/2020-02-25_PDA。

创建新批处理：

cd /storage1/fs1/path/to/rna_expression_pipeline/

# 创建批处理名称
export BATCH="2020-02-25_PDA"

# 创建批处理文件夹
mkdir $BATCH
cd $BATCH

# 创建samples.list。例如：
cut -f2 CPTAC3.catalog/BamMap/compute1.BamMap.dat | tail +2 > samples.list

# 缓存BAM map，建议在文件名中添加CPTAC3.catalog的当前提交哈希
cp CPTAC3.catalog/BamMap/compute1.BamMap.dat compute1.BamMap.cptac3_catalog_commit_6509b69.dat

# 创建snakemake_config.json
{
    "sample_list": "samples.list",
    "bam_map": "compute1.BamMap.cptac3_catalog_commit_6509b69.dat",
    "gdc_gtf": "../annotations/gencode.v22.annotation.gtf",
    "gdc_gene_info": "../annotations/gencode.gene.info.v22.tsv",
    "workflow_root": "/storage1/fs1/lding/Active/CPTAC3/Analysis/rna_expression_pipeline/pipeline_workflow"
}

# 设置日志结构
mkdir logs
mkdir logs/cluster

# 复制bash脚本run.sh和run_master_job.sh

设置snakemake配置文件。默认配置文件位于/storage1/fs1/lding/Active/CPTAC3/Analysis/rna_expression_pipeline/ris_lsf。必要时修改以下文件：

• config.yaml控制作业失败时的重启次数和提交的总作业数。

• lsf-submit.py设置额外的LSF选项。注意作业组是必需的。

  python
  def generate_lsf_command(job_properties: dict, cluster: dict) -> str:
      # ...
      # RIS compute1命令
      queue = "general"
      job_grp = "/liang-bo.wang/cptac_rna"
      docker_img = "lbwang/cptac_rna_expression"
      # ...
  

若使用新的snakemake配置文件，修改run_master_job.sh以指向新配置文件。

运行pipeline：

bash run.sh

将pipeline应用于其他数据集

BAM文件应对齐到hg38（最好与https://gdc.cancer.gov/about-data/data-harmonization-and-generation/gdc-reference-files一致）。bam_map参数提供的BAM map应模仿CPTAC3目录的BAM map列，必须包含以下TSV格式的列：

• # sample_name：注意# 是列名的一部分；样本名称在BAM map中必须唯一
• case
• disease
• UUID：若UUID不可用，可使用任何标识符
• data_pth：样本BAM的绝对路径；文件名必须为*.rna_seq.genomic.gdc_realn.bam
• experimental_strategy：必须为RNA-Seq
• data_format：必须为BAM
• reference：必须为hg38