improve annotation in pipeline

[slsevilla]

Currently annotation calling is one of the largest bottlenecks of the pipeline. It is currently split into several rules and accompanying scripts.

Rules

• peak_Transcripts
• peak_ExonIntron
• peak_RMSK
• peak_Transcripts
• peak_junctions
• peak_process
• project_annotations

Scripts

• 05_Anno_ExonIntron.R 05_Anno_Process.R
• 05_Anno_RMSK.R
• 05_Anno_Transcript.R
• 05_Anno_junctions.R
• 05_get_site2peak_lookup.sh
• 05_jcounts2peakconnections.py
• 05_peak_annotation.R
• 05_peak_annotation_functions.R
• 06_annotation.Rmd

The general workflow is to run each annotation type separately before merging into one RMD file. This requires a significant amount of time, and is generating individual jobs per sample per rule, which also utilizes more Biowulf resources than maybe necessary.

Goals for the re-write

1. Speed up performance
2. Reduce the number of input/output files required for execution
3. Transfer all file creation from R files to snakemake
4. Reduce the number of rules required without sacrificing speed considerably

[slsevilla]

Rule ExonIntron

Project info

Three projects were created from previous runs to complete benchmarking analysis

• project1: mESC_clip_4_v2.0
• project2: 8-09-21-HaCaT_fCLIP_v2.0
• project3: mES_fclip_1_YL_011622_v2.0

File info

• all projects are set-up with the following structure

  ├── proj_number
  │   └── exp_output
  │   └── input

• Required inputs for one sample

  
  └── input
  └── 04_annotation
      ├── 01_project
      │   ├── 7SKRNA_Repeatmasker.bed
      │   ├── annotations.txt
      │   ├── DNA_Repeatmasker.bed
      │   ├── lincRNA_Gencode.bed
      │   ├── LINE\ SINE_Repeatmasker.bed
      │   ├── lncRNA_Gencode.bed
      │   ├── lncRNA_Gencode.txt
      │   ├── Low_complexity_Repeatmasker.bed
      │   ├── LTR_Repeatmasker.bed
      │   ├── miRNA_Gencode.bed
      │   ├── ncRNA_annotations.txt
      │   ├── Other_Repeatmasker.bed
      │   ├── ref_gencode.txt
      │   ├── rRNA_Custom.bed
      │   ├── rRNA_Gencode.bed
      │   ├── rRNA_Repeatmasker.bed
      │   ├── Satellite_Repeatmasker.bed
      │   ├── scRNA_Repeatmasker.bed
      │   ├── Simple_repeat_Repeatmasker.bed
      │   ├── sncRNA_Custom.bed
      │   ├── snoRNA_Gencode.bed
      │   ├── snRNA_Gencode.bed
      │   ├── srpRNA_Repeatmasker.bed
      │   ├── tRNA_Custom.bed
      │   ├── Unknown_Repeatmasker.bed
      │   └── yRNA_Repeatmasker.bed
      └── 02_peaks
          ├── Control1hr_Clip_ALLreadPeaks_AllRegions.txt
          └──  Control7hr_Clip_ALLreadPeaks_AllRegions.txt
  └── config
     └── annotation_config.txt

- Expected outputs for one sample (Control1hr)

├── exp_output │   └── 04_annotation │   └── 02_peaks │   ├── Control1hr_Clip_ALLreadPeaks_AllRegions_IntronExon_OppoStrand.txt │   ├── Control1hr_Clip_ALLreadPeaks_AllRegions_IntronExon_SameStrand.txt │   ├── Control7hr_Clip_ALLreadPeaks_AllRegions_IntronExon_OppoStrand.txt │   └──Control7hr_Clip_ALLreadPeaks_AllRegions_IntronExon_SameStrand.txt

Project input/output files are located here:

/data/RBL_NCI/Wolin/Sam/annotation_testing

## Script calling
Script location:

/data/RBL_NCI/Pipelines/iCLIP/v2.0/workflow/scripts/

Example R script (SameStrand, proj_1):

Rscript /data/RBL_NCI/Pipelines/iCLIP/v2.0/workflow/scripts/05_Anno_ExonIntron.R \ --rscript /data/RBL_NCI/Pipelines/iCLIP/v2.0/workflow/scripts/05_peak_annotation_functions.R \ --peak_type ALL \ --anno_anchor max_total \ --read_depth 3 \ --sample_id Control1hr_Clip \ --ref_species mm10 \ --anno_dir /data/RBL_NCI/Wolin/Sam/annotation_testing/proj_1/input/04_annotation/01_project/ \ --reftable_path /data/RBL_NCI/Wolin/Sam/annotation_testing/proj_1/config/annotation_config.txt \ --gencode_path /data/CCBR_Pipeliner/iCLIP/ref/annotations/mm10/Gencode_VM23/fromGencode/gencode.vM23.annotation.gtf.txt \ --intron_path /data/CCBR_Pipeliner/iCLIP/ref/annotations/mm10/Gencode_VM23/fromUCSC/KnownGene/KnownGene_GRCm38_introns.bed \ --rmsk_path /data/CCBR_Pipeliner/iCLIP/ref/annotations/mm10/repeatmasker/rmsk_GRCm38.txt \ --tmp_dir /data/RBL_NCI/Wolin/Sam/annotation_testing/tmp \ --out_dir /data/RBL_NCI/Wolin/Sam/annotation_testing/proj_1/testing_output/04_annotation/02_peaks/ \ --out_file /data/RBL_NCI/Wolin/Sam/annotation_testing/proj_1/testing_output/04_annotation/02_peaks/Control1hr_Clip_ALLreadPeaks_AllRegions_IntronExon_SameStrand.txt \ --anno_strand "SameStrand"

Example R script (OppoStrand, proj_1):

Rscript /data/RBL_NCI/Pipelines/iCLIP/v2.0/workflow/scripts/05_Anno_ExonIntron.R \ --rscript /data/RBL_NCI/Pipelines/iCLIP/v2.0/workflow/scripts/05_peak_annotation_functions.R \ --peak_type ALL \ --anno_anchor max_total \ --read_depth 3 \ --sample_id Control1hr_Clip \ --ref_species mm10 \ --anno_dir /data/RBL_NCI/Wolin/Sam/annotation_testing/proj_1/input/04_annotation/01_project/ \ --reftable_path /data/RBL_NCI/Wolin/Sam/annotation_testing/proj_1/config/annotation_config.txt \ --gencode_path /data/CCBR_Pipeliner/iCLIP/ref/annotations/mm10/Gencode_VM23/fromGencode/gencode.vM23.annotation.gtf.txt \ --intron_path /data/CCBR_Pipeliner/iCLIP/ref/annotations/mm10/Gencode_VM23/fromUCSC/KnownGene/KnownGene_GRCm38_introns.bed \ --rmsk_path /data/CCBR_Pipeliner/iCLIP/ref/annotations/mm10/repeatmasker/rmsk_GRCm38.txt \ --tmp_dir /data/RBL_NCI/Wolin/Sam/annotation_testing/tmp \ --out_dir /data/RBL_NCI/Wolin/Sam/annotation_testing/proj_1/testing_output/04_annotation/02_peaks/ \ --out_file /data/RBL_NCI/Wolin/Sam/annotation_testing/proj_1/testing_output/04_annotation/02_peaks/Control1hr_Clip_ALLreadPeaks_AllRegions_IntronExon_OppoStrand.txt \ --anno_strand "OppoStrand"

[slsevilla]

Create DAG of pipeline v2.0 for review

dag.pdf

[wilfriedguiblet]

Improved IE_calling speed in 05_peak_annotation_functions.R.