I computed linkage disequilibrium (r2) in a bacterial species using the --hap-r2 method in vcftools. The manual says it assumes that the VCF input file has phased haplotypes, which makes me wonder if the method is applicable for the haploid bacteria. However, the results I got seem to makes sense: r2 decreases with distance between the positions and approaches a baseline at around 1000 bp distance, which is in the range of what I'd expect from other analyses inferring the sequence length of homologous recombination events.
Now my question is if anyone knows roughly how vcftools computes r2 and if it's reasonable to use it for haploid organisms, as I couldn't find this information in the documentation.
My VCF file was obtained from mapping the reads of metagenomes from 40 different lake water samples against a reference genome and calling the SNPs. So the VCF file contains 40 samples, but I saw the hap-r2 output calls them 40 individuals. Yet, each sample represents a population and not one individual haplotype. For instance, at position 611 there are 0 A, 6 T, 9 C, 0 G in sample1, 0 A, 0 T, 15 C, 0 G in sample2, 0 A, 10 T, 5 C, 0 G in sample3 etc. If each sample was an individual (haplotype) I suppose it would only have one allele at each position (e.g. 0 A, 0 T, 1 C, 0 G), right? If that's the case my VCF file doesn't conform to the usual input for r2-hap...
Does anyone know how vcftools might have computed r2 from my VCF file then?
Below the first lines of my VCF file, the log file and first lines of the results I get from:
vcftools --vcf myfile.vcf --hap-r2 --ld-window-bp 10000 --out ld_window_10000
First lines of myfile.vcf (excluding ##INFO and ##FORMAT lines):
##fileformat=VCFv4.2
##fileDate=20220728
##source=freeBayes v1.3.1-dirty
##reference=RAW_DATA/Genomes/Ppan/UB-Kaiv-W7.fa
##contig=<ID=UB-Kaiv-W7,length=1830921>
##phasing=none
##commandline="freebayes -f RAW_DATA/Genomes/Ppan/UB-Kaiv-W7.fa 05_BAM_merged/Ppan/params_cov_30_bdth_50_subsamp_TRUE_mpq_20_bq_15/UB-Kaiv-W7_merged_sorted_position.bam -C 4 -p 1 --pooled-continuous --read-max-mismatch-fraction 0.05 --min-alternate-fraction 0.01 -q 15 --report-monomorphic"
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT EMTI_200616 AlpFen_Clou_10July2018 AM-lipid-02-D3 EMTI_200818 EMOF_200619 EMOF_200714 AM-lipid-01-D2 EMTI_200714 Umea2p1 EMOF_200617 EMOF_200615 EMTI_200615 AlpFen_Aci1_07Aug2018 EMTI_200618 AlpFen_Clou_11June2018 SRR3657436 EMOF_200616 AM-lipid-02-D2 AlpFen_EnzMain_19July2018 EMOF_180719 EMOF_180528 SRR3657429 AlpFen_Gerlos4_19July2018 EMTI_201023 AlpFen_Hutt_28May2018 EMTI_200619 EMTI_200617 EMOF_200913 EMOF_200618 AlpFen_Troeg7_08July2018 EMTI_200913 AlpFen_Latsch_28May2018 AM-lipid-02-D1 AlpFen_Roav4_06Aug2018 MJ130820B AlpFen_Pas2_05Aug2018 AlpFen_Gerlos5_19July2018 EMOF_200818 AlpFen_EnzMain_28May2018 Ki1-0-7m
UB-Kaiv-W7 162 . G A 4843.8 . AB=0;ABP=0;AC=19;AF=0.475;AN=40;AO=520;CIGAR=1X;DP=1045;DPB=1045;DPRA=2.23904;EPP=46.4565;EPPR=102.761;GTI=0;LEN=1;MEANALT=1.02632;MQM=30.6731;MQMR=34.0497;NS=40;NUMALT=1;ODDS=0.650087;PAIRED=1;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;QA=18895;QR=19081;RO=523;RPL=263;RPP=3.16063;RPPR=8.09644;RPR=257;RUN=1;SAF=372;SAP=212.541;SAR=148;SRF=357;SRP=154.478;SRR=166;TYPE=snp;technology.illumina=1 GT:DP:AD:RO:QR:AO:QA:GL 1:44:17,27:17:629:27:953:-17.1181,0 1:25:11,14:11:407:14:518:-18.5895,0 1:28:9,19:9:321:19:679:-20.4704,0 0:20:12,8:12:432:8:284:0,-4.1718 1:35:15,20:15:555:20:740:-8.98318,0 1:20:1,19:1:37:19:691:-44.0818,0 1:8:0,8:0:0:8:296:-24.6777,0 1:30:1,29:1:37:29:1049:-70.9385,0 1:46:19,27:19:703:27:987:-20.4106,0 0:42:22,20:22:802:20:740:0,-10.693 0:41:26,15:26:950:15:543:0,-40.4758 1:28:10,18:10:370:18:642:-6.14721,0 0:15:8,7:8:296:7:259:0,-6.5424 1:17:6,11:6:222:11:392:-7.92128,0 0:26:22,4:22:802:4:148:0,-45.2113 0:33:22,11:22:810:11:379:0,-42.2545 0:32:21,11:21:765:11:407:0,-28.051 1:21:8,13:8:296:13:469:-11.5018,0 1:17:4,13:4:148:13:481:-22.0559,0 1:17:4,13:4:148:13:481:-22.0559,0 1:32:1,31:1:37:31:1123:-77.1867,0 0:32:26,6:26:912:6:218:0,-63.4818 0:17:16,1:16:580:1:37:0,-38.8954 0:27:25,2:25:913:2:62:0,-64.5382 1:22:11,11:11:383:11:407:-7.97025,0 1:18:6,12:6:222:12:444:-13.5026,0 0:30:16,14:16:580:14:518:0,-13.0582 0:20:14,6:14:518:6:222:0,-18.7604 0:49:22,25:22:814:25:913:0,-9.58021 0:15:15,0:15:555:0:0:0,-39.339 0:27:23,4:23:839:4:148:0,-48.9739 1:15:4,11:4:148:11:383:-15.8752,0 1:21:4,17:4:148:17:605:-34.6146,0 0:9:9,0:9:321:0:0:0,-22.5805 0:35:16,19:16:580:19:703:0,-0.35746 0:15:12,3:12:432:3:111:0,-20.9927 0:17:16,1:16:592:1:37:0,-34.5289 0:37:22,15:22:802:15:555:0,-22.3924 1:32:1,31:1:37:31:1123:-77.1867,0 0:30:26,4:26:938:4:148:0,-60.3239
UB-Kaiv-W7 174 . T C 5283.81 . AB=0;ABP=0;AC=22;AF=0.55;AN=40;AO=516;CIGAR=1X;DP=994;DPB=994;DPRA=1.6735;EPP=40.1947;EPPR=122.232;GTI=0;LEN=1;MEANALT=1;MQM=30.3973;MQMR=33.91;NS=40;NUMALT=1;ODDS=5.39987;PAIRED=1;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;QA=18948;QR=17028;RO=478;RPL=247;RPP=5.04711;RPPR=12.6229;RPR=269;RUN=1;SAF=356;SAP=164.676;SAR=160;SRF=298;SRP=66.2647;SRR=180;TYPE=snp;technology.illumina=1 GT:DP:AD:RO:QR:AO:QA:GL 1:38:16,22:16:580:22:802:-6.56128,0 1:21:8,13:8:284:13:481:-21.9281,0 1:32:10,22:10:370:22:802:-25.7499,0 1:19:10,9:10:346:9:333:-3.18881,0 1:29:11,18:11:383:18:666:-17.0081,0 1:20:0,20:0:0:20:740:-49.4428,0 1:8:0,8:0:0:8:296:-24.6777,0 1:25:1,24:1:37:24:888:-57.9259,0 1:41:16,25:16:580:25:925:-24.4418,0 1:44:20,24:20:680:24:882:-4.11793,0 0:37:25,12:25:913:12:444:0,-47.4209 1:28:9,19:9:321:19:691:-12.4895,0 0:16:10,6:10:358:6:222:0,-11.69 1:18:5,13:5:173:13:481:-15.4243,0 0:22:18,4:18:606:4:148:0,-32.7204 0:34:25,9:25:920:9:327:0,-56.5553 0:30:19,11:19:667:11:407:0,-23.1971 1:21:7,14:7:259:14:518:-17.8002,0 1:13:3,10:3:111:10:370:-16.7831,0 1:13:3,10:3:111:10:370:-16.7831,0 1:32:1,31:1:37:31:1123:-78.2206,0 0:35:29,6:29:1056:6:234:0,-72.3261 0:16:14,2:14:494:2:74:0,-30.6728 0:23:21,2:21:753:2:74:0,-53.8767 1:21:9,12:9:333:12:444:-13.225,0 1:18:7,11:7:235:11:395:-10.4866,0 0:32:15,17:15:531:17:629:0,-1.70167 0:19:12,7:12:444:7:259:0,-8.08624 0:48:23,25:23:791:25:913:0,-8.9099 0:15:15,0:15:531:0:0:0,-35.3309 0:19:15,4:15:543:4:148:0,-24.4005 1:16:3,13:3:111:13:469:-26.3921,0 1:22:4,18:4:148:18:654:-37.0133,0 0:9:8,1:8:272:1:37:0,-15.957 1:36:16,20:16:568:20:740:-1.95719,0 0:11:8,3:8:248:3:111:0,-8.72046 0:16:15,1:15:531:1:37:0,-32.1763 0:35:21,14:21:777:14:506:0,-21.5744 1:32:1,31:1:37:31:1123:-78.2206,0 0:30:25,5:25:889:5:185:0,-51.9286
UB-Kaiv-W7 224 . GCGC ACGT,ACGC 1094.34 . AB=0,0;ABP=0,0;AC=0,10;AF=0,0.25;AN=40;AO=10,283;CIGAR=1X2M1X,1X3M;DP=954;DPB=964.5;DPRA=0.988038,0.841856;EPP=3.0103,66.5508;EPPR=7.77566;GTI=0;LEN=4,1;MEANALT=2,1.24242;MQM=24.8,29.5371;MQMR=32.6337;NS=40;NUMALT=2;ODDS=2.71402;PAIRED=1,1;PAIREDR=1;PAO=9,7;PQA=331,257;PQR=438;PRO=12;QA=368,7928;QR=23962;RO=658;RPL=4,192;RPP=3.87889,81.2831;RPPR=14.1119;RPR=6,91;RUN=1,1;SAF=5,136;SAP=3.0103,3.93874;SAR=5,147;SRF=396;SRP=62.2671;SRR=262;TYPE=complex,snp;technology.illumina=1,1 GT:DP:AD:RO:QR:AO:QA:GL 0:38:24,0,13:24:879:0,13:0,346:0,-58.478,-30.4761 0:24:13,0,11:13:467:0,11:0,272:0,-40.6429,-19.079 0:40:40,0,0:40:1456:0,0:0,0:0,-112.342,-112.342 0:24:14,0,10:14:506:0,10:0,246:0,-36.6228,-16.8426 0:27:19,0,8:19:700:0,8:0,214:0,-43.9381,-26.6227 0:28:28,0,0:28:1030:0,0:0,0:0,-65.7736,-65.7736 0:14:14,0,0:14:515:0,0:0,0:0,-41.3471,-41.3471 0:38:37,0,1:37:1337:0,1:0,23:0,-97.7006,-95.903 0:33:18,0,15:18:657:0,15:0,457:0,-48.4436,-10.4808 0:33:21,0,12:21:768:0,12:0,306:0,-52.4849,-27.9823 0:29:17,0,12:17:626:0,12:0,323:0,-46.2814,-19.8353 0:25:19,0,6:19:697:0,6:0,208:0,-42.9153,-25.1377 0:12:7,1,4:7:249:1,4:37,135:0,-17.2967,-7.93454 0:17:12,0,5:12:438:0,5:0,157:0,-29.4101,-16.2907 2:21:4,0,17:4:131:0,17:0,439:-25.1111,-35.2009,0 0:27:21,0,6:21:791:0,6:0,200:0,-59.1629,-42.5118 2:27:12,0,15:12:444:0,15:0,418:-6.2714,-34.6423,0 0:26:26,0,0:26:953:0,0:0,0:0,-77.4054,-77.4054 0:14:13,0,1:13:468:0,1:0,23:0,-32.3154,-30.5207 0:14:13,0,1:13:468:0,1:0,23:0,-32.3154,-30.5207 0:29:29,0,0:29:1035:0,0:0,0:0,-76.0292,-76.0292 0:32:29,1,2:29:1080:1,2:38,63:0,-83.9901,-80.8386 2:14:4,0,10:4:148:0,10:0,276:-9.675,-22.4908,0 2:20:0,0,20:0:0:0,20:0,615:-51.2788,-51.2788,0 0:14:7,0,6:7:246:0,6:0,182:0,-19.5651,-4.48008 0:16:12,0,4:12:441:0,4:0,94:0,-28.1946,-20.8287 0:33:26,0,7:26:947:0,7:0,177:0,-65.3137,-51.1235 2:19:8,0,11:8:293:0,11:0,310:-3.48755,-25.5017,0 0:37:24,0,13:24:862:0,13:0,443:0,-56.6479,-18.3713 2:16:4,0,12:4:145:0,12:0,296:-12.2196,-23.6953,0 2:16:5,0,10:5:182:0,10:0,261:-7.35545,-20.9727,0 0:14:12,0,2:12:431:0,2:0,60:0,-34.0302,-29.0644 0:26:26,0,0:26:956:0,0:0,0:0,-77.2535,-77.2535 2:7:1,1,5:1:37:1,5:37,144:-8.68771,-9.87173,0 0:36:27,4,5:27:980:4,5:145,159:0,-64.5264,-59.5158 2:10:3,0,7:3:111:0,7:0,175:-4.22069,-13.8775,0 2:12:3,0,9:3:111:0,9:0,236:-10.522,-19.0251,0 0:32:21,0,11:21:750:0,11:0,340:0,-53.4893,-25.0882 0:29:29,0,0:29:1035:0,0:0,0:0,-76.0292,-76.0292 0:31:16,3,12:16:592:3,12:111,307:0,-41.7276,-24.433
ld_window_10000.log:
VCFtools - 0.1.16
(C) Adam Auton and Anthony Marcketta 2009
Parameters as interpreted:
--vcf myfile.vcf
--ld-window-bp 10000
--max-alleles 2
--min-alleles 2
--hap-r2
--out ld_window_10000
--phased
Warning: Expected at least 2 parts in INFO entry: ID=AF,Number=A,Type=Float,Description="Estimated allele frequency in the range (0,1]">
Warning: Expected at least 2 parts in INFO entry: ID=PRO,Number=1,Type=Float,Description="Reference allele observation count, with partial observations recorded fractionally">
Warning: Expected at least 2 parts in INFO entry: ID=PAO,Number=A,Type=Float,Description="Alternate allele observations, with partial observations recorded fractionally">
Warning: Expected at least 2 parts in INFO entry: ID=SRP,Number=1,Type=Float,Description="Strand balance probability for the reference allele: Phred-scaled upper-bounds estimate of the probability of observing the deviation between SRF and SRR given E(SRF/SRR) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=SAP,Number=A,Type=Float,Description="Strand balance probability for the alternate allele: Phred-scaled upper-bounds estimate of the probability of observing the deviation between SAF and SAR given E(SAF/SAR) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=AB,Number=A,Type=Float,Description="Allele balance at heterozygous sites: a number between 0 and 1 representing the ratio of reads showing the reference allele to all reads, considering only reads from individuals called as heterozygous">
Warning: Expected at least 2 parts in INFO entry: ID=ABP,Number=A,Type=Float,Description="Allele balance probability at heterozygous sites: Phred-scaled upper-bounds estimate of the probability of observing the deviation between ABR and ABA given E(ABR/ABA) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=RPP,Number=A,Type=Float,Description="Read Placement Probability: Phred-scaled upper-bounds estimate of the probability of observing the deviation between RPL and RPR given E(RPL/RPR) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=RPPR,Number=1,Type=Float,Description="Read Placement Probability for reference observations: Phred-scaled upper-bounds estimate of the probability of observing the deviation between RPL and RPR given E(RPL/RPR) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=EPP,Number=A,Type=Float,Description="End Placement Probability: Phred-scaled upper-bounds estimate of the probability of observing the deviation between EL and ER given E(EL/ER) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=EPPR,Number=1,Type=Float,Description="End Placement Probability for reference observations: Phred-scaled upper-bounds estimate of the probability of observing the deviation between EL and ER given E(EL/ER) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
Warning: Expected at least 2 parts in INFO entry: ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
Warning: Expected at least 2 parts in INFO entry: ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
Warning: Expected at least 2 parts in INFO entry: ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
Warning: Expected at least 2 parts in INFO entry: ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
Warning: Expected at least 2 parts in INFO entry: ID=CIGAR,Number=A,Type=String,Description="The extended CIGAR representation of each alternate allele, with the exception that '=' is replaced by 'M' to ease VCF parsing. Note that INDEL alleles do not have the first matched base (which is provided by default, per the spec) referred to by the CIGAR.">
Warning: Expected at least 2 parts in FORMAT entry: ID=GQ,Number=1,Type=Float,Description="Genotype Quality, the Phred-scaled marginal (or unconditional) probability of the called genotype">
Warning: Expected at least 2 parts in FORMAT entry: ID=GL,Number=G,Type=Float,Description="Genotype Likelihood, log10-scaled likelihoods of the data given the called genotype for each possible genotype generated from the reference and alternate alleles given the sample ploidy">
After filtering, kept 40 out of 40 Individuals
Outputting Pairwise LD (phased bi-allelic only)
After filtering, kept 22848 out of a possible 25742 Sites
Run Time = 33.00 seconds
First lines of ld_window_10000.hap.ld:
CHR POS1 POS2 N_CHR R^2 D Dprime
UB-Kaiv-W7 162 174 40 0.74026 0.21375 1
UB-Kaiv-W7 162 253 40 0.047619 -0.02375 -1
UB-Kaiv-W7 162 269 40 0.159664 -0.07125 -1
UB-Kaiv-W7 162 288 40 0.129252 -0.059375 -1
0 answers
No answers yet.
Log in to answer this question.