fasta file in polarize command

[smallfishcui]

Hi Andrew,

I would like to add Ancestral allele info to my vcf file, so I would be able to get the derived allele for my site frequency spectrum. I was running a command:
python /homeappl/home/cuiwang/.local/bin/vcfdo polarize -i SNP_NOmissing_withOG.recode.vcf.gz -f SNP_NOmissing_withOG.recode.min4.fasta -n 10
INFO:vcfutils:Task: polarize
INFO:polarize:Connecting to VCF file <SNP_NOmissing_withOG.recode.vcf.gz>
INFO:polarize:Ancestral genome: <Y56.Nomissing.snps.fasta>
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=20190605
##source="Stacks v2.3e"
##INFO=<ID=AD,Number=R,Type=Integer,Description="Total Depth for Each Allele">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allele Depth">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
##FORMAT=<ID=GL,Number=G,Type=Float,Description="Genotype Likelihood">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##INFO=<ID=loc_strand,Number=1,Type=Character,Description="Genomic strand the corresponding Stacks locus aligns on">
##contig=<ID=36240>
##contig=<ID=36608>
##contig=<ID=36676>
##contig=<ID=38047>
##contig=<ID=42142>
##contig=<ID=55999>
##contig=<ID=57521>
##contig=<ID=59460>
##contig=<ID=68348>
##contig=<ID=77537>
##contig=<ID=78266>
##contig=<ID=79842>
##contig=<ID=88637>
##contig=<ID=91051>
##contig=<ID=101109>
##contig=<ID=109425>
##contig=<ID=109947>
##contig=<ID=133565>
##contig=<ID=145117>
##contig=<ID=150642>
##contig=<ID=152613>
##contig=<ID=176997>
##contig=<ID=180466>
##contig=<ID=180823>
##contig=<ID=187069>
##contig=<ID=192399>
##contig=<ID=194525>
##contig=<ID=205540>
##contig=<ID=207712>
##contig=<ID=221630>
##contig=<ID=226714>
##contig=<ID=231824>
##contig=<ID=245655>
##contig=<ID=251295>
##contig=<ID=279854>
##contig=<ID=280550>
##contig=<ID=290678>
##contig=<ID=299392>
##contig=<ID=313458>
##contig=<ID=331289>
##contig=<ID=342201>
##contig=<ID=353296>
##contig=<ID=358991>
##contig=<ID=370919>
##contig=<ID=379970>
##contig=<ID=387164>
##contig=<ID=401049>
##contig=<ID=413420>
##contig=<ID=419404>
##contig=<ID=427387>
##contig=<ID=430678>
##contig=<ID=430954>
##contig=<ID=431130>
##contig=<ID=443976>
##contig=<ID=446500>
##contig=<ID=452581>
##contig=<ID=456380>
##contig=<ID=459837>
##contig=<ID=466950>
##contig=<ID=2326531>
##vcfdo_command=polarize -i SNP_NOmissing_withOG.recode.vcf.gz -f Y56.Nomissing.snps.fasta -n 10; pwd=/wrk/cuiwang/DONOTREMOVE/Phrag_RAD/denovo/ustacks/M4m4/popu_K6; timestamp=2019-10-16 03:29:11
##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral allele inferred from outgroup">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Y1 Y3 Y4 Y6 Y7 Y9 Y10 Y11 Y12 Y13 Y14 Y15 Y16 Y17_1 Y18 Y19 Y20 Y21 Y22 Y23 Y24 Y25 Y26 Y27 Y28 Y30 Y31 Y32 Y33 Y34 Y35 Y36 Y37 Y38 Y39 Y40 Y41 Y43 Y44 Y45 Y46 Y47 Y48 Y49 Y51_1 Y53 Y55 Y56 Y57
However I get a error message:
INFO:polarize:Processing VCF file ...
Traceback (most recent call last):
File "/homeappl/home/cuiwang/.local/lib/python3.6/site-packages/pyfaidx/init.py", line 1014, in getitem
return self.records[rname]
KeyError: '36240'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/homeappl/home/cuiwang/.local/bin/vcfdo", line 2152, in
args.func(vcf = vcf, **vars(args))
File "/homeappl/home/cuiwang/.local/bin/vcfdo", line 810, in polarize_anc_vs_der
og = str(genome[site.CHROM][start:end])
File "/homeappl/home/cuiwang/.local/lib/python3.6/site-packages/pyfaidx/init.py", line 1016, in getitem
raise KeyError("{0} not in {1}.".format(rname, self.filename))
KeyError: '36240 not in Y56.Nomissing.snps.fasta.'

I set the fasta file to be the outgroup sequence, and 36240 was a name of a cotig, why the package is complaining the contig is not in the fasta file? Could you please give a hint how should I format the fasta file here?

Thanks,
Cui

[ChenqCQ]

Hi smallfishcui:
Sorry to bother you. I also want to add Ancestral allele info to my vcf file and get the derived allele site frequency spectrum. I meet a problem similar to yous when I run vcfdo, How did you solve it? My command is: vcfdo polarize -i /home/cg/data/vcf/SCGCg.vcf -f /home/cg/data/vcf/GC.fasta
The error is：‘Traceback (most recent call last):
File "/home/cg/miniconda3/envs/minicondainstall/bin/vcfdo", line 2219, in
args.func(vcf = vcf, **vars(args))
File "/home/cg/miniconda3/envs/minicondainstall/bin/vcfdo", line 804, in polarize_anc_vs_der
genome = Fasta(args.fasta)
File "/home/cg/miniconda3/envs/minicondainstall/lib/python3.8/site-packages/pyfaidx/init.py", line 998, in init
self.faidx = Faidx(
File "/home/cg/miniconda3/envs/minicondainstall/lib/python3.8/site-packages/pyfaidx/init.py", line 426, in init
self.read_fai()
File "/home/cg/miniconda3/envs/minicondainstall/lib/python3.8/site-packages/pyfaidx/init.py", line 493, in read_fai
raise ValueError('Duplicate key "%s"' % key)
ValueError: Duplicate key "SRR952792.1"’
Can you give me some suggestions? Thank you so much.

[andrewparkermorgan]

It looks like you have duplicate sequence names in the fasta file you are using as the ancestral sequence.

…

On Oct 25, 2021, at 4:24 AM, chenqing-1996 ***@***.***> wrote:  Hi smallfishcui: Sorry to bother you. I also want to add Ancestral allele info to my vcf file and get the derived allele site frequency spectrum. I meet a problem similar to yous when I run vcfdo, How did you solve it? My command is: vcfdo polarize -i /home/cg/data/vcf/SCGCg.vcf -f /home/cg/data/vcf/GC.fasta The error is：‘Traceback (most recent call last): File "/home/cg/miniconda3/envs/minicondainstall/bin/vcfdo", line 2219, in args.func(vcf = vcf, **vars(args)) File "/home/cg/miniconda3/envs/minicondainstall/bin/vcfdo", line 804, in polarize_anc_vs_der genome = Fasta(args.fasta) File "/home/cg/miniconda3/envs/minicondainstall/lib/python3.8/site-packages/pyfaidx/init.py", line 998, in init self.faidx = Faidx( File "/home/cg/miniconda3/envs/minicondainstall/lib/python3.8/site-packages/pyfaidx/init.py", line 426, in init self.read_fai() File "/home/cg/miniconda3/envs/minicondainstall/lib/python3.8/site-packages/pyfaidx/init.py", line 493, in read_fai raise ValueError('Duplicate key "%s"' % key) ValueError: Duplicate key "SRR952792.1"’ Can you give me some suggestions? Thank you so much. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android.

[ChenqCQ]

Hi andrew,
Thank you for your reply. My fasta file is converted by paired-sequencing fastq file. Do I need to remove duplicate names in this case?

[smallfishcui]

Hi, Sorry, it has been a long time that I tried this that I couldn't remember what happened later on. But I m afraid I wont be of any help here because I didn't go through with vcfdo eventually.... best, smallfishcui chenqing-1996 ***@***.***> 于2021年10月25日周一 下午2:42写道：

…

Hi andrew, Thank you for your reply. My fasta file is converted by paired-sequencing fastq file. Do I need to remove duplicate names in this case? — You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub <#6 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AJZNAZJQKNRZ2QIGRQBBX23UIU7DPANCNFSM4JBEV6CA> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[andrewparkermorgan]

vcfdo expects the ancestral and references sequences to be 'parallel' (ie to share the same coordinate system), and requires that sequence names match. (The assumption that the ancestral and reference genomes differ only by single-nucleotide changes is obviously not realistic, but relaxing that assumption is outside the intended scope of vcfdo.) Consider an example of this VCF entry:

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT ancestor
chr1 4 . A G . PASS . GT 1/1 

The fasta files for reference and ancestral would be expected to look like this:

reference.fa:
>chr1
...A...

ancestral.fa:
>chr1
...G...

In your case, it looks like you have used a sample called SRR952792.1 to generate an ancestral genome file, but you will need to edit the resulting file so that the sequence names (that is, the chromosome or contig names) match those in the reference genome.

[ChenqCQ]

Hi andrew,
Thank you very much. I think I know how to do that. Thank you again.

[ChenqCQ]

Hi andrew,
According to your suggestion yesterday, I converted the VCF file of my outgroup into fasta file to make the outgroup consistent with the reference genome. This is the command I used:

bgzip -c my.vcf > my.vcf.gz
tabix -p vcf my.vcf.gz
cat ref.fasta | bcftools consensus my.vcf.gz > my.fa

Then, I polarized my vcf file, but the final results I got showed that:

0 sites have undefined ancestral allele
0 sites have REF == ancestral
2245523 sites have REF == derived (100.000%)

Ancestral loci should not be zero, which is strange,.What do you think?

[andrewparkermorgan]

Can you show the command you used to polarize the VCF, and maybe also paste in the header of the VCF? If your VCF file has only one sample (representing an outgroup) with variants called against the reference genome and you use that to define ancestral alleles, then run the same file through vcfdo polarize, then the ancestral alleles will all be ALT by definition, which means derived alleles will be REF.

[ChenqCQ]

I check my vcf file, it is ok. its header contains samples : S1, S2, S3, S4, S5, S6, S7,S8,S9, S10, S11, S12. and my command is:
vcfdo polarize -i /home/cg/data/vcf/SW.vcf -f /home/chenqing/data/vcf/GC.fasta
Is taht a problem?

[andrewparkermorgan]

I think there is still some misunderstanding about what your "ancestral" fasta represents, and that vcfdo polarize is working as I intend but maybe not giving you what you need. I am going to try to lay out as explicitly as I can the scenario(s) I had in mind when I designed the tool. You probably understand this already so I apologize for being pedantic but I want to put it out there so we have a shared understanding of the problem.

In general it is not easy to get the true ancestral genome sequence. There are (at least) two options: one is a formal ancestral state reconstruction from a whole-genome multiple sequence alignment. This is ideal but time-consuming and assumes the presence of such an alignment. Depending on the species of interest, such data may not be readily available. Second option is to use WGS data from an outgroup, aligned to the reference genome for the ingroup, and use that to assign likely ancestral states at sites that are polymorphic in the ingroup. The underlying assumption is this: if the TMRCA between outgroup and ingroup is relatively short, the probability that the same site has mutated twice independently since the MRCA is very small. So if a site is polymorphic in the ingroup, it is unlikely also polymorphic in the outgroup, and genotype any random individual from the outgroup should represent the ancestral state. To use this method, the conditioning on being polymorphic in the ingroup is key.

To see what I mean, consider these three example trees. Let O be the outgroup and {A1, A2, A3} be individuals from the ingroup, and let * represent the mutation of interest. In scenario (1), the site is polymorphic in the ingroup and the genotype of O gives the ancestral state.

(1)
   |
.-----.
|     |
|  .-----.
|  |     |
|  |     *
|  |     |
|  |   .---.
|  |   |   |
O  A1  A2  A3

In scenarios (2) and (3), all ingroup individuals have the same allele. In scenario (2) the ingroup is fixed for the derived allele whereas in scenario (3) the ingroup is fixed for the derived allele. When aligning to a reference genome, the ingroup genotypes could be all homozygous ALT or REF depending on whether the reference genome represents ancestral or derived state. Without a third species to use as a tie-breaker, the ancestral state is ambiguous. Evening having multiple outgroup individuals doesn't necessarily break the tie.

(2)
   |
.-----.
|     |
|     *
|     |
|  .-----.
|  |     |
|  |   .---.
|  |   |   |
O  A1  A2  A3

(3)
   |
.-----.
|     |
|  .-----.
|  |     |
|  |     |
*  |     |
|  |   .---.
|  |   |   |
O  A1  A2  A3

Taking all this into account, one expected use of vcfdo polarize is as follows. First, ascertain variant sites using WGS data from the ingroup. Then genotype these sites in WGS data from the outgroup. Use those calls to polarize alleles in the ingroup.

Getting back to your own data, what is the overlap between the samples in your outgroup VCF and the VCF you are trying to polarize? I think that is where the confusion arises. And maybe from the behavior of bcftools consensus, which is to just insert all ALT alleles into the consensus fasta, which may not be what you want.

[Neato-Nick]

@andrewparkermorgan I've bookmarked your comment and come back to it multiple times. I used angsd, which uses parsimony among multiple outgroups, to infer an ancestral genome fasta with SNP data. Then vcfdo polarize on that fasta.
I'm commenting because I have one correction to make for others who find this issue helpful:

In scenarios (2) and (3), all ingroup individuals have the same allele. In scenario (2) the ingroup is fixed for the derived allele whereas in scenario (3) the ingroup is fixed for the derived allele.

Scenario 3 should read either "ingroup is fixed for ancestral allele" or "outgroup is fixed for derived allele," depending on where you want the emphasis. Thank you for this fantastic explanation of pop gen theory. IMO you should consider adding this to the README somewhere.