search

raphael-group / hatchet

HATCHet (Holistic Allele-specific Tumor Copy-number Heterogeneity) is an algorithm that infers allele and clone-specific CNAs and WGDs jointly across multiple tumor samples from the same patient, and that leverages the relationships between clones in these samples.
BSD 3-Clause "New" or "Revised" License
66 stars 31 forks source link

Issue running count_alleles #150

Open wir963 opened 2 years ago

wir963 commented 2 years ago

Hey,

I'm trying to convert from running HATCHet v0 to HATCHet v1 and I'm running into an issue with the count_alleles step. Here's the traceback message. I'm using hatchet 1.0.1 downloaded via bioconda. 

hatchet count-alleles 
--tumors tumor_bams
--normal normal_bams
--reference REFERENCE_FASTQ_FILE
--snps {params.snps} 
--outputnormal {output[0]} output/4214N_PBL_November_15_2016/normal.1bed
--outputtumors {output[1]} output/4214N_PBL_November_15_2016/tumor.1bed
--outputsnps output/germline_heterozygous_SNPs
--samples 4214N_PBL_November_15_2016 4214-1Met_FrTu_November_15_2016 "
--mincov 8
--maxcov 300
--processes 8 
ESC[0mESC[1mESC[95m[2022-Jul-18 11:56:24]# Counting SNPs alleles from the matched-normal sample
ESC[0mESC[2K^MESC[96mProgress:ESC[0m |----------------------------------------| ESC[96m0.0% CompleteESC[0m ESC[96m[[2022-Jul-18 11:56:24]AlleleCounter-1 starts on 4214N_PBL_November_15_2016 for 1]ESC[0mProcess AlleleCounter-1:
ESC[2K^MESC[96mProgress:ESC[0m |----------------------------------------| ESC[96m0.0% CompleteESC[0m ESC[96m[[2022-Jul-18 11:56:24]AlleleCounter-2 starts on 4214N_PBL_November_15_2016 for 2]ESC[0mProcess AlleleCounter-2:
ESC[2K^MESC[96mProgress:ESC[0m |----------------------------------------| ESC[96m0.0% CompleteESC[0m ESC[96m[[2022-Jul-18 11:56:24]AlleleCounter-3 starts on 4214N_PBL_November_15_2016 for 3]ESC[0mProcess AlleleCounter-3:
Traceback (most recent call last):
Traceback (most recent call last):
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/84674ea22fe618b33f6df301957543a6/lib/python3.9/multiprocessing/process.py", line 315, in _bootstrap
    self.run()
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/84674ea22fe618b33f6df301957543a6/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 291, in run
    snps = self.countAlleles(
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/84674ea22fe618b33f6df301957543a6/lib/python3.9/multiprocessing/process.py", line 315, in _bootstrap
    self.run()
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/84674ea22fe618b33f6df301957543a6/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 321, in countAlleles
    with open(errname, 'w') as err:
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/84674ea22fe618b33f6df301957543a6/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 291, in run
    snps = self.countAlleles(
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/84674ea22fe618b33f6df301957543a6/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 321, in countAlleles
    with open(errname, 'w') as err:
FileNotFoundError: [Errno 2] No such file or directory: 'output/germline_heterozygous_SNPs/4214N_PBL_November_15_2016_1_bcftools.log'
FileNotFoundError: [Errno 2] No such file or directory: 'output/germline_heterozygous_SNPs/4214N_PBL_November_15_2016_1_bcftools.log'

It also doesn't terminate but it just keeps doing nothing as far as I can tell (I eventually let it time out after 4 days).

I'm not sure at all what's going on and would appreciate some help. I was able to get HATCHet v0 to work quite easily.

Best, Welles

vineetbansal commented 2 years ago

It looks like the folder output/germline_heterozygous_SNPs may not exist in your case. The count_alleles step will try to run bcftools for allele counting using this folder as a working directory of sorts. Pre-creating this folder if it doesn't exist as a snakemake step should fix this issue.

This behavior in this step may have changed from the previous HATCHet versions.

wir963 commented 2 years ago

Hey @vineetbansal,

That worked - thanks for the help! One note - after the above error, the job kept running and did not exit (from my perspective, I think the correct behavior would have been to exit).

Best, Welles

wir963 commented 2 years ago

Hey @vineetbansal ,

I'm still getting the weird behavior where the count-alleles step just keep going indefinitely. After fixing the issue above, now it just keeps running (it's been running for over two days now) but the last progress message was from 2022-Jul-25 12:26:15. I've added the last part of the log message at the end.

Maybe it just takes three plus days to finish but it's weird that it hasn't logged anything in almost two days.

Best, Welles

ESC[0mESC[1mESC[95m[2022-Jul-25 12:05:09]# Counting SNPs alleles from the matched-normal sample
ESC[0mESC[2K^MESC[96mProgress:ESC[0m |----------------------------------------| ESC[96m0.0% CompleteESC[0m ESC[96m[[2022-Jul-25 12:05:09]AlleleCounter-1 starts on 4423N_PBL_April_25_2022 for 1]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |----------------------------------------| ESC[96m0.0% CompleteESC[0m ESC[96m[[2022-Jul-25 12:05:09]AlleleCounter-2 starts on 4423N_PBL_April_25_2022 for 2]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |----------------------------------------| ESC[96m0.0% CompleteESC[0m ESC[96m[[2022-Jul-25 12:05:09]AlleleCounter-3 starts on 4423N_PBL_April_25_2022 for 3]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |----------------------------------------| ESC[96m0.0% CompleteESC[0m ESC[96m[[2022-Jul-25 12:05:09]AlleleCounter-4 starts on 4423N_PBL_April_25_2022 for 4]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |----------------------------------------| ESC[96m0.0% CompleteESC[0m ESC[96m[[2022-Jul-25 12:05:09]AlleleCounter-5 starts on 4423N_PBL_April_25_2022 for 5]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |----------------------------------------| ESC[96m0.0% CompleteESC[0m ESC[96m[[2022-Jul-25 12:05:09]AlleleCounter-6 starts on 4423N_PBL_April_25_2022 for 6]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |----------------------------------------| ESC[96m0.0% CompleteESC[0m ESC[96m[[2022-Jul-25 12:05:09]AlleleCounter-7 starts on 4423N_PBL_April_25_2022 for 7]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |----------------------------------------| ESC[96m0.0% CompleteESC[0m ESC[96m[[2022-Jul-25 12:05:09]AlleleCounter-8 starts on 4423N_PBL_April_25_2022 for 8]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |█---------------------------------------| ESC[96m4.2% CompleteESC[0m ESC[96m[[2022-Jul-25 12:11:03]AlleleCounter-2 ends on 4423N_PBL_April_25_2022 for 2]
ESC[0mESC[2K^MESC[96mProgress:ESC[0m |█---------------------------------------| ESC[96m4.2% CompleteESC[0m ESC[96m[[2022-Jul-25 12:11:03]AlleleCounter-2 starts on 4423N_PBL_April_25_2022 for 9]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |███-------------------------------------| ESC[96m8.3% CompleteESC[0m ESC[96m[[2022-Jul-25 12:11:16]AlleleCounter-6 ends on 4423N_PBL_April_25_2022 for 6]
ESC[0mESC[2K^MESC[96mProgress:ESC[0m |███-------------------------------------| ESC[96m8.3% CompleteESC[0m ESC[96m[[2022-Jul-25 12:11:16]AlleleCounter-6 starts on 4423N_PBL_April_25_2022 for 10]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |█████-----------------------------------| ESC[96m12.5% CompleteESC[0m ESC[96m[[2022-Jul-25 12:11:52]AlleleCounter-8 ends on 4423N_PBL_April_25_2022 for 8]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |█████-----------------------------------| ESC[96m12.5% CompleteESC[0m ESC[96m[[2022-Jul-25 12:11:52]AlleleCounter-8 starts on 4423N_PBL_April_25_2022 for 11]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |██████----------------------------------| ESC[96m16.7% CompleteESC[0m ESC[96m[[2022-Jul-25 12:13:41]AlleleCounter-5 ends on 4423N_PBL_April_25_2022 for 5]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |██████----------------------------------| ESC[96m16.7% CompleteESC[0m ESC[96m[[2022-Jul-25 12:13:41]AlleleCounter-5 starts on 4423N_PBL_April_25_2022 for 12]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |████████--------------------------------| ESC[96m20.8% CompleteESC[0m ESC[96m[[2022-Jul-25 12:13:41]AlleleCounter-1 ends on 4423N_PBL_April_25_2022 for 1]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |████████--------------------------------| ESC[96m20.8% CompleteESC[0m ESC[96m[[2022-Jul-25 12:13:41]AlleleCounter-1 starts on 4423N_PBL_April_25_2022 for 13]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |██████████------------------------------| ESC[96m25.0% CompleteESC[0m ESC[96m[[2022-Jul-25 12:14:07]AlleleCounter-3 ends on 4423N_PBL_April_25_2022 for 3]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |██████████------------------------------| ESC[96m25.0% CompleteESC[0m ESC[96m[[2022-Jul-25 12:14:07]AlleleCounter-3 starts on 4423N_PBL_April_25_2022 for 14]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |███████████-----------------------------| ESC[96m29.2% CompleteESC[0m ESC[96m[[2022-Jul-25 12:14:10]AlleleCounter-7 ends on 4423N_PBL_April_25_2022 for 7]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |███████████-----------------------------| ESC[96m29.2% CompleteESC[0m ESC[96m[[2022-Jul-25 12:14:10]AlleleCounter-7 starts on 4423N_PBL_April_25_2022 for 15]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |█████████████---------------------------| ESC[96m33.3% CompleteESC[0m ESC[96m[[2022-Jul-25 12:14:35]AlleleCounter-4 ends on 4423N_PBL_April_25_2022 for 4]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |█████████████---------------------------| ESC[96m33.3% CompleteESC[0m ESC[96m[[2022-Jul-25 12:14:35]AlleleCounter-4 starts on 4423N_PBL_April_25_2022 for 16]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |███████████████-------------------------| ESC[96m37.5% CompleteESC[0m ESC[96m[[2022-Jul-25 12:17:19]AlleleCounter-2 ends on 4423N_PBL_April_25_2022 for 9]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |███████████████-------------------------| ESC[96m37.5% CompleteESC[0m ESC[96m[[2022-Jul-25 12:17:19]AlleleCounter-2 starts on 4423N_PBL_April_25_2022 for 17]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |████████████████------------------------| ESC[96m41.7% CompleteESC[0m ESC[96m[[2022-Jul-25 12:17:35]AlleleCounter-6 ends on 4423N_PBL_April_25_2022 for 10]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |████████████████------------------------| ESC[96m41.7% CompleteESC[0m ESC[96m[[2022-Jul-25 12:17:35]AlleleCounter-6 starts on 4423N_PBL_April_25_2022 for 18]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |██████████████████----------------------| ESC[96m45.8% CompleteESC[0m ESC[96m[[2022-Jul-25 12:18:18]AlleleCounter-8 ends on 4423N_PBL_April_25_2022 for 11]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |██████████████████----------------------| ESC[96m45.8% CompleteESC[0m ESC[96m[[2022-Jul-25 12:18:18]AlleleCounter-8 starts on 4423N_PBL_April_25_2022 for 19]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |████████████████████--------------------| ESC[96m50.0% CompleteESC[0m ESC[96m[[2022-Jul-25 12:20:22]AlleleCounter-5 ends on 4423N_PBL_April_25_2022 for 12]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |████████████████████--------------------| ESC[96m50.0% CompleteESC[0m ESC[96m[[2022-Jul-25 12:20:22]AlleleCounter-5 starts on 4423N_PBL_April_25_2022 for 20]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |█████████████████████-------------------| ESC[96m54.2% CompleteESC[0m ESC[96m[[2022-Jul-25 12:20:36]AlleleCounter-1 ends on 4423N_PBL_April_25_2022 for 13]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |█████████████████████-------------------| ESC[96m54.2% CompleteESC[0m ESC[96m[[2022-Jul-25 12:20:36]AlleleCounter-1 starts on 4423N_PBL_April_25_2022 for 21]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |███████████████████████-----------------| ESC[96m58.3% CompleteESC[0m ESC[96m[[2022-Jul-25 12:20:57]AlleleCounter-3 ends on 4423N_PBL_April_25_2022 for 14]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |███████████████████████-----------------| ESC[96m58.3% CompleteESC[0m ESC[96m[[2022-Jul-25 12:20:57]AlleleCounter-3 starts on 4423N_PBL_April_25_2022 for 22]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |█████████████████████████---------------| ESC[96m62.5% CompleteESC[0m ESC[96m[[2022-Jul-25 12:21:25]AlleleCounter-7 ends on 4423N_PBL_April_25_2022 for 15]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |█████████████████████████---------------| ESC[96m62.5% CompleteESC[0m ESC[96m[[2022-Jul-25 12:21:25]AlleleCounter-7 starts on 4423N_PBL_April_25_2022 for X]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |██████████████████████████--------------| ESC[96m66.7% CompleteESC[0m ESC[96m[[2022-Jul-25 12:23:01]AlleleCounter-4 ends on 4423N_PBL_April_25_2022 for 16]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |██████████████████████████--------------| ESC[96m66.7% CompleteESC[0m ESC[96m[[2022-Jul-25 12:23:01]AlleleCounter-4 starts on 4423N_PBL_April_25_2022 for Y]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |████████████████████████████------------| ESC[96m70.8% CompleteESC[0m ESC[96m[[2022-Jul-25 12:23:38]AlleleCounter-2 ends on 4423N_PBL_April_25_2022 for 17]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |██████████████████████████████----------| ESC[96m75.0% CompleteESC[0m ESC[96m[[2022-Jul-25 12:23:52]AlleleCounter-6 ends on 4423N_PBL_April_25_2022 for 18]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |███████████████████████████████---------| ESC[96m79.2% CompleteESC[0m ESC[96m[[2022-Jul-25 12:24:33]AlleleCounter-8 ends on 4423N_PBL_April_25_2022 for 19]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |█████████████████████████████████-------| ESC[96m83.3% CompleteESC[0m ESC[96m[[2022-Jul-25 12:25:54]AlleleCounter-5 ends on 4423N_PBL_April_25_2022 for 20]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |███████████████████████████████████-----| ESC[96m87.5% CompleteESC[0m ESC[96m[[2022-Jul-25 12:26:02]AlleleCounter-1 ends on 4423N_PBL_April_25_2022 for 21]ESC[0mESC[2K^MESC[96mProgress:ESC[0m |████████████████████████████████████----| ESC[96m91.7% CompleteESC[0m ESC[96m[[2022-Jul-25 12:26:15]AlleleCounter-3 ends on 4423N_PBL_April_25_2022 for 22]ESC[0m
wir963 commented 2 years ago

Hey @vineetbansal ,

I'm bumping this thread because I'd like to figure out what's happening with the count-alleles step in my data so we can run the updated version of HATCHet.

Best, Welles

vineetbansal commented 2 years ago

Hi @wir963 - so I guess your progress bar still doesn't go forward? One of the issues with doing multiprocessing in HATCHet is that if something goes wrong, it has the potential to stall progress (like you saw in your original case with the missing directory). Later today we're releasing a 1.0.2 version with slightly modified code for count_alleles that will make any exceptions in the step obvious (snakemake should error out with an appropriate message).

So while this won't solve your problem, it will hopefully make it more obvious as to what's going on, so we can proceed further. If you can try your workflow with this new version after it's available on bioconda later today, you should get these changes.

wir963 commented 2 years ago

Hey @vineetbansal ,

You were indeed correct. I just pulled down and re-ran my code using hatchet 1.0.2 and this time it fails (which is progress).

Here's the error message - any idea what happened?

Traceback (most recent call last):
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/80845879bd9da5c9c796cfbbd23e868d/bin/hatchet", line 33, in <module>
    sys.exit(load_entry_point('hatchet==1.0.2', 'console_scripts', 'hatchet')())
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/80845879bd9da5c9c796cfbbd23e868d/lib/python3.9/site-packages/hatchet/__main__.py", line 61, in main
    globals()[command](args)
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/80845879bd9da5c9c796cfbbd23e868d/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 29, in main
    snps = counting(
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/80845879bd9da5c9c796cfbbd23e868d/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 193, in counting
    results = worker.run(work=work, n_instances=num_workers)
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/80845879bd9da5c9c796cfbbd23e868d/lib/python3.9/site-packages/hatchet/utils/multiprocessing.py", line 97, in run
    raise handler_result.__class__(f'HANDLER {handler_id} FAILED\n\n{error_string}')
RuntimeError: HANDLER 2 FAILED

Traceback (most recent call last):
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/80845879bd9da5c9c796cfbbd23e868d/lib/python3.9/site-packages/hatchet/utils/multiprocessing.py", line 31, in run
    result = self.worker.work(*args)
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/80845879bd9da5c9c796cfbbd23e868d/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 225, in work
    snps = self.countAlleles(bamfile=bamfile, samplename=samplename, chromosome=chromosome)
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/80845879bd9da5c9c796cfbbd23e868d/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 263, in countAlleles
    error('Allele counting failed on {} of {}, please check errors in {}!').format(
AttributeError: 'NoneType' object has no attribute 'format'

Best, Welles

alvinwt commented 2 years ago

Hi @vineetbansal, I have been having problems with count_alleles as well and I think there is something wrong with the workers writing out the bcftools output. I ran it with 12 parallel processes and it seems like only the last 12 files had output written to them. I have tried this with 8 processes and it had 8 files written out as well. I'm not familiar with how the parallel processing is carried out with bcftools and perhaps that might be the issue.

https://github.com/raphael-group/hatchet/blob/9832b016db07e5a6cbfa5234bd28c4b6ebdcc90d/src/hatchet/utils/count_alleles.py#L59 Best wishes, Alvin

-rw-r--r-- 1 ng02 stlab    0 Aug  8 14:54 TMP_chr10.tsv
-rw-r--r-- 1 ng02 stlab    0 Aug  8 14:54 TMP_chr11.tsv
-rw-r--r-- 1 ng02 stlab    0 Aug  8 14:54 TMP_chr12.tsv
-rw-r--r-- 1 ng02 stlab  76K Aug  8 14:54 TMP_chr13.tsv
-rw-r--r-- 1 ng02 stlab  70K Aug  8 14:54 TMP_chr14.tsv
-rw-r--r-- 1 ng02 stlab  51K Aug  8 14:54 TMP_chr15.tsv
-rw-r--r-- 1 ng02 stlab  67K Aug  8 14:54 TMP_chr16.tsv
-rw-r--r-- 1 ng02 stlab  64K Aug  8 14:54 TMP_chr17.tsv
-rw-r--r-- 1 ng02 stlab  60K Aug  8 14:54 TMP_chr18.tsv
-rw-r--r-- 1 ng02 stlab  51K Aug  8 14:54 TMP_chr19.tsv
-rw-r--r-- 1 ng02 stlab    0 Aug  8 14:54 TMP_chr1.tsv
-rw-r--r-- 1 ng02 stlab  49K Aug  8 14:54 TMP_chr20.tsv
-rw-r--r-- 1 ng02 stlab  29K Aug  8 14:54 TMP_chr21.tsv
-rw-r--r-- 1 ng02 stlab  33K Aug  8 14:54 TMP_chr22.tsv
-rw-r--r-- 1 ng02 stlab    0 Aug  8 14:54 TMP_chr2.tsv
-rw-r--r-- 1 ng02 stlab    0 Aug  8 14:54 TMP_chr3.tsv
-rw-r--r-- 1 ng02 stlab    0 Aug  8 14:54 TMP_chr4.tsv
-rw-r--r-- 1 ng02 stlab    0 Aug  8 14:54 TMP_chr5.tsv
-rw-r--r-- 1 ng02 stlab    0 Aug  8 14:54 TMP_chr6.tsv
-rw-r--r-- 1 ng02 stlab    0 Aug  8 14:54 TMP_chr7.tsv
-rw-r--r-- 1 ng02 stlab    0 Aug  8 14:54 TMP_chr8.tsv
-rw-r--r-- 1 ng02 stlab    0 Aug  8 14:54 TMP_chr9.tsv
-rw-r--r-- 1 ng02 stlab 6.1K Aug  8 14:54 TMP_chrX.tsv
-rw-r--r-- 1 ng02 stlab 1.3K Aug  8 14:54 TMP_chrY.tsv
wir963 commented 2 years ago

Hey @vineetbansal,

I tried to adapt the previous demo to the updated version (and incorporating snakemake) and I also ran into a different issue with count-alleles. The codebase is here (https://github.com/wir963/HATCHet-demo) and the error message is below. 

Building DAG of jobs...
Using shell: /usr/bin/bash
Provided cores: 4
Rules claiming more threads will be scaled down.
Select jobs to execute...

[Fri Aug 12 19:46:46 2022]
rule hatchet_count_alleles:
    input: data/normal.bam, data/bulk_03clone1_06clone0_01normal.sorted.bam, data/bulk_08clone1_Noneclone0_02normal.sorted.bam, data/bulk_Noneclone1_09clone0_01normal.sorted.bam
    output: new_output/baf/normal.1bed, new_output/baf/tumor.1bed, new_output/count_alleles/snps
    jobid: 0
    resources: mem_mb=7451, disk_mb=7451, tmpdir=/tmp

Activating conda environment: .snakemake/conda/62886ca61eec4c5786f8f5bdb83ae223
[2022-Aug-12 19:47:16]# Parsing the input arguments, checking the consistency of given files, and extracting required information

    normal: ('data/normal.bam', 'normal')
    samples: {('data/bulk_Noneclone1_09clone0_01normal.sorted.bam', 'tumor3'), ('data/bulk_03clone1_06clone0_01normal.sorted.bam', 'tumor1'), ('data/bulk_08clone1_Noneclone0_02normal.sorted.bam', 'tumor2')}
    chromosomes: ['chr1', 'chr2', 'chr3', 'chr4', 'chr5', 'chr6', 'chr7', 'chr8', 'chr9', 'chr10', 'chr11', 'chr12', 'chr13', 'chr14', 'chr15', 'chr16', 'chr17', 'chr18', 'chr19', 'chr20', 'chr21', 'chr22', 'chrX', 'chrY']
    samtools: samtools
    bcftools: bcftools
    snps: {'chr1': 'new_output/snps/chr1.vcf.gz', 'chr2': 'new_output/snps/chr2.vcf.gz', 'chr3': 'new_output/snps/chr3.vcf.gz', 'chr4': 'new_output/snps/chr4.vcf.gz', 'chr5': 'new_output/snps/chr5.vcf.gz', 'chr6': 'new_output/snps/chr6.vcf.gz', 'chr7': 'new_output/snps/chr7.vcf.gz', 'chr8': 'new_output/snps/chr8.vcf.gz', 'chr9': 'new_output/snps/chr9.vcf.gz', 'chr10': 'new_output/snps/chr10.vcf.gz', 'chr11': 'new_output/snps/chr11.vcf.gz', 'chr12': 'new_output/snps/chr12.vcf.gz', 'chr13': 'new_output/snps/chr13.vcf.gz', 'chr14': 'new_output/snps/chr14.vcf.gz', 'chr15': 'new_output/snps/chr15.vcf.gz', 'chr16': 'new_output/snps/chr16.vcf.gz', 'chr17': 'new_output/snps/chr17.vcf.gz', 'chr18': 'new_output/snps/chr18.vcf.gz', 'chr19': 'new_output/snps/chr19.vcf.gz', 'chr20': 'new_output/snps/chr20.vcf.gz', 'chr21': 'new_output/snps/chr21.vcf.gz', 'chr22': 'new_output/snps/chr22.vcf.gz', 'chrX': 'new_output/snps/chrX.vcf.gz', 'chrY': 'new_output/snps/chrY.vcf.gz'}
    reference: data/hg19.fa
    j: 2
    q: 0
    Q: 11
    qual: 11
    E: False
    gamma: 0.05
    maxshift: 0.5
    mincov: 8
    maxcov: 300
    outputNormal: new_output/baf/normal.1bed
    outputTumors: new_output/baf/tumor.1bed
    outputSnps: new_output/count_alleles/snps
    verbose: False
[2022-Aug-12 19:47:28]# Counting SNPs alleles from the matched-normal sample
[2022-Aug-12 19:47:48]
[2022-Aug-12 19:47:48]
[2022-Aug-12 19:48:08]
[2022-Aug-12 19:48:08]
[2022-Aug-12 19:48:27]
[2022-Aug-12 19:48:27]
[2022-Aug-12 19:48:40]
[2022-Aug-12 19:48:40]
[2022-Aug-12 19:48:54]
[2022-Aug-12 19:48:54]
[2022-Aug-12 19:49:13]
[2022-Aug-12 19:49:13]
[2022-Aug-12 19:49:31]
[2022-Aug-12 19:49:31]
[2022-Aug-12 19:49:43]
[2022-Aug-12 19:49:43]
[2022-Aug-12 19:49:55]
[2022-Aug-12 19:49:55]
[2022-Aug-12 19:50:06]
[2022-Aug-12 19:50:06]
[2022-Aug-12 19:50:18]
[2022-Aug-12 19:50:30]
[2022-Aug-12 19:50:41]
[2022-Aug-12 19:50:53]
Progress: |████████████████████████████████████████| 100.0% Complete
Traceback (most recent call last):
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-demo/.snakemake/conda/62886ca61eec4c5786f8f5bdb83ae223/bin/hatchet", line 33, in <module>
    sys.exit(load_entry_point('hatchet==1.0.2', 'console_scripts', 'hatchet')())
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-demo/.snakemake/conda/62886ca61eec4c5786f8f5bdb83ae223/lib/python3.9/site-packages/hatchet/__main__.py", line 61, in main
    globals()[command](args)
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-demo/.snakemake/conda/62886ca61eec4c5786f8f5bdb83ae223/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 29, in main
    snps = counting(
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-demo/.snakemake/conda/62886ca61eec4c5786f8f5bdb83ae223/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 194, in counting
    results = {(v[0][0], v[0][1]): v for v in results}
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-demo/.snakemake/conda/62886ca61eec4c5786f8f5bdb83ae223/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 194, in <dictcomp>
    results = {(v[0][0], v[0][1]): v for v in results}
IndexError: list index out of range
[Fri Aug 12 19:50:57 2022]
Error in rule hatchet_count_alleles:
    jobid: 0
    output: new_output/baf/normal.1bed, new_output/baf/tumor.1bed, new_output/count_alleles/snps
    conda-env: /gpfs/gsfs9/users/Robinson-SB/HATCHet-demo/.snakemake/conda/62886ca61eec4c5786f8f5bdb83ae223
    shell:
        mkdir new_output/count_alleles/snps && hatchet count-alleles --tumors data/bulk_03clone1_06clone0_01normal.sorted.bam data/bulk_08clone1_Noneclone0_02normal.sorted.bam data/bulk_Noneclone1_09clone0_01normal.sorted.bam --normal data/normal.bam --samples normal tumor1 tumor2 tumor3 --reference data/hg19.fa --snps new_output/snps/*.vcf.gz --outputnormal new_output/baf/normal.1bed --outputtumors new_output/baf/tumor.1bed --outputsnps new_output/count_alleles/snps --mincov 8 --maxcov 300 
        (one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)

Removing output files of failed job hatchet_count_alleles since they might be corrupted:
new_output/count_alleles/snps
Shutting down, this might take some time.
Exiting because a job execution failed. Look above for error message

One extra bug is that the argument --outputsnps is described as the output file for the list of identified heterozygous germline SNPs here but from our previous discussion in this thread, it should actually be a directory?

I'd love to keep using HATCHet but it's unusable in its current form so please let me know when you've got a working demo or can help me understand why it isn't working on either my dataset or on the provided demo data.

Best, Welles

vineetbansal commented 2 years ago

Hi @alvinwt , @wir963 - sorry I was away for a couple of weeks. Thanks for all the comments and diagnostic information. I'll get together with the team this week and see how we can proceed to fix these issues. I think we should be able to fix these issue with count_alleles by this weekend. Stay tuned!

vineetbansal commented 2 years ago

Hello @wir963 - both your earlier message with the error:

    error('Allele counting failed on {} of {}, please check errors in {}!').format(
AttributeError: 'NoneType' object has no attribute 'format'

and the latter one with the error:

    results = {(v[0][0], v[0][1]): v for v in results}
IndexError: list index out of range

indicate that the count-alleles step failed for whatever reason (as you may have guessed already). The problem here is that there's no clear indication of what went wrong, and what one can do about it (due to another unrelated bug we found when debugging this). In HATCHet v1.0.3 released today (and available here and on bioconda), this error message would be replaced by what you should have seen all along - something like:

Worker 1 failed - please see Normal_chr22:30256931-32622323_bcftools.log for details.

I'm not sure if you are able to re-run your snakemake step that is failing, on short notice; if you can, I'd urge you to install HATCHet v1.0.3 and try again. Again, while this won't fix the error, this will give a better diagnostic message so you can proceed further.

However, you don't strictly need to do this if re-running the step would be time-consuming for you, and if you already see an error log file with the prefix _bcftools.log in your running directory. This file should provide additional clues as to what's going on.

Finally, if the *_bcftools.log file doesn't help and if you're willing, can you send me an email at vineetb [at] princeton.edu and we can schedule a zoom call to go through this? We'd love for you to continue using HATCHet, and once we understand the error better we can fix it relatively quickly

wir963 commented 2 years ago

Hey @vineetbansal,

I can re-run everything with a single script so it's a pretty low lift for me to re-run everything so happy to do whenever it would help. It took a little while because snakemake deletes the output directory if the run fails. I updated conda to use HATCHet v1.0.3 and here are the results.

For the homemade demo (https://github.com/wir963/HATCHet-demo), here's the updated error message (I'm not sure your fix worked for this error). 

Traceback (most recent call last):
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-demo/.snakemake/conda/87bb8bdff7b742fe7bfd8cce9b3dd3ae/bin/hatchet", line 33, in <module>
    sys.exit(load_entry_point('hatchet==1.0.3', 'console_scripts', 'hatchet')())
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-demo/.snakemake/conda/87bb8bdff7b742fe7bfd8cce9b3dd3ae/lib/python3.9/site-packages/hatchet/__main__.py", line 61, in main
    globals()[command](args)
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-demo/.snakemake/conda/87bb8bdff7b742fe7bfd8cce9b3dd3ae/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 29, in main
    snps = counting(
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-demo/.snakemake/conda/87bb8bdff7b742fe7bfd8cce9b3dd3ae/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 194, in counting
    results = {(v[0][0], v[0][1]): v for v in results}
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-demo/.snakemake/conda/87bb8bdff7b742fe7bfd8cce9b3dd3ae/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 194, in <dictcomp>
    results = {(v[0][0], v[0][1]): v for v in results}
IndexError: list index out of range
[Sat Aug 20 13:38:47 2022]
Error in rule hatchet_count_alleles:
    jobid: 0
    output: new_output/baf/normal.1bed, new_output/baf/tumor.1bed, new_output/count_alleles/snps
    conda-env: /gpfs/gsfs9/users/Robinson-SB/HATCHet-demo/.snakemake/conda/87bb8bdff7b742fe7bfd8cce9b3dd3ae
    shell:
        mkdir new_output/count_alleles/snps && hatchet count-alleles --tumors data/bulk_03clone1_06clone0_01normal.sorted.bam data/bulk_08clone1_Noneclone0_02normal.sorted.bam data/bulk_Noneclone1_09clone0_01normal.sorted.bam --normal data/normal.bam --samples normal tumor1 tumor2 tumor3 --reference data/hg19.fa --snps new_output/snps/*.vcf.gz --outputnormal new_output/baf/normal.1bed --outputtumors new_output/baf/tumor.1bed --outputsnps new_output/count_alleles/snps --mincov 8 --maxcov 300 
        (one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)

While trying to run our samples, here's the output

Traceback (most recent call last):
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/385bca84b190ae60a795776b13f0c436/bin/hatchet", line 33, in <module>
    sys.exit(load_entry_point('hatchet==1.0.3', 'console_scripts', 'hatchet')())
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/385bca84b190ae60a795776b13f0c436/lib/python3.9/site-packages/hatchet/__main__.py", line 61, in main
    globals()[command](args)
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/385bca84b190ae60a795776b13f0c436/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 29, in main
    snps = counting(
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/385bca84b190ae60a795776b13f0c436/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 193, in counting
    results = worker.run(work=work, n_instances=num_workers)
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/385bca84b190ae60a795776b13f0c436/lib/python3.9/site-packages/hatchet/utils/multiprocessing.py", line 97, in run
    raise handler_result.__class__(f'HANDLER {handler_id} FAILED\n\n{error_string}')
RuntimeError: HANDLER 0 FAILED

Traceback (most recent call last):
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/385bca84b190ae60a795776b13f0c436/lib/python3.9/site-packages/hatchet/utils/multiprocessing.py", line 31, in run
    result = self.worker.work(*args)
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/385bca84b190ae60a795776b13f0c436/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 225, in work
    snps = self.countAlleles(bamfile=bamfile, samplename=samplename, chromosome=chromosome)
  File "/gpfs/gsfs9/users/Robinson-SB/HATCHet-analysis/.snakemake/conda/385bca84b190ae60a795776b13f0c436/lib/python3.9/site-packages/hatchet/utils/count_alleles.py", line 262, in countAlleles
    raise ValueError(
ValueError: ESC[91mAllele counting failed on X of 4423N_PBL_October_02_2020, please check errors in output/HATCHet_v1/count-alleles/minreads_5_max_1000_min_20/germline_heterozygous_SNPs/4423N_PBL_October_02_2020_X_bcftools.log!ESC[0m

The output of 4423N_PBL_October_02_2020_X_bcftools.log is

[E::hts_open_format] Failed to open file "[]" : No such file or directory
bcftools mpileup: Could not read file "[]"Failed to read from standard input: unknown file type

I'm only passing VCF files for chromosomes 1-22 so I'm guessing using the sex chromosomes is an update for HATCHet v1? How does this work? Do you pass sex in as a parameter somewhere?

Best, Welles

wir963 commented 2 years ago

Hey @vineetbansal,

When I generate and pass VCF files for chromosomes X and Y (even for female patients), the count_alleles step works. I'm still getting a strange and unexplained error for count_alleles for the demo.

Should the outputSnps directory be empty after the run if there were no errors? Because mine is.

Best, Welles

vineetbansal commented 2 years ago

Hello @wir963, @alvinwt - thanks for your patience. We found the bug that was causing the IndexError in the count-alleles step and have fixed it in version HATCHet 1.1.0. This version should be available on bioconda soon, but I'm mentioning this here in case you want to check it out sooner from the master branch or releases.

This version also fixes the problems running the end-end demo of HATCHet. See the updated documentation at: https://compbio.cs.brown.edu/hatchet/script/README.html#command-for-running-the-hatchet-workflow

Essentally, setting the key [run].chromosomes (or the env var HATCHET_RUN_CHROMOSOMES) to chr22 should allow you to run the demo to completion. @wir963 - this should allow you to run your pipeline by specifying all but the sex chromosomes (space separated) as the value for this key, but I'd encourage you to try it out on a single chromosome first.

You were right that outputsnps specifies a folder, not a file - I've updated this in the documentation as well.

wir963 commented 2 years ago

Hey @vineetbansal ,

Thanks for getting back to me. I was able to run the updated demo using HATCHet 1.1.0. I'll start breaking down each step to make sure I understand the input and output and I'll try to run on my data and keep you posted.

One issue that you may want to fix. I got an error when I ran it initially - I had to change phase_snps = True to phase_snps = False, which wasn't mentioned in the demo. If we want to use phased snps, should we still use something like the Michigan Imputation server or is there a new update to HATCHet that does this?

Best, Welles

vineetbansal commented 2 years ago

Hi @wir963 - to avoid further frustrations with the count_alleles step, you will want to try out HATCHet 1.1.1, and skip 1.1.0. We're trying to get this released on bioconda for the past 2 days, but keep running into transient build errors (not related to HATCHet). If you have a way of cloning the repo/installing or installing from the release .zip, you may want to do that, or wait till tomorrow (hopefully) till HATCHET v1.1.1 is available.

Thanks for the other tip on phase_snps. I'll open an issue on it so we can keep track of it.

wir963 commented 2 years ago

Hey @vineetbansal,

I'm now running HATCHet 1.1.1. Two more things for you.

When I run count_alleles (https://github.com/wir963/HATCHet-demo/blob/main/Snakefile#L20) and it works, there are no files present in the --outputsnps directory. Is that okay? Shouldn't count_alleles output a file containing the location of heterozygous SNPs that is passed to count-reads using the --baffile argument?

In the count-alleles page, the argument --snps is listed once under Input and once under Optional parameters with very different descriptions. Which is correct?

count-reads lists --tumors but I get an error saying that I'm missing the argument --tumor.

Best, Welles

wir963 commented 1 year ago

Hey @vineetbansal ,

All of my remaining questions on this issue are more about documentation and clarity so I'll go ahead and close.

Best, Welles

vineetbansal commented 1 year ago

Hi @wir963 - let me look into these 2 concerns while we have your attention here. I hope to close this by the weekend. Thanks again for exercising HATCHet and letting us know of these issues!