perf(CNVkit): genemetrics

tasks/CNVkit.wdl

     
     # Derive each segment's absolute integer copy number, ploidy must be int value
     PURITY=`awk -F'\t' '{print $6}' ${hrd} | sed -n '2p'`
-    cnvkit.py call ${sample}.cns -y -m clonal --purity $PURITY -o ${sample}.call.cns
-    cnvkit.py call ${sample}.cnr -y -m clonal --purity $PURITY -o ${sample}.call.cnr
-
+    cnvkit.py segmetrics ${sample}.cnr -s ${sample}.cns --ci -o ${sample}.segmetrics.cns
+    cnvkit.py call ${sample}.segmetrics.cns --drop-low-coverage --filter ci -m threshold --purity $PURITY -o ${sample}.call.cns
+    
     # Plot the results
     cnvkit.py scatter ${sample}.cnr -s ${sample}.call.cns -o ${sample}.scatter.pdf
     cnvkit.py diagram ${sample}.cnr -s ${sample}.call.cns -o ${sample}.diagram.pdf
     cnvkit.py heatmap ${sample}.cnr ${sample}.call.cns -o ${sample}.heatmap.pdf
     
     # Genemetrics
-    mkdir gainloss
-    cnvkit.py genemetrics ${sample}.call.cnr -t 0.2 -m 3 -o ${sample}.ratio_cnv.txt
-    cnvkit.py genemetrics ${sample}.call.cnr -s ${sample}.call.cns -t 0.2 -m 3 -o ${sample}.segment_cnv.txt
+    cnvkit.py genemetrics ${sample}.cnr -t 0.2 -m 3 -o ${sample}.ratio_cnv.txt
+    cnvkit.py genemetrics ${sample}.cnr -s ${sample}.call.cns -t 0.2 -m 3 -o ${sample}.segment_cnv.txt
     # Filter genes
     cat ${sample}.ratio_cnv.txt | tail -n+2 | cut -f1 | sort | uniq > ratio_cnv.txt
     cat ${sample}.segment_cnv.txt | tail -n+2 | cut -f1 | sort | uniq > segment_cnv.txt
     comm -12 ratio_cnv.txt segment_cnv.txt > ${sample}.trusted_genes.txt
-    # Scatter plot for each gene
-    touch failed_genes.txt
+    
+    # # Scatter plot for each gene
+    # mkdir gainloss
+    # touch failed_genes.txt
     # for gene in `cat ${sample}.trusted_genes.txt`
     # do
-    #   cnvkit.py scatter ${sample}.call.cnr -s ${sample}.call.cns -g $gene -o ./gainloss/${sample}.$gene.scatter.pdf || echo $gene >> failed_genes.txt
+    #   cnvkit.py scatter ${sample}.cnr -s ${sample}.call.cns -g $gene -o ./gainloss/${sample}.$gene.scatter.pdf || echo $gene >> failed_genes.txt
     # done
+    
+    # Filter by trusted_genes
+    awk 'NR==FNR {a[$1]=$2;next} NR!=FNR {if(FNR == 1 || (FNR != 1 && $1 in a)) print $0}' ${sample}.trusted_genes.txt ${sample}.ratio_cnv.txt > ${sample}.ratio_cnv.trusted.txt
+    # Infer absolute CN (not adjust by purity)
+    cnvkit.py call ${sample}.ratio_cnv.trusted.txt -m threshold -o ${sample}.ratio_cnv.call.txt
+    awk '{if ($6 != 2) print $0}' ${sample}.ratio_cnv.call.txt > ${sample}.ratio_cnv.call.filter.txt
   >>>
   
   runtime {
     File cnr = "${sample}.cnr"
     File cns = "${sample}.cns"
     File stats = "${sample}.stats"
-    File call_cnr = "${sample}.call.cnr"
     File call_cns = "${sample}.call.cns"
     File ratio_cnv = "${sample}.ratio_cnv.txt"
     File segment_cnv = "${sample}.segment_cnv.txt"
     File trusted_genes = "${sample}.trusted_genes.txt"
-    File failed_genes = "${sample}.failed_genes.txt"
-    Array[File] gainloss = glob("./gainloss/*")
+    File? failed_genes = "${sample}.failed_genes.txt"
+    Array[File]? gainloss = glob("./gainloss/*")
+    File ratio_cnv_trusted = "${sample}.ratio_cnv.trusted.txt"
+    File ratio_cnv_trusted_call = "${sample}.ratio_cnv.call.txt"
+    File ratio_cnv_trusted_call_filter = "${sample}.ratio_cnv.call.filter.txt"
   }
 }