This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Linkage disequilibrium from VCF file of variants in bacterial metagenomes

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
snps metagenomes ld bacteria vcf

0 answers

No answers yet.

Log in to answer this question.