Issue Running ProHap on All-of-Us Dataset

[mengysun]

Hi @vasicek58 ,

I am the collaborator of @bschilder, who opened a couple of issues before. I'm currently running ProHap on the All-of-Us dataset (414,000 individuals) and encountered an issue that's been a bit puzzling. For context, the pipeline exits with status code 1 and provides minimal error information.

Here’s the relevant error message from the log (from /home/jupyter/workspaces/gpmaps2/ProHap/Snakefile_chr17, line 277):

CalledProcessError:
Command 'source /home/jupyter/miniconda3/bin/activate '/home/jupyter/workspaces/gpmaps2/ProHap/.snakemake/conda/8176726ab2db4e11e088d31f55112ddf_'; set -euo pipefail; mkdir -p tmp/transcript_vcf_haplo; mkdir -p log; mkdir -p results; python3 src/prohap.py -i "data/vcf/phased_0000014706/chr17_phased_filtered.vcf.gz" -db "data/gtf/Homo_sapiens.GRCh38.113.chr_patch_hapl_scaff_chr17.db" -transcripts "data/transcripts_selected.csv" -cdna "data/fasta/total_cdnas_113.fa" -s "data/AoU_meta.tsv" -chr 17 -min_hap_freq 0 -min_hap_count 0 -acc_prefix enshap_17 -id_prefix haplo_chr17 -require_start 1 -ignore_UTR 1 -skip_start_lost 1 -x_par1_to 2781479 -x_par2_from 155701383 -threads 50 -log "log/prohap_chr17.log" -tmp_dir "tmp/transcript_vcf_haplo" -output_csv "results/All_of_Us_haplotypes_0000014706_tmp/haplo_chr17.tsv.gz" -output_fasta "results/All_of_Us_haplotypes_0000014706_tmp/haplo_chr17.fa" -output_cdna_fasta "results/All_of_Us_haplotypes_0000014706_tmp/haplo_cdna_chr17.fa"'
returned non-zero exit status 1.

Interestingly, a .tsv file is successfully generated in the tmp/transcript_vcf_haplo directory before the error occurs, so the failure seems to happen right after that step.

I’m wondering if you might have any insight into potential causes? One thought I had was whether the issue could stem from a mismatch between the expected file extension or format (e.g., .csv vs .tsv). However, I didn’t see anywhere in the Snakemake file where the output format is explicitly specified or could be misconfigured.

Would appreciate any guidance or suggestions you might have!

Best regards

[vasicek58]

Hi @mengysun,

Thank you for raising the two issues, and so sorry for my delayed response! Things are quite busy on my end until the end of the month - would it be okay if I look closely into this in early June?

All the best,
Jakub

[mengysun]

Sure. Thanks!
Best,
Mengyi

[mengysun]

Hi Jakub,

Just following up — I’ve realized that part of the issue stems from needing to specify the gender in the metadata, not just the population code. It’s been a bit tricky to find this information for the "All of Us" sample.

Would it be possible to run the code without specifying gender? Alternatively, would it cause any problems if I temporarily assign a placeholder value like "All" or even randomly specify the gender, just to test the run?

Best,
Mengyi

[mengysun]

Ok I sort of found the file of gender information. But probably if there's an option that can run without that would be good for fast testing.

[vasicek58]

Hi! For fast testing, you could just assign everyone as male. It would only affect the X chromosome, where the second haplotype will be ignored for female genotypes. Do you have Y chromosomes too in the dataset?

[mengysun]

Thanks! I think I fix that issue by identify the files related to the gender information of the samples. Now I run into a new issue:

"File "/home/jupyter/workspaces/gpmaps2/ProHap/src/prohap.py", line 146, in
gene_haplo_df = get_gene_haplotypes(all_transcripts, sample_ids, args.tmp_dir, args.log_file, args.threads, (args.chromosome == 'X'), args.x_par1_to, args.x_par2_from, samples_df)
File "/home/jupyter/workspaces/gpmaps2/ProHap/src/modules/get_haplotypes.py", line 316, in get_gene_haplotypes
result_df['Frequency_population'] = result_df.apply(count_freq_pop, axis=1)
File "/home/jupyter/.local/lib/python3.10/site-packages/pandas/core/frame.py", line 4301, in setitem
self._set_item_frame_value(key, value)
File "/home/jupyter/.local/lib/python3.10/site-packages/pandas/core/frame.py", line 4459, in _set_item_frame_value
raise ValueError(
ValueError: Cannot set a DataFrame with multiple columns to the single column Frequency_population"

I got a bit confused by this issue, but I wonder if I skip the population analysis (for example, by setting the population code and superpopulation code to 'ALL'), whether or not I can bypass it.

[mengysun]

I believe the error mentioned above might be resolved by modifying the get_haplotypes.py script. Specifically, update the following lines:

result_df['Frequency'] = result_df.apply(count_freq, axis=1)
result_df['Frequency_population'] = result_df.apply(count_freq_pop, axis=1)
result_df['Frequency_superpopulation'] = result_df.apply(count_freq_superpop, axis=1)

to:

result_df['Frequency'] = result_df.apply(count_freq, axis=1)

freq_df = result_df.apply(count_freq_pop, axis=1, result_type="expand")
freq_df = freq_df.add_prefix("Frequency_")
result_df = pd.concat([result_df, freq_df], axis=1)

freq_df = result_df.apply(count_freq_superpop, axis=1, result_type="expand")
freq_df = freq_df.add_prefix("Frequency_")
result_df = pd.concat([result_df, freq_df], axis=1)

At least in my case, this change appears to work around the error—I was able to run the entire pipeline to completion with outputs. However, I tested the above under the specification where both the population code and super-population code are set to 'all'. I’ll test again with standard population code settings to confirm if it still works. That said, this might already provide some clues as to what’s going wrong.
@vasicek58

[vasicek58]

Hi @mengysun, apologies for the delay once again, and thank you for looking into this and suggesting a solution!

This error is still puzzling me a little bit, the count_freq_pop function should be returning a single string value per row, and that shouldn't be interpreted as a multi-column DataFrame by pandas. Did you by any chance make any changes to the count_freq_pop and / or count_freq_superpop functions?