ichorCNA working module

demo/Snakefile

@@ -37,6 +37,7 @@ configfile: "../modules/utils/2.0/config/default.yaml"
 configfile: "../modules/bwa_mem/1.0/config/default.yaml"
 configfile: "../modules/liftover/1.0/config/default.yaml"
 configfile: "../modules/controlfreec/1.1/config/default.yaml"
+configfile: "../modules/ichorcna/1.0/config/default.yaml"
 
 
 # Load project-specific config, which includes the shared 
@@ -67,7 +68,7 @@ include: "../modules/utils/2.0/utils.smk"
 include: "../modules/bwa_mem/1.0/bwa_mem.smk"
 include: "../modules/liftover/1.0/liftover.smk"
 include: "../modules/controlfreec/1.1/controlfreec.smk"
-
+include: "../modules/ichorcna/1.0/ichorcna.smk"
 
 
 ##### TARGETS ######
@@ -83,4 +84,5 @@ rule all:
         rules._strelka_all.input,
         rules._bwa_mem_all.input,
         rules._liftover_all.input,
-        rules._controlfreec_all.input
+        rules._controlfreec_all.input,
+        rules._ichorcna_all.input

demo/config.yaml

@@ -77,5 +77,8 @@ lcr-modules:
         tool: "battenberg"
         inputs:
             sample_seg: "data/{tool}/hg38/{tumour_sample_id}--{normal_sample_id}_subclones.igv.seg"
-
-
+
+    ichorcna:
+        inputs:
+            sample_bam: "data/{sample_id}.bam"
+            sample_bai: "data/{sample_id}.bam.bai"

modules/ichorcna/1.0/config/default.yaml

@@ -16,6 +16,7 @@ lcr-modules:
                 chrs:
                     hg19: "1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,X,Y"
                     grch37: "1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,X,Y"
+                    hs37d5: "1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,X,Y"
                     hg38: "chr1,chr2,chr3,chr4,chr5,chr6,chr7,chr8,chr9,chr10,chr11,chr12,chr13,chr14,chr15,chr16,chr17,chr18,chr19,chr20,chr21,chr22,chrX,chrY"
                     grch38: "chr1,chr2,chr3,chr4,chr5,chr6,chr7,chr8,chr9,chr10,chr11,chr12,chr13,chr14,chr15,chr16,chr17,chr18,chr19,chr20,chr21,chr22,chrX,chrY"
                 qual: 20
@@ -45,24 +46,28 @@ lcr-modules:
                 ichorCNA_centromere:  
                     grch37: "{MODSDIR}/src/inst/extdata/GRCh37.p13_centromere_UCSC-gapTable.txt"
                     hg19: "{MODSDIR}/src/inst/extdata/GRCh37.p13_centromere_UCSC-gapTable.txt"
+                    hs37d5: "{MODSDIR}/src/inst/extdata/GRCh37.p13_centromere_UCSC-gapTable.txt"
                     grch38: "{MODSDIR}/src/inst/extdata/GRCh38.GCA_000001405.2_centromere_acen.txt"
                     hg38: "{MODSDIR}/src/inst/extdata/GRCh38.GCA_000001405.2_centromere_acen.txt"
                 ichorCNA_minMapScore: 0.75
                 ichorCNA_chrs:  
                     grch37: "c('1', '2', '3','4','5','6','7','8','9','10','11','12','13','14','15','16','17','18','19','20','21','22','X')"
                     hg19: "c('1', '2', '3','4','5','6','7','8','9','10','11','12','13','14','15','16','17','18','19','20','21','22','X')"
+                    hs37d5: "c('1', '2', '3','4','5','6','7','8','9','10','11','12','13','14','15','16','17','18','19','20','21','22','X')"
                     grch38: "c('chr1', 'chr2', 'chr3','chr4','chr5','chr6','chr7','chr8','chr9','chr10','chr11','chr12','chr13','chr14','chr15','chr16','chr17','chr18','chr19','chr20','chr21','chr22','chrX')"
                     hg38: "c('chr1', 'chr2', 'chr3','chr4','chr5','chr6','chr7','chr8','chr9','chr10','chr11','chr12','chr13','chr14','chr15','chr16','chr17','chr18','chr19','chr20','chr21','chr22','chrX')"
                 ichorCNA_fracReadsInChrYForMale: 0.002
                 ichorCNA_genomeStyle:  # can set this to UCSC or NCBI
                     grch37: "NCBI"
                     hg19: "NCBI"
+                    hs37d5: "NCBI"
                     grch38: "UCSC"
                     hg38: "UCSC"
                 # chrs used for training ichorCNA parameters, e.g. tumor fraction. 
                 ichorCNA_chrTrain:  
                     grch37: "c(1:22)"
                     hg19: "c(1:22)"
+                    hs37d5: "c(1:22)"
                     grch38: "paste0('chr', c(1:22))"
                     hg38: "paste0('chr', c(1:22))"
                 # non-tumor fraction parameter restart values; higher values should be included for cfDNA
@@ -103,8 +108,10 @@ lcr-modules:
         resources:
             readcounter:
                 mem_mb: 2000
+                bam: 1
             run: 
                 mem_mb: 2000
+                bam: 1
 
         pairing_config:
             genome:

modules/ichorcna/1.0/ichorcna.smk

@@ -44,6 +44,7 @@ localrules:
 ##### RULES #####
 # ---------------------------------------------------------------------------- #
 
+CFG["runs"]["binSize"] = str(CFG["options"]["readcounter"]["binSize"])
 
 # Symlinks the input files into the module results directory (under '00-inputs/')
 rule _ichorcna_input_bam:
@@ -52,10 +53,12 @@ rule _ichorcna_input_bam:
         bai = CFG["inputs"]["sample_bai"]
     output:
         bam = CFG["dirs"]["inputs"] + "bam/{seq_type}--{genome_build}/{sample_id}.bam",
-        bai = CFG["dirs"]["inputs"] + "bam/{seq_type}--{genome_build}/{sample_id}.bam.bai" # specific to readCounter
+        bai = CFG["dirs"]["inputs"] + "bam/{seq_type}--{genome_build}/{sample_id}.bam.bai", # specific to readCounter
+        crai = CFG["dirs"]["inputs"] + "bam/{seq_type}--{genome_build}/{sample_id}.bam.crai" 
     run:
         op.relative_symlink(input.bam, output.bam)
         op.relative_symlink(input.bai, output.bai)
+        op.relative_symlink(input.bai, output.crai)
 
 
 rule _ichorcna_read_counter:
@@ -173,7 +176,7 @@ rule _ichorcna_all:
             pair_status=CFG["runs"]["pair_status"],
             tumour_id=CFG["runs"]["tumour_sample_id"],
             normal_id=CFG["runs"]["normal_sample_id"],
-            binSize=str(CFG["options"]["readcounter"]["binSize"]))
+            binSize=CFG["runs"]["binSize"])
 
 
 

modules/ichorcna/1.0/src/inst/extdata/HD_ULP_PoN_hs37d5_1Mb_median_normAutosome_median.rds

@@ -0,0 +1 @@
+HD_ULP_PoN_hg19_1Mb_median_normAutosome_median.rds

modules/ichorcna/1.0/src/inst/extdata/HD_ULP_PoN_hs37d5_500kb_normAutosome_median.rds

@@ -0,0 +1 @@
+HD_ULP_PoN_hg19_500kb_normAutosome_median.rds

modules/ichorcna/1.0/src/inst/extdata/gc_hs37d5_1000kb.wig

@@ -0,0 +1 @@
+gc_hg19_1000kb.wig

modules/ichorcna/1.0/src/inst/extdata/gc_hs37d5_10kb.wig

@@ -0,0 +1 @@
+gc_hg19_10kb.wig

modules/ichorcna/1.0/src/inst/extdata/gc_hs37d5_500kb.wig

@@ -0,0 +1 @@
+gc_hg19_500kb.wig

modules/ichorcna/1.0/src/inst/extdata/gc_hs37d5_50kb.wig

@@ -0,0 +1 @@
+gc_hg19_50kb.wig

modules/ichorcna/1.0/src/inst/extdata/map_hs37d5_1000kb.wig

@@ -0,0 +1 @@
+map_hg19_1000kb.wig

modules/ichorcna/1.0/src/inst/extdata/map_hs37d5_10kb.wig

@@ -0,0 +1 @@
+map_hg19_10kb.wig

modules/ichorcna/1.0/src/inst/extdata/map_hs37d5_500kb.wig

@@ -0,0 +1 @@
+map_hg19_500kb.wig

modules/ichorcna/1.0/src/inst/extdata/map_hs37d5_50kb.wig

@@ -0,0 +1 @@
+map_hg19_50kb.wig

modules/ichorcna/1.0/src/runIchorCNA.R

@@ -125,7 +125,8 @@ if (!is.null(libdir) && libdir != "None"){
 }
 
 ## load seqinfo 
-seqinfo <- getSeqInfo(genomeBuild, genomeStyle)
+# seqinfo <- getSeqInfo(genomeBuild, genomeStyle)
+seqinfo <- NULL
 
 if (substr(tumour_file,nchar(tumour_file)-2,nchar(tumour_file)) == "wig") {
   wigFiles <- data.frame(cbind(patientID, tumour_file))

modules/ichorcna/CHANGELOG.md

@@ -5,10 +5,19 @@ All notable changes to the `ichorcna` module will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [1.0] - 2020-06-28
+## [1.0] - 2021-03-31
 
 This release was authored by jawong.
 
-<!-- TODO: Explain each important module design decision below. -->
+IchorCNA requires indices to be ".bam.bai" format. Modified the initial symlink rule appropriately.
+Modified runIchorCNA.R - line 128: seqinfo <- NULL - to allow for flexibility in genome builds (i.e. hs37d5).
 
-- No module design decisions explained here yet.
+IchorCNA is currently set up to run in unpaired mode using a panel of normals. It currently does not work with CRAM files.
+
+The output files in 99-outputs include:
+- .cna.seg - per-bin CNA state and log ratio
+- .seg - segments called by the Viterbi algorithm - IGV compatible
+- .segTxt - same as .seg but also includes subclonal status of segments (0 = clonal; 1 = subclonal) - not IGV compatible
+- .genomeWide.pdf - CNA profile
+- .param.txt - parameters including ploidy, tumour fraction, sex
+- .corrDepth.txt - per-bin size corrected depth