Hi all, I posted a few days ago about my nf-core/sarek run using controlfreec failing in makeGraph2.0.R. Error was complaining about certain lines in my ratio file not having 12 columns. After some more digging I think I've found a possible solution and am sharing it here. Turns out the ratio file has a few lines that read like "Pathogenic (nssv579289)" - the space is being interpreted by read.table (makeGraph2.0.R, line 31) as a new column. I suggested a possible fix to makeGraph2.0.R (https://github.com/BoevaLab/FREEC/issues/163) - think specifying sep="\t" should sort the issue. Line 31:

ratio <-read.table(args[ratioFileInd], header=TRUE);

Suggested fix:

ratio <-read.table(args[ratioFileInd], header=TRUE, sep="\t");

Details:

Command executed:

  cat $(which makeGraph2.0.R) | R --slave --args  sample89.tumor.mpileup.gz_ratio.txt sample89.tumor.mpileup.gz_BAF.txt

  mv *_BAF.txt.png sample89_BAF.png
  mv *_ratio.txt.log2.png sample89_ratio.log2.png
  mv *_ratio.txt.png sample89_ratio.png

  cat <<-END_VERSIONS > versions.yml
  "NFCORE_SAREK:SAREK:BAM_VARIANT_CALLING_TUMOR_ONLY_ALL:BAM_VARIANT_CALLING_TUMOR_ONLY_CONTROLFREEC:MAKEGRAPH2":
      controlfreec: 11.6b
  END_VERSIONS

Command exit status:
  1

Command output:
  (empty)

Command error:
  Error in scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :
    line 224918 did not have 12 elements
  Calls: read.table -> scan
  Execution halted

(base) [ec2-user@ip-172-31-20-221 b158d738818e329f11c7e47dd9eef0]$ awk 'NF != 12' sample89.tumor.mpileup.gz_ratio.txt 
Y	17877785	0.272727	0.554471	1	-1	1	A	-1	 Pathogenic (nssv579289)	0	0
Y	17877835	0.0181818	0.554471	1	-1	1	A	-1	 Pathogenic (nssv579289)	0	0
Y	18028145	0.327273	0.554471	1	-1	1	A	-1	 Pathogenic (nssv579289)	0	0
Y	18028178	0.4	0.554471	1	-1	1	A	-1	 Pathogenic (nssv579289)	0	0
Y	21381697	0.290909	0.554471	1	-1	1	A	-1	 Pathogenic (nssv579289)	0	0
(base) [ec2-user@ip-172-31-20-221 b158d738818e329f11c7e47dd9eef0]$ head -n 1 sample89.tumor.mpileup.gz_ratio.txt 
Chromosome	Start	Ratio	MedianRatio	CopyNumber	BAF	estimatedBAF	Genotype	UncertaintyOfGT	Gene	Subclone_CN	Subclone_Population
(base) [ec2-user@ip-172-31-20-221 b158d738818e329f11c7e47dd9eef0]$ awk 'NF == 12' sample89.tumor.mpileup.gz_ratio.txt | head -n 3
Chromosome	Start	Ratio	MedianRatio	CopyNumber	BAF	estimatedBAF	Genotype	UncertaintyOfGT	Gene	Subclone_CN	Subclone_Population
1	65565	0.0545455	1.23636	2	-1	2	-	-1	1:65564-65573	0	0
1	358067	0.0181818	1.23636	2	-1	2	-	-1	1:358066-358183	0	0

Yup, that makes sense, clearly a case of sloppy programming. (BTW, I always use

read.delim()

in my R scripts, that already has

sep = "\t"

as default.) Not sure if/when they'll fix it, given that the last commit was in 2023 and they have >90 open issues. Where does the "Pathogenic (nssv579289)" come from? Is there a way you can replace the space (e.g. with an underscore) as a work-around?

Thanks for the tip on read.delim()! I couldn't figure out where the "Pathogenic (nssv579289)" was coming from - searching the FREEC source for "Pathogenic" doesn't turn up anything. Thought about using a config override to convert to underscore, but wondering if it's worth the effort over switching to something more actively maintained like CNVkit...

It certainly wouldn't come from the source code, but presumably some reference file or database (ClinVar?), so if you can control the reference file, you could change it (not sure how you would do that with a "config"). I suspect this an exception, since ControlFREEC runs just fine most of the time, so it must have been tested either with a different reference file or the variants that trigger these lines in the file are pretty rare ...

... but wondering if it's worth the effort over switching to something more actively maintained ...

Fortuantely sarek has options. I have used CNVkit with WES (to do that properly you should run it on the duplicate-marked BAMs, not the recalibrated ones) a bit, and it worked alright, but I haven't dealt much with copy number variation and don't feel qualified to make recommendations.

Your comments inspired me to persevere just a bit more. I traced it back to the reference BED file that I was feeding into sarek using the --intervals option -- a BED file from the Twist Exome WES kit that was used in the sequencing (https://www.twistbioscience.com/resources/data-files/twist-exome-20-bed-files). Checking this bed file picked up the "Pathogenic (nssv...)" ..!

[alpha@alpha null]$ grep "Pathogenic" bedfile.bed 
chrY	17877784	17877829	Deprecated_content:ClinGen CNVs: Pathogenic (nssv579289)
chrY	17877834	17877862	Deprecated_content:ClinGen CNVs: Pathogenic (nssv579289)
chrY	18028144	18028172	Deprecated_content:ClinGen CNVs: Pathogenic (nssv579289)
chrY	18028177	18028222	Deprecated_content:ClinGen CNVs: Pathogenic (nssv579289)
chrY	21381696	21381743	Deprecated_content:ClinGen CNVs: Pathogenic (nssv579289)

Fixing this bed reference file, replacing

Pathogenic (nssv579289)

with

Pathogenic_(nssv579289)

should solve the problem. I haven't re-run this yet to verify, but I think can safely consider this one solved. Hope it helps anyone else working with the Twist Exome WES kit and controlfreec.

Actually, what I normally do is

cut -f 1-3 old.bed > new.bed

to avoid sarek misinterpreting an extra column as timing info and going off the rails by spawning tens of thousands of processes (that should be in the 5th column, i.e. "score", I think, but I had bad experience with a WES targets file from Agilent). But CNVkit (haven't tried ControlFREEC) actually use the the 4th column ("name") to assign gene names to the loci it find (that's optional, but often helpful). Assuming ControlFREE does the same, you have to think about what is useful. If you want those

nssv

numbers and

Pathogenic

designations in the output, then yes, something like

sed s/ /_/g

. would do it. But if not,

cut -f 1-3

may be easier and cleaner, or -- if you have that information -- you could put gene names into the 4th column (they shouldn't have spaces).

Thanks for the tip on just going with the first 3 columns! Sounds like the right way to go here...