gavinha / titancna Goto Github PK
View Code? Open in Web Editor NEWAnalysis of subclonal copy number alterations (CNA) and loss of heterozygosity (LOH) in cancer
License: GNU General Public License v3.0
Analysis of subclonal copy number alterations (CNA) and loss of heterozygosity (LOH) in cancer
License: GNU General Public License v3.0
Hi!
I'm wondering how to use the new GATK4 copy number segmentation as input to TITAN. I guess it should be possible (see eg. https://software.broadinstitute.org/gatk/documentation/article?id=11088) but I'm not sure how. I was able to use a recent pipeline to convert its output to a format that mirrors the Allelic CapSeg output, which could be used for eg. ABSOLUTE. But I'd like to use TITAN.
Any advice?
Thanks!
Using the bioconda release of TitanCNA (https://anaconda.org/bioconda/bioconductor-titancna) I am getting an error if the --libdir
parameter is omitted from the command line:
$ Rscript `which titanCNA.R` --id tmp2 --hetFile tmp.hets.tsv --cnFile tmp2.cr.tsv --numClusters 1 --numCores 1 --normal_0 0.5 --ploidy_0 2 --chrs "c(1:22, \"X\")" --estimatePloidy TRUE --outDir ./
There were 11 warnings (use warnings() to see them)
Error in if (!is.null(libdir) & libdir != "None") { :
argument is of length zero
Execution halted
If I specify --libdir
it works fine.
I suppose that
if (!is.null(libdir) & libdir != "None")
should be changed to
if (!is.null(libdir) && libdir != "None") { :
EDIT: looks also to be the case in current repo version https://github.com/gavinha/TitanCNA/blob/master/scripts/R_scripts/titanCNA.R#L111
Thanks for releasing this great software!
I've been running TitanCNA on bams with coverage 1x, 5x, 10x, and 20x.
There's consistently an error for the low coverage bam files of 1-5x coverage. As an example, here is the output I'm seeing for the low coverage bam files:
Running TITAN...
titan: Loading data input.hets.txt
titan: Loading GC content and mappability corrected log2 ratios...
titan: Extracting read depth...
Removed 0 centromeric positions
titan: Loading default parameters
Loading required package: iterators
Using 4 cores.
titan: Running HMM...
fwdBack: Iteration 1 chr: 1 2 Error in { :
task 1 failed - "INTEGER() can only be applied to a 'integer', not a 'NULL'"
Calls: runEMclonalCN -> %dopar% -> <Anonymous>
Execution halted
TITAN works perfectly for the 20x and 40x coverage files though:
Running TITAN...
titan: Loading data input.hets.txt
titan: Loading GC content and mappability corrected log2 ratios...
titan: Extracting read depth...
Removed 0 centromeric positions
titan: Loading default parameters
Loading required package: iterators
Using 4 cores.
titan: Running HMM...
fwdBack: Iteration 1 chr: 1 2 3 4 8 7 6 5 12 11 9 10 16 15 13 14 20 19 23 18 17 21 22
Coord Descent[1]=-561231 n=0.5 s=[0.001,0.200,0.300] phi=4
Coord Descent[2]=-274262 n=0.2084 s=[0.051,0.204,0.344] phi=4.016
Coord Descent[3]=-272542 n=0.2943 s=[0.05,0.20,0.34] phi=4.014
Coord Descent[4]=-272522 n=0.2948 s=[0.052,0.196,0.341] phi=4.014
Coord Descent[5]=-272522 n=0.2946 s=[0.052,0.196,0.341] phi=4.014
...
Is it possible that runEMclonalCN
is getting stuck at gaps in the coverage? Have you noticed there's a "minimal coverage" necessary for TITAN to work properly?
Why does titanCNA.R
exclude sex chromosomes for males?
TitanCNA/scripts/R_scripts/titanCNA.R
Lines 191 to 194 in 2a7e268
It seems that correctIntegerCN
already adjusts for sex, so no need to remove X for males?
https://github.com/gavinha/TitanCNA/blob/master/R/utils.R#L1213-L1227
# Adjust chrX copy number if purity is sufficiently high
# males - all data points in chrX is corrected
# females - only
if (purity >= minPurityToCorrect){
if (gender == "male" & length(chrXStr) > 0){
segs[Chromosome == chrXStr, Corrected_Copy_Number := as.integer(round(logR_Copy_Number))]
segs[Chromosome == chrXStr, Corrected_Call := names[Corrected_Copy_Number + 2]]
cn[Chr == chrXStr, Corrected_Copy_Number := as.integer(round(logR_Copy_Number))]
cn[Chr == chrXStr, Corrected_Call := names[Corrected_Copy_Number + 2]]
}else if (gender == "female"){
segs[Chromosome == chrXStr & Copy_Number >= maxCNtoCorrect.X, Corrected_Copy_Number := as.integer(round(logR_Copy_Number))]
segs[Chromosome == chrXStr & Copy_Number >= maxCNtoCorrect.X, Corrected_Call := "HLAMP"]
cn[Chr == chrXStr & CopyNumber >= maxCNtoCorrect.X, Corrected_Copy_Number := as.integer(round(logR_Copy_Number))]
cn[Chr == chrXStr & CopyNumber >= maxCNtoCorrect.X, Corrected_Call := "HLAMP"]
}
In a test sample I'm finding that selectSolution
is not selecting "what should biologically speaking" be the best solution.
For some reason, it prefers this ploidy=4, clusters=3 solution:
Normal contamination estimate: 0.6033
Average tumour ploidy estimate: 3.869
Clonal cluster cellular prevalence Z=2: 1 0.7426
AllelicRatio binomial means for clonal cluster Z=1: 0.5 0.6237 0.6983 0.5198 0.7483 0.5828 0.784 0.642 0.5284 0.8109 0.6865 0.5622 0.8318 0.7212 0.6106 0.5332 0.8485 0.749 0.6494 0.5498 0.8623 0.7717 0.6811 0.5906 0.5362
logRatio Gaussian means for clonal cluster Z=1: -1.184 -0.7739 -0.455 -0.455 -0.194 -0.194 0.02699 0.02699 0.02699 0.2186 0.2186 0.2186 0.3877 0.3877 0.3877 0.3877 0.539 0.539 0.539 0.539 0.6759 0.6759 0.6759 0.6759 0.6759
AllelicRatio binomial means for clonal cluster Z=2: 0.5 0.5864 0.6473 0.5147 0.6926 0.5642 0.7276 0.6138 0.5228 0.7554 0.6532 0.5511 0.7781 0.6854 0.5927 0.5278 0.7969 0.7121 0.6272 0.5424 0.8128 0.7346 0.6564 0.5782 0.5313
logRatio Gaussian means for clonal cluster Z=2: -0.9585 -0.6849 -0.455 -0.455 -0.2568 -0.2568 -0.08252 -0.08252 -0.08252 0.07294 0.07294 0.07294 0.2133 0.2133 0.2133 0.2133 0.3411 0.3411 0.3411 0.3411 0.4586 0.4586 0.4586 0.4586 0.4586
logRatio Gaussian variance: 0.009996 0.009996 0.01096 0.01096 0.009909 0.009909 0.01014 0.01014 0.01014 0.01051 0.01051 0.01051 0.009999 0.009999 0.009999 0.009999 0.01001 0.01001 0.01001 0.01001 0.01007 0.01007 0.01007 0.01007 0.01007
Number of iterations: 5
Log likelihood: -52150
S_Dbw dens.bw (LogRatio): 0.1997
S_Dbw scat (LogRatio): 1.0000
S_Dbw validity index (LogRatio): 1.1997
S_Dbw dens.bw (AllelicRatio): 0.4919
S_Dbw scat (AllelicRatio): 1.0000
S_Dbw validity index (AllelicRatio): 1.4919
S_Dbw dens.bw (Both): 0.6916
S_Dbw scat (Both): 2.0000
S_Dbw validity index (Both): 2.6916
Over a more meaningful ploidy = 2 solution like this one:
Normal contamination estimate: 0.401
Average tumour ploidy estimate: 1.911
Clonal cluster cellular prevalence Z=1: 1
AllelicRatio binomial means for clonal cluster Z=1: 0.5 0.7138 0.7995 0.53 0.8457 0.6152 0.8746 0.6873 0.5375 0.8944 0.7366 0.5789 0.9088 0.7725 0.6363 0.5409 0.9197 0.7998 0.6799 0.56 0.9283 0.8212 0.7142 0.6071 0.5428
logRatio Gaussian means for clonal cluster Z=1: -1.279 -0.4744 0.03912 0.03912 0.4171 0.4171 0.7163 0.7163 0.7163 0.964 0.964 0.964 1.175 1.175 1.175 1.175 1.36 1.36 1.36 1.36 1.523 1.523 1.523 1.523 1.523
logRatio Gaussian variance: 0.009996 0.01122 0.01062 0.01062 0.01015 0.01015 0.01004 0.01004 0.01004 0.01 0.01 0.01 0.009996 0.009996 0.009996 0.009996 0.009996 0.009996 0.009996 0.009996 0.009996 0.009996 0.009996 0.009996 0.009996
Number of iterations: 5
Log likelihood: -60370
S_Dbw dens.bw (LogRatio): 0.1482
S_Dbw scat (LogRatio): 1.0000
S_Dbw validity index (LogRatio): 1.1482
S_Dbw dens.bw (AllelicRatio): 0.4143
S_Dbw scat (AllelicRatio): 1.0000
S_Dbw validity index (AllelicRatio): 1.4143
S_Dbw dens.bw (Both): 0.5626
S_Dbw scat (Both): 2.0000
S_Dbw validity index (Both): 2.5626
What parameters are used to select the optimal variant and how can we adjust this? LOH of chromosome arms 1p/19q (as in the second, non-selected solution) is a well recognized marker of this cancer type.
UPDATE: I guess it looks at the log-likelihood. Why do you think the first solution was scored higher than the second solution?
Hi there,
I've git clone's TitanCNA and updated my samples.yaml to reflect the samples I'll be using however, when trying to run through the snakemake I get this error from the first step:
Terminal command:
snakemake -s TitanCNA.snakefile -np
Output:
SyntaxError:
Not all output, log and benchmark files of rule read_counter contain the same wildcards. This is crucial though, in order to avoid that two or more jobs write to the same file.
File "/home/pchamely/TitanCNA/scripts/snakemake/TitanCNA.snakefile", line 4, in
File "/home/pchamely/TitanCNA/scripts/snakemake/ichorCNA.snakefile", line 25, in
This is my first time using snakemake so I'm not sure what this means and was wondering if you could walk me through it?
I had some samples called foo_CAP
and foo_CAP-DEEP
.
The code in selectSolution.R
does a grep on the file list (also valid for ploidy of 3 and 4):
phi2Samples <- grep(id, phi2Files, value=T)
However in this case it will match both files and cause inconsistent results to be written. As the pattern in default cases is id plus an underscore, a solution like this may work:
phi2Samples <- grep(paste0(id, "_"), phi2Files, value=T)
Exact string matching won't help as it will still catch the shorter string contained in the long one.
Just wondering if you can shed some light on the segmentation results.
I know this sample has an amplification of EGFR, chr7:55086714-55324313
However when I examine the .segs.txt file,
chr7 | 15726066 | 50175681 | 164
chr7 | 51096036 | 53103673 | 4
chr7 | 56059591 | 56174140 | 6
The whole gene seems to have been excluded.
Where do the start and stop come from for these segs?
Why are certain parts of the genome excluded?
I'm running PhyloWGS with Titan and I'm getting an error that says I have CNVs with the same start and end locations. When I look as the segs.txt file for that sample, I see that there are indeed several CNVs with very short lengths. This seems strange to me, can anyone suggest a reason this might be happening?
Hello,
The .segs.txt file lists both Copy number, Minor CN, Major CN, and a Corrected_Copy_Number
For the CN, I see that there is a max of CN 8 in the titanCNA script, however the Corrected copy number at one point is up to 84. Given that I've seen massive increase in a few other copy number callers and not just in Titan, so I have reason to believe it's a realistic integer, I was wondering if it would be fair to use this Corrected_Copy_Number in my next analysis, specifically using this as the total copy number and assigning the Major CN and Minor CN to be:
Major_Cn = Corrected_Copy_Number - MinorCN
otherwise, could you maybe elaborate on the Corrected_Copy_Number derivation and how it differs from the Max/Min CN in the other columns.
Thank you!
Hi Gavin,
Apologies for raising so many issues. Completely understand if you don't have time to address everything. I hope that I contribute to improving this excellent software in a way.
I've found some samples that show a homozygous deletion in CDKNA, but that TITAN classifies as HLAMP, despite getting the corrected_copy_number and logr_copy_number correct. Do you know why this could be happening?
pair_barcode | chrom | pos | num_snp | median_ratio | median_logr | titan_state | titan_call | copy_number | minor_cn | major_cn | clonal_cluster | cellular_prevalence | logr_copy_number | corrected_copy_number | corrected_call |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
xxx | 9 | [21900599,21997016) | 26 | 0.854816 | -0.778442 | 20 | ALOH | 8 | 0 | 8 | 3 | 0.21310998 | 0.015625 | 0 | HLAMP |
TCGA-14-1402 | 9 | [211762,34639475) | 22534 | 0.862745 | -0.770525 | 20 | ALOH | 8 | 0 | 8 | 1 | 0.99649755 | 0.24458808 | 0 | HLAMP |
xxx | 9 | [334016,35884107) | 690 | 0.871795 | -0.68648 | 20 | ALOH | 8 | 0 | 8 | 3 | 0.42151919 | 0.55352253 | 1 | HLAMP |
The second sample is a TCGA sample and I can share it as needed.
Floris
P.s. This issue raises another one in my mind: how good is TITAN at detecting homozygous deletions, as are common in eg. CDKN2A? I wonder because I imagine that the read depth at the normal het site in the matching tumors in this area may not meet the minimum threshold set by TITAN. Any thoughts on this?
Line 50 in 49081fb
Line 358 in 49081fb
Hi,
I've decided to move ahead in hg38 with a smaller bin size of mapability and GC% as the mapability file is much quicker to make at smaller BED intervals. The first TenX script doesn't appear to be pulling in the haplotype.R script. Any advice on what I'm doing wrong?
Rscript getMoleculeCoverage.R
--id TMI40424
-l TitanCNA/scripts/
-d Final_GCandMAP/
-t TMI40424_Tumor.10kb.bxTile.chrX.bed
-n TMI40424_Germline.10kb.bxTile.chrX.bed
--minReadsPerBX=4
--chrs "X"
--outDir TMI40424/
--centromere GRCh38.p13_centromere_UCSC-gapTable.txt
Error:
Loading required package: grid
Error: could not find function "loadBXcountsFromBEDDir"
Execution halted
Best,
Jake
As the subject says, having to hunt a supplementary table in a manuscript to see what an abbreviation indicates is awkward (and, IIRC, there is no "supplementary table 2" in the TitanCNA paper).
Ideally this should be in the README.md or in the repo (or in the help files).
In file https://github.com/gavinha/TitanCNA/blob/master/R/utils.R
in function computeSDbwIndex code:
else if (S_Dbw.method == "Tong"){
###### NOTE #######
## The authors originally used ((N-ni)/N), but this is incorrect ##
## We want a weighted-sum ##
scat.Ci[i] <- ((ni) / N) * (var.Ci/var.D)
Are your sure?
In paper (Youngok Kim and Soowon Lee. A clustering validity assessment Index., 2003,, 602–608):
If two clusters have the same variance but different numbers of tuples, the variance
of cluster with more tuples has more effect on the scatter function than the one with
less tuples. By minimizing the influence of small clusters or noisy clusters, Scat* can
be more robust than Scat of S_Dbw index.
Alexander Lashkov
Hi Gavin,
Thanks for your previous suggestions and I can run TitanCNA using snakemake pipeline successfully now. I have some questions about the parameters in snakemake pipeline:
First is about the parameter: numMaxClonalClusters. Does it mean the max number of clonal clusters to be considered by TitanCNA? In my understanding, there should be only one clonal cluster and all CNAs in that cluster have a cellular prevalence = 1. Or does it just mean the number of clusters to be considered?
I'm not familiar with ichorCNA but I noticed that there are many parameters to be set for it, some of which seems similar to TitanCNA (e.g. ichorCNA_normal/ichorCNA_ploidy and TitanCNA_normalInit/TitanCNA_maxPloidy). I don't know how much influence will the result of ichorCNA have on downstream analysis by TitanCNA. Is the default setting appropriate for my data?
Information for my data:
WXS
depth: 300X
purity: 0.2-0.6
ploidy: 2-4, most samples are 2, about 15% are 3, less than 5% are 4.
(The purity and ploidy are estimated by two other softwares: Sequenza, Facets, and their consistency is good, so I think the estimation should be credible.)
I'm running TitanCNA on a cluster and want to use 20 cores. But our job submission system is different from yours. Our rule is something like:
Run -node_type -num_of_node -num_of_cpus_per_node [command].
So will it work if I set the TitanCNA_numCores in config.yaml to 20 and use the command below?
Run -fat4way -1 -20 snakemake -s TitanCNA.snakefile --cores 20
If I have a SNV with CP calculated by those frequently-used softwares (PyClone, PhyloWGS...) and a CNA with CP calculated by TitanCNA, could I simply compare the CPs to judge their relative order during tumor evolution. For example, CP(SNV) > CP(CNA) means the SNV occurs previous to the CNA? I know it maybe unreasonable because the CP estimated by different softwares might be uncomparable, but I don't know other way to achieve my goal. Do you have some suggestions?
Thank you very much!
Yang
Hi Gavin,
I was wondering whether it was possible to specify the size and orientation of the gene labels in the plotCNlogRByChr() function. I am using the geneAnnot parameter to pass in gene labels. It appears the code calls plotGeneAnnotation, and the values seem to be hardcoded.
Thanks,
Error in data.table(Subclone1 = subc1) : object 'subc1' not found
Calls: outputTitanResults ... eval -> eval -> eval -> getSubcloneProfiles -> data.table
This is likely related to an incompatible solution (i.e. ploidy and number of clusters). I will put in a few checks to allow it to complete.
Would you please suggest how to set some parameters?
I am working on 26 samples from 13 ovarian cancer patients but not sure if my choice of parameters was good enough. Currently, I set the maxPloidy at 4 but I found that 9 out of 26 samples were estimated to have Ploidy at 4. Do I need to increase the max Ploidy greater?
To break my questions down further,
What are your suggestions on these anomalies? Do I need to dump them?
Lastly, I am using "snakemake" as you suggested before and it's fantastic!
Thank you again your answers in advance (please!) and also for developing and maintaining this wonderful application with us!
Bests,
Lijin
Snakemake files are not ready for UCSC chromosome naming. Simple change of CHRS is not enough to succesful execution.
If it is explained in the paper, feel free to close this issue. The definition is clear when you have a single clonal cluster, but when there are more, it gets a little counter-intuitive:
selectSolution
indicates the cellular prevalence in a syntethic way, in order of larger -> smaller (example: 1, 0.86, 0.68). In this case, what does the prevalence signify? The percentage of general alterations in the specific tumor cluster?Hello,
I'm trying to run TitanCNA on some WGS samples. The input data format as per the vignette is
The format of the file must be 6 columns: chromosome, position, reference base, reference read counts, non-reference base, non-reference read counts.
So I've constructed a data frame that looks like
> head(data)
# A tibble: 6 x 6
chr posn ref refOriginal nonRef tumDepth
<chr> <int> <chr> <int> <chr> <int>
1 1 10109 A 172 T 59
2 1 10255 A 36 T 10
3 1 10321 C 8 T 8
4 1 12383 G 11 A 19
5 1 13302 C 14 T 8
6 1 14699 C 25 G 7
Using the naming convention in the vignette. The documentation under ?loadAlleleCounts
also states
The reference and non-reference base columns can be any arbitrary character; it is not used by TitanCNA.
However, when I run this I get the error
Error in data$ref/data$tumDepth : non-numeric argument to binary operator
which traceback()
shows is caused by
cor(data$ref/data$tumDepth, data$logR, use = "pairwise.complete.obs")
ie TitanCNA is trying to take the correlation of ref
which according to the documentation is the "reference base" and isn't used by TitanCNA.
Any help with where I'm going wrong with this analysis much appreciated.
Thanks
Kieran
Hi
I am trying to run snakemake version of pipeline using human_g1k_v37.fasta genome using
snakemake -s TitanCNA.snakefile --cores 5
The run was Okay till this command is executed:
Rscript /home/a/asna2/TitanCNA/scripts/R_scripts/titanCNA.R --hetFile results/titan/tumCounts/tumor_sample_1.tumCounts.txt --cnFile results/ichorCNA/tumor_sample_1/tumor_sample_1.correctedDepth.txt --outFile results/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster1.titan.txt --outSeg results/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster1.segs.txt --outParam results/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster1.params.txt --outIGV results/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster1.seg --outPlotDir results/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster1/ --libdir /home/a/asna2/TitanCNA/ --id tumor_sample_1 --numClusters 1 --numCores 1 --normal_0 0.5 --ploidy_0 3 --genomeStyle NCBI --genomeBuild hg19 --cytobandFile None --chrs "c(1:22, "X")" --estimateNormal map --estimatePloidy True --estimateClonality True --centromere /home/a/asna2/miniconda3/pkgs/r-ichorcna-0.1.0.20180710-r341_0/lib/R/library/ichorCNA/extdata/GRCh37.p13_centromere_UCSC-gapTable.txt --alphaK 10000 --txnExpLen 1e15 --plotYlim "c(-2,4)" > logs/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster1.log 2> logs/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster1.log
returned non-zero exit status 1.
File "/home/a/asna2/TitanCNA/scripts/snakemake/TitanCNA.snakefile", line 62, in __rule_runTitanCNA
File "/home/a/asna2/miniconda3/lib/python3.6/concurrent/futures/thread.py", line 56, in run
The tumor_sample_1_cluster1.log file showed
$dens.bw
[1] 0.379763
$scat
[1] 0.08379496
$S_Dbw
[1] 0.463558
nd in "RSQLite" for requests: dbGetQuery
Running TITAN...
titan: Loading data results/titan/tumCounts/tumor_sample_1.tumCounts.txt
titan: Loading GC content and mappability corrected log2 ratios...
titan: Extracting read depth...
Removed 79 centromeric positions
Removed Chrs:
titan: Loading default parameters
titan: Using 1 cores.
titan: Parameter estimation
Optimal state path computation: Using 1 cores.
Writing results to results/titan/hmm/titanCNA_ploidy2/tumor_sample_1_cluster1.titan.txt, results/titan/hmm/titanCNA_ploidy2/tumor_sample_1_cluster1.segs.txt, results/titan/hmm/titanCNA_ploidy2/tumor_sample_1_cluster1.params.txt
titan: Saving parameters to results/titan/hmm/titanCNA_ploidy2/tumor_sample_1_cluster1.params.txt
Error in cytoband[, "chrom"] : incorrect number of dimensions
Calls: plotIdiogram
In addition: Warning messages:
1: In .replace_seqlevels_style(x_seqlevels, value) :
found more than one best sequence renaming map compatible with seqname style "NCBI" for this object, using the first one
2: In .replace_seqlevels_style(x_seqlevels, value) :
found more than one best sequence renaming map compatible with seqname style "NCBI" for this object, using the first one
3: In .replace_seqlevels_style(x_seqlevels, value) :
found more than one best sequence renaming map compatible with seqname style "NCBI" for this object, using the first one
4: In .replace_seqlevels_style(x_seqlevels, value) :
found more than one best sequence renaming map compatible with seqname style "NCBI" for this object, using the first one
5: In .replace_seqlevels_style(x_seqlevels, value) :
found more than one best sequence renaming map compatible with seqname style "NCBI" for this object, using the first one
Execution halted
I am not sure about what causing the problem but I showed here the the reference settings that have used
reference settings and paths to reference files
genomeBuild: hg19
genomeStyle: NCBI
refFasta: /home/a/asna2/Scrach/WES/human_g1k_v37.fasta
snpVCF: /home/a/asna2/Scrach/WES/hapmap_3.3.b37.vcf.gz
ichorCNA_exons: NULL
cytobandFile: None
centromere: /home/a/asna2/miniconda3/pkgs/r-ichorcna-0.1.0.20180710-r341_0/lib/R/library/ichorCNA/extdata/GRCh37.p13_centromere_UCSC-gapTable.txt
sex: females # use None if both females and males are in sample set
The TitanCNV parameters are
use this for NCBI chr naming
chrs: 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,X,Y
TitanCNA_maxNumClonalClusters: 5
TitanCNA_chrs: c(1:22, "X")
TitanCNA_normalInit: 0.5
TitanCNA_maxPloidy: 3
TitanCNA_estimateNormal: map
TitanCNA_estimatePloidy: TRUE
TitanCNA_estimateClonality: TRUE
TitanCNA_alleleModel: binomial
TitanCNA_alphaK: 10000
TitanCNA_alphaR: 10000
TitanCNA_txnExpLen: 1e15
TitanCNA_plotYlim: c(-2,4)
TitanCNA_solutionThreshold: 0.05
TitanCNA_numCores: 1
TitanCNA_mem: 16G
TitanCNA_runtime: "300:00:00"
TitanCNA_pe: -pe smp 1 -binding linear:1
Any assistance is appreciated.
Hello!
When i trying to run snakemake version of pipeline using hg38 genome, when this command is executed:
Rscript ../R_scripts/titanCNA.R --hetFile results/titan/tumCounts/tumor_sample_1.tumCounts.chrom_fix.txt --cnFile results/ichorCNA/tumor_sample_1/tumor_sample_1.correctedDepth.txt --outFile results/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster1.titan.txt --outSeg results/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster1.segs.txt --outParam results/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster1.params.txt --outIGV results/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster1.seg --outPlotDir results/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster1/ --libdir ../../ --id tumor_sample_1 --numClusters 1 --numCores 1 --normal_0 0.5 --ploidy_0 3 --chrs "c(1:22)" --estimateNormal map --estimatePloidy True --estimateClonality True --centromere /home/vsvekolkin/projects/ichorCNA/inst/extdata/GRCh38.GCA_000001405.2_centromere_acen.txt --alphaK 2500 --txnExpLen 1e15 --plotYlim "c(-2,4)"
Following error is produced:
Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘seqinfo’ for signature ‘"integer"’
Calls: seqlevelsStyle<- ... seqlevels -> seqlevels -> seqlevels -> seqinfo -> <Anonymous>
Execution halted
I believe this is error is caused by parameter --chrs "c(1:22)", which in turn affects GenomeInfoDb , which is used in following line of runTitanCNA script:
seqinfo <- Seqinfo(genome=genomeBuild)
seqlevelsStyle(chrs) <- genomeStyle
Parameter creates integer array, i tried to convert it into character array using: --chrs "paste(c(1:22))"
And i get following error:
Running TITAN...
titan: Loading data results/titan/tumCounts/tumor_sample_1.tumCounts.chrom_fix.txt
titan: Loading GC content and mappability corrected log2 ratios...
titan: Extracting read depth...
Removed 168 centromeric positions
Removed Chrs:
titan: Loading default parameters
titan: Using 1 cores.
titan: Parameter estimation
Error in chrsI[[i]] <- which(data$chr == chrs[i]) :
attempt to select less than one element in integerOneIndex
Calls: runEMclonalCN
Execution halted
I believe this is caused by incorrect filtering in removeCentromere function, but i cannot fix cause of error. Any help is appreciated.
Hi Gavin,
I was able to successfully run TITAN on my BAM files, however, I am trying to run the CNV information through phyloWGS (which says that it is compatible with titan's cnv file format), but it says that my input cnv file should include the following fields:
chromosome
Start_Position
End_Position
MajorCN
MinorCN
I saw in the documentation of the R package that by doing
segs <- outputTitanSegments(results, id = "test", convergeParams, filename = NULL, igvfilename = NULL)
head(segs)
you can get a output with the following sections:
From running the snakemake I found some output files that had the Chr, start, end and copy.number (specifically my tumor_sample_1.seg file in the results/ichorCNA/sample1 folder) however, I can't seem to find them in any of my output files with MinorCN or MajorCN. I've never used R before and was wondering if this was the only way to get this output or if I can get it from the snakemake output?
Thanks,
Paulina
Line 543 in 49081fb
Hi Gavin,
I recently found a little bug to be fixed in selectSolution.R. I run TitanCNA using snakemake. I found some patients disappeared from the final optimalClusterSolution.txt but some patients appeared twice. For example, 07T disappeared but 107T appeared twice. I read the code and found the bug is in this line: phi2Samples <- grep(id, phi2Files, value=T)
If my id is 07T, then grep will catch files of both 07T and 107T, and compare them together. You can imagine that in my case, 107T is always the winner in all conditions, so 07T disappeared from the optimalClusterSolution.txt and 107T appeared twice.
I change this line to:
phi2Samples <- grep(paste('/', id, sep = ''), phi2Files, value=T)
to make sure the id is in the beginning of the filename. It solves my problem. I'm not good at coding and I think you should have a better way to fix this bug in next version of TitanCNA.
Thanks again for making such a convenient pipeline!
Yang
Dear Gavin,
I am trying to run TitanCNA with WGS data generated from a matched tumor/normal pair. After preparing the data and reference using the HHMcopy suite according to the vignette on Bioconductor, I tried to run the function "extractAlleleReadCounts" to generate the allele counts file but I get the error message 'could not find function "extractAlleleReadCounts"'. Do I have to load additional packages? Or do I have to generate the file outside from R as well? I do find the help page for the function but cannot use the function itself. I have installed the package (TitanCNA_1.16.0) using Bioconductor.
Many thanks in advance
Marcel
After numerous changes and tweaks to the default parameter settings over the years, it appears from my experience and from others that TitanCNA tends to overestimate genome-doubling.
I will adjust the default parameters to help correct for this.
Write a Parameter Settings and Guidelines Section in the Wiki
Hi Gavin,
I had a general question about TitanCNA, I saw that there is a maxCN value under 'ichorCNA_maxCN' and this has a default value of 8. For the tumor sample that I'm working with, I know that the CN's can be up to 100 but I am not sure what the maxCN value would be. I saw that the newest version of TitanCNA has the function correctIntegerCN(). Does the new function correctIntegerCN() take care of this issue of not knowing what your maxNC should be and does it simply replace the maxCN when it is higher than inputed value?
I also read that this new function re-calculates the CN by "correcting log ratio based on purity and ploidy and then convert to decimal CN value" and so are the purity and ploidy also values that I should be specifying myself (since the CN value is calculated off of that) or are these also estimated?
Thanks,
Paulina
Hi Gavin,
Your paper that you shared in a previous post (https://www.ncbi.nlm.nih.gov/pubmed/29909985) uses PyClone on top of TITAN, using the TITAN major and minor copy number (I guess) to provide the correct input to PyClone. I'm wondering why did you chose this approach?
For my data, I am interested in comparing clonal clusters across multiple related samples. TITAN I guess is a single-sample approach: "cluster 2" in two related samples may not refer to the same cluster. PyClone on the other hand can take a multi-sample approach, and would be able to track the same clusters between multiple related samples. Is this a sound strategy? Do you have any special suggestions for using TITAN output as PyClone input?
Thanks!
Floris
Line 266 in 49081fb
Hi,
Is there any chance that you have an hg38 mappability file fo rthe 10x Whole Genome Sequencing pipeline (map_hg19_10kb.wig)? I can make them via the gem-mappability tool, but a 10kb length is super inefficient in processing time.
Line 176 in 49081fb
Hi
I want to use TitanCNA and PhyloWGS on my WES data. I noticed that these tools are designed for WGS, so I wonder that have you tried TitanCNA on WES data and does it work well?
Many thanks!
Yang
Hi,
I was able to run TitanCNA succesfully before. Now when I try to run with hg19 aligned samples, it gives error at rule runTitanCNA step;
Error in job runTitanCNA while creating output files ... RuleException: CalledProcessError in line 54 of /TitanCNAmaster/scripts/snakemake/TitanCNA.snakefile:
and the log files say;
Error in fread(cytobandFile) : object 'cytobandFile' not found
It says cytobandFile is only needed if ref.genome is hg38. I don't know how to get around this error. Any help would be much appreciated. Thanks!
Irem
Hi,
I'm currently trying to get the snakemake pipeline up and running, I got to the point where the piepline is attempting to run the rule runTitanCNA:
rule runTitanCNA:
input: results/titan/tumCounts/tumor_sample_1.tumCounts.txt, results/ichorCNA/tumor_sample_1/tumor_sample_1.correctedDepth.txt
output: results/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster2.params.txt, results/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster2.titan.txt, results/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster2/, results/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster2.segs.txt, results/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster2.seg
log: logs/titan/hmm/titanCNA_ploidy3/tumor_sample_1_cluster2.log
jobid: 4
wildcards: clustNum=2, tumor=tumor_sample_1, ploidy=3
The error I get is:
In storage.mode(peek.optstring) = spec[current.flag, col.mode] :
NAs introduced by coercion
Running TITAN...
titan: Loading data results/titan/tumCounts/tumor_sample_1.tumCounts.txt
titan: Loading GC content and mappability corrected log2 ratios...
titan: Extracting read depth...
Removed 46 centromeric positions
Removed Chrs:
titan: Loading default parameters
Model selection: Using 1 cores.
Error in uniroot(funcP, intervalPhi, tol = 1e-15) :
f.lower = f(lower) is NA
Calls: runEMclonalCN ... estimateClonalCNParamsMap -> updateParameters -> uniroot
Execution halted
I've noticed there are some NAs in the file results/ichorCNA/tumor_sample_1/tumor_sample_1.correctedDepth.txt. see below. Could this be the problem?
chr start end log2_TNratio_corrected
1 1000001 2000000 NA
1 2000001 3000000 NA
1 4000001 5000000 -0.115543015089316
1 5000001 6000000 -0.0223702614590051
1 6000001 7000000 -0.106587597655531
1 7000001 8000000 0.0909244472332157
1 8000001 9000000 0.141146096230084
1 9000001 10000000 -0.0804800930951752
Any ideas on what's going wrong?
Thanks
Rob
I'm running TitanCNA v1.4.0 and have the following issue:
>
> convergeParams <- runEMclonalCN(workdata,gParams=params$genotypeParams,nParams=params$normalParams,
+ pParams=params$ploidyParams,sParams=params$cellPrevParams,
+ maxiter=20,maxiterUpdate=1500,txnExpLen=1e15,txnZstrength=1e5,
+ useOutlierState=FALSE,
+ normalEstimateMethod="map",estimateS=TRUE,estimatePloidy=TRUE)
titan: Running HMM...
fwdBack: Iteration 1 chr: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 24 *** Error in `/usr/local/R/R-3.1.1/lib/R/bin/exec/R': free(): invalid next size (fast): 0x0000000003a5c8a0 ***
How can I solve this issue? I'm working with UCSC style chromosome names.
Gavin and all;
I'm hitting an issue with build 37 test runs and not sure about the cause. These are small runs that I don't expect real results on, so I'm not sure if it's due to something with that are being build 37. Using the latest development script it errors out with object 'genes' not found
from plotCNlogRByChr
plotting:
TitanCNA/scripts/R_scripts/titanCNA.R
Line 346 in eb288ea
I'm not sure where genes
gets injected into the environment and available. My guess is it comes from some import I'm missing or is out of date. Do you have any pointers so I can try to track down further?
Thanks so much for any ideas.
set -o pipefail; unset R_HOME && unset R_LIBS && export PATH=/usr/local/share/bcbio/anaconda/bin:$PATH && titanCNA.R --id c-tumor --hetFile /home/chapmanb/bio/bcbio-nextgen/tests/test_automated_output/structural/c-tumor/titancna/PairedBatch-effects-annotated-damage-ann-vardict-prep.csv --cnFile /home/chapmanb/bio/bcbio-nextgen/tests/test_automated_output/structural/c-tumor/titancna/c-tumor-normalized.cn --numClusters 1 --ploidy 2 --numCores 1 --outDir /home/chapmanb/bio/bcbio-nextgen/tests/test_automated_output/bcbiotx/tmprS3mKs --libdir None --chrs "c('22','X')" --genomeBuild hg19 --genomeStyle UCSC
There were 11 warnings (use warnings() to see them)
Attaching package: 'Biostrings'
The following object is masked from 'package:base':
strsplit
Running TITAN...
titan: Loading data /home/chapmanb/bio/bcbio-nextgen/tests/test_automated_output/structural/c-tumor/titancna/PairedBatch-effects-annotated-damage-ann-vardict-prep.csv
titan: Loading GC content and mappability corrected log2 ratios...
titan: Extracting read depth...
Removed Chrs:
titan: Loading default parameters
Warning message:
In cor(data$ref/data$tumDepth, data$logR, use = "pairwise.complete.obs") :
the standard deviation is zero
titan: Using 1 cores.
titan: Parameter estimation
Optimal state path computation: Using 1 cores.
Writing results to /home/chapmanb/bio/bcbio-nextgen/tests/test_automated_output/bcbiotx/tmprS3mKs/c-tumor_cluster01.titan.txt, /home/chapmanb/bio/bcbio-nextgen/tests/test_automated_output/bcbiotx/tmprS3mKs/c-tumor_cluster01.segs.txt, /home/chapmanb/bio/bcbio-nextgen/tests/test_automated_output/bcbiotx/tmprS3mKs/c-tumor_cluster01.params.txt
titan: Saving parameters to /home/chapmanb/bio/bcbio-nextgen/tests/test_automated_output/bcbiotx/tmprS3mKs/c-tumor_cluster01.params.txt
$dens.bw
[1] NaN
$scat
[1] NaN
$S_Dbw
[1] NaN
Error in plotCNlogRByChr(dataIn = results, chr = chrs, segs = segs, ploidy = ploidy, :
object 'genes' not found
Hi,
I am trying to run titanCNA on matched tumor-normal WES samples. I am following snakemake workflow but I am also running it inside of a HPC cloud. I downloaded both titanCNA and ichorCNA.
-My first question is does the parameters of both config files (titan and ichor) need to be matched? I know titan snakefile invokes the others but in titan config file there is still a path for ichor script and libdir so I thought it migth affect.
-When I changed the binsize from 10.000 to 100.000 or 1.000.000 the workflow gives error on earlier stages of the jobs. When I keep it at 10.00 I still get an error at runTitanCNA which is:
Error in job runTitanCNA while creating output files results/titan/hmm/titanCNA_ploidy2/RTF8_D_cluster2/, results/titan/hmm/titanCNA_ploidy2/RTF8_D_cluster2.titan.txt, results/titan/hmm/titanCNA_ploidy2/RTF8_D_cluster2.params.txt, results/titan/hmm/titanCNA_ploidy2/RTF8_D_cluster2.segs.txt, results/titan/hmm/titanCNA_ploidy2/RTF8_D_cluster2.seg.
RuleException:
CalledProcessError in line 53 of /mnt/lustre/sarihae/HNSCC_WES/titanCNA/scripts/snakemake/TitanCNA.snakefile: returned non-zero exit status 1.
File "/mnt/lustre/sarihae/HNSCC_WES/titanCNA/scripts/snakemake/TitanCNA.snakefile", line 53, in __rule_runTitanCNA
File "/cm/shared/apps/python/3.6.2/lib/python3.6/concurrent/futures/thread.py", line 55, in run
Removing output files of failed job runTitanCNA since they might be corrupted:
results/titan/hmm/titanCNA_ploidy2/RTF8_D_cluster2/
Will exit after finishing currently running jobs.
Exiting because a job execution failed. Look above for error message
I looked at the config file and changed libdirs from path/to/titanCNA/R/ path/to/ichorCNA-master/R/
to path/to/titanCNA/ path/to/ichorCNA-master/
and now I get an error at selectSolution which looks like:
Error in job selectSolution while creating output file results/titan/hmm/optimalClusterSolution.txt.
RuleException:
CalledProcessError in line 70 of /mnt/lustre/sarihae/HNSCC_WES/titanCNA/scripts/snakemake/TitanCNA.snakefile: returned non-zero exit status 2.
File "/mnt/lustre/sarihae/HNSCC_WES/titanCNA/scripts/snakemake/TitanCNA.snakefile", line 70, in __rule_selectSolution
File "/cm/shared/apps/python/3.6.2/lib/python3.6/concurrent/futures/thread.py", line 55, in run
Will exit after finishing currently running jobs.
Exiting because a job execution failed. Look above for error message
What would you think it's causing my problem, is it related to input parameters, or sth else? I'd really appreciate any kind of help and guidance cause I am a little bit confused with this workflow.
Thanks in advance,
Irem
snakemake -s ichorCNA.snakefile -np
Building DAG of jobs...
ChildIOException:
File/directory is a child to another output:
/spin1/home/TitanCNA/scripts/snakemake/results/ichorCNA/tumor_X
/spin1/home/TitanCNA/scripts/snakemake/results/ichorCNA/tumor_X/tumor_X.cna.seg
As per error #26 , switching back from snakemake 5.3 to 3.8 seemed to resolve the problem, but I wanted to make the problem visible.
In my case, I had 4 clonal clusters. In the code for the function, clusters
expands to
> clusters
ClonalCluster CellularPrevalence
1: 2 0.7714099
2: 3 0.6370672
3: 1 1.0000000
4: 4 0.5077241
When plotting, the code does
if (nrow(clusters) > 0) {
for (j in 1:length(clusters[, 1])) {
However:
> 1:length(clusters[, 1])
[1] 1
So only whatever is in the first line (and not necessarily the most abundant cluster) is plotted on the graph. I will make a PR rewriting the part using iteration on the whole data.table (which is small).
The main TitanCNA script generates file with param
in their names for parameters. However, selectSolution
looks (using list.files
) for files using params
as pattern, and therefore does not run to completion because all file lists are empty.
Hello Gavin,
I am following your tutorial but got into a trouble from the very beginning step.
Would you mind to take a closer look at "getPositionOverlap" function?
I installed this package a few weeks ago.
Thank you for your help.
Lijin
logR <- getPositionOverlap(chr=data$chr, posn=data$posn, dataVal=cnData)
head(logR)
[1] NA NA NA NA NA NA
tail(logR)
[1] NA NA NA NA NA NA
As subject says.
Hello,
I need to download TitanCNA 1.5.8 to use it with Titanrunner pipeline. I couldn't find a link to download this version. Is it available for download?
Thank you
Correcting reads with targeted regions failed because I had genomeStyle
set to NCBI
while my data was in UCSC
style, therefore no overlaps would be found.
Hi, I'm having trouble having the snakemake output results. I'm attempting to run the command "snakemake -s ichorCNA.snakefile -np" and get the following output (see below). I unfortunately can't tell clearly the error from the output, as the output does not explicitly throw an exception or state which process/dependency might be improperly linked to in config.yaml. It would seem that the output of the command would be in the "results" directory, but there is none to be found.
I'm using the following versions of software:
TitanCNA - 1.17.0 (I believe the current version in the Titan snakemake repo)
ichorCNA - 0.1.0
HMMcopy - 1.24.0
optparse - 1.6.0
SNPchip - 2.28.0
stringr - 1.3.1
Python-3.5.2
snakemake - 3.12.0
PySAM-0.15.1
PyYAML-3.13
samtools-1.3.1
HMM_Copy_Utils-1.0
rule read_counter:
input: /home/ubuntu/Sample_PM1043_Z11_1_p2_Case_HALO.md.bam
output: results/readDepth/tumor_sample_1.bin10000.wig
log: logs/readDepth/tumor_sample_1.bin10000.log
jobid: 1
wildcards: samples=tumor_sample_1, binSize=10000
resources: mem=4
/home/ubuntu/hmmcopy_utils/bin/readCounter /home/ubuntu/Sample_PM1043_Z11_1_p2_Case_HALO.md.bam -c 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,X,Y -w 10000 -q 20 > results/readDepth/tumor_sample_1.bin10000.wig 2> logs/readDepth/tumor_sample_1.bin10000.log
rule read_counter:
input: /home/ubuntu/Sample_PM1043_EBC4_2_p2_Ctrl_HALO.md.bam
output: results/readDepth/normal_sample_1.bin10000.wig
log: logs/readDepth/normal_sample_1.bin10000.log
jobid: 3
wildcards: samples=normal_sample_1, binSize=10000
resources: mem=4
/home/ubuntu/hmmcopy_utils/bin/readCounter /home/ubuntu/Sample_PM1043_EBC4_2_p2_Ctrl_HALO.md.bam -c 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,X,Y -w 10000 -q 20 > results/readDepth/normal_sample_1.bin10000.wig 2> logs/readDepth/normal_sample_1.bin10000.log
rule ichorCNA:
input: results/readDepth/tumor_sample_1.bin10000.wig, results/readDepth/normal_sample_1.bin10000.wig
output: results/ichorCNA/tumor_sample_1/tumor_sample_1.correctedDepth.txt, results/ichorCNA/tumor_sample_1/, results/ichorCNA/tumor_sample_1/tumor_sample_1.cna.seg
log: logs/ichorCNA/tumor_sample_1.log
jobid: 2
wildcards: tumor=tumor_sample_1
resources: mem=4
Rscript /home/ubuntu/ichorCNA/scripts/runIchorCNA.R --libdir /home/ubuntu/ichorCNA/R/ --id tumor_sample_1 --WIG results/readDepth/tumor_sample_1.bin10000.wig --gcWig /home/ubuntu/ichorCNA/inst/extdata/gc_hg19_10kb.wig --mapWig /home/ubuntu/ichorCNA/inst/extdata/map_hg19_10kb.wig --NORMWIG results/readDepth/normal_sample_1.bin10000.wig --ploidy "c(2,3)" --normal "c(0.5)" --maxCN 8 --includeHOMD False --genomeStyle NCBI --chrs "c(1:22, "X")" --estimateNormal True --estimatePloidy True --estimateScPrevalence True --scStates "c(1,3)" --centromere /home/ubuntu/ichorCNA/inst/extdata/GRCh37.p13_centromere_UCSC-gapTable.txt --exons.bed None --txnE 0.9999 --txnStrength 10000 --fracReadsInChrYForMale 0.001 --plotFileType png --plotYLim "c(-2,4)" --outDir results/ichorCNA/tumor_sample_1/ > logs/ichorCNA/tumor_sample_1.log 2> logs/ichorCNA/tumor_sample_1.log
localrule correctDepth:
input: results/ichorCNA/tumor_sample_1/tumor_sample_1.cna.seg, results/ichorCNA/tumor_sample_1/tumor_sample_1.correctedDepth.txt, results/readDepth/tumor_sample_1.bin10000.wig, results/readDepth/normal_sample_1.bin10000.wig
jobid: 0
Job counts:
count jobs
1 correctDepth
1 ichorCNA
2 read_counter
4
Any help would be appreciated! Thank you!
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.