extract_phase_bcftools.sh

[Ptero64]

Hello,
I am trying to use your script to phase targeted genes form some samples. I followed your tutorial to generate the different files with GATK then whatshap. However I am not sure to really understand your script using bcftools consensus. (For now I am running the different .sh .py for 1 sample).
The script extract_phase_bcftools.sh is suppose to generate 1 fasta file per gene including 2 sequences (1 per allele; e.g. Sample-gene01.phased.fasta)?
But to run haplonerate.py I need the gtf file form whatshap and 2 fasta files containing sequences for 1 or more genes. So these input fasta are not the same that the output of extract_phase_bcftools.sh?

Thanks for the help.

Nicolas

[Ptero64]

Hi,
Also notice that you must be careful using the script extract_phase_bcftools.sh. Starting line 8:
mkdir -p $prefix
cd $prefix
rm -r * #! don't know why there is this line? and if there is an issue this you $prefix (typo or not founded) it's turn very bad :(

[mossmatters]

Hello,

The extract_phase_bcftools.sh script will generate two FASTA files. Each FASTA file will have one "allele" for each gene, so if you have 100 genes, there will be 100 sequences in each file. The script is following this protocol on the WhatsHap user guide

You then give these two files and the GTF from Whatshap to haplonerate.py.

That is a good point about the rm -r *! I recently learned that if you add:

set -eo pipefail

at the beginning of the bash script, it will exit whenever any commands return an error. I will update all the bash scripts to have this.

[Ptero64]

Hello,

Thank you very much for your reply and explanations.

Ok that's what I thought for the output of bcftools.sh. But currently the code save the two "allele" in a same file (and the 2 seqs have the same name >$prefix-gene).
line 19:
bcftools consensus -H 2 $prefix.supercontigs.fasta.snps.whatshap.vcf.gz >> phased_bcftools/$prefix-{1}.phased.fasta" :::: $genelist

Is it rather > instead of >>? (and changing output name?).

Thanks

[rubimeza]

Hi Nicolas and Matt
I did exactly what Nicolas suggested:

Is it rather > instead of >>? (and changing output name?)

And worked well, I think.

Now I'm trying to use haplonerate.py, but nothing happens.
Apparently haplonerate recognizes the files and commands, but it produces only an empty outfile.
I'm out of ideas here, so I'll attach the files I'm using.
And really hope you can help me out.
Thanks in advance
Rubi

para probar script.zip