Page 3 of 5 FirstFirst 12345 LastLast
Results 21 to 30 of 49

Thread: Updated BAM Analysis Kit. Any interest?

  1. #21
    Moderator
    Posts
    5,035
    Sex
    Location
    Normandy
    Ethnicity
    northwesterner
    Y-DNA
    U152>L2>Z367
    mtDNA
    H5a1

    Normandie Netherlands Friesland Finland Orkney
    Generalissimo, teepean47, I had completely forgotten this possible issue with BAM Analysis Kit. The results of such an experiment may be very interesting. Very useful job in any case.
    En North alom, de North venom
    En North fum naiz, en North manom

    (Roman de Rou, Wace, 1160-1170)

  2. The Following User Says Thank You to anglesqueville For This Useful Post:

     lukaszM (11-05-2018)

  3. #22
    Registered Users
    Posts
    20
    Sex

    Quote Originally Posted by anglesqueville View Post
    Generalissimo, teepean47, I had completely forgotten this possible issue with BAM Analysis Kit. The results of such an experiment may be very interesting. Very useful job in any case.
    I agree. The only comparison I have done so far is to see if haplogroups match with those from the scientific papers.

  4. #23
    Gold Member Class
    Posts
    5,164
    Sex

    Quote Originally Posted by Generalissimo View Post
    Download some samples from the Reich Lab, and see how they compare to your versions. Not very difficult.

    If you can't go through such a basic step to validate your pipeline output, then don't bother with it.

    Of course you're aware that calling accurate diploid genotypes from low quality (even less than 20x coverage) data is very difficult, hence, in academic works low coverage data are called as pseudo-haploid genotypes.

    Here's an example of one way around that. But if you don't have a way around it, then there's no point using your pipeline.



    https://github.com/kveeramah/aDNA_GenoCaller
    In our world ( aDNA) 20X average read depth is certainly not low quality, in fact it's exceptional. I have diploid genotyped aDNA with even 2.5X coverage (Scyth-I0247-Volga) successfully using various GATK pipelines (both Unified Genotyper and Haplotype Caller) and of course the ATLAS pipeline.

    I'll take a minute to clear up alot of the confusion out there regarding GQ and low coverage aDNA, but 1st outline the disadvantage of using pseudo-haploid genotypes, not that they should never be used, because sometimes there is isn't much of a choice such as with aDNA with say 1X coverage where you have a very large number of positions supported by 1 read or less.

    First a word about pseudo haploids:

    Pseudo-haploids work fine in programs such as ADMIXTURE and I believe ADMIXTOOLS (it's been a while since I compared dipoloids and pseudo-haploids using ADMIXTOOLS) because ADMIXTURE relies on only 1 allele in the pair, the one it determines to be the minor allele. For example, say you are looking at a position with an actual genotype of AG in the aDNA, but is pseudo-haploid genotyped as GG. As long as G is in fact the minor allele, the inferred admixture assignment at that site will be correct. That is the reason why you don't need correct phase information either in ADMIXTURE. So for example say you have a sample with incorrect phase of GA, the ADMIXTURE output will still be correct. Of course not the case for IBD or IBS analysis because here pseudo-haploids will give you very inaccurate results.

    PROOF:

    Take a grandparent and grandchild genome. Leave the grandparent as diploid, and transform the grandchild into pseudo-haploid by splicing its genome (one can easily write a script in Linux to do this). Upload the grandparent and both versions of the grandchild (pseudo-haploid and diploid) to GEDMATCH and do 1 to 1 comparison. I already did this. Here is the result:

    1- Grandparent - grandchild -diploid - largest seg 188.9cM, total seg 1707cM
    2- Grandparent - grandchild -pseudohaploid - largest seg 38.7cM, total seg 68.7cM

    Pseudo-Haploids are about 30% wrong because we are hetrozygous at about 30% of our genome.



    ANOTHER PROOF:

    I compared IBS on Scythian-I0247-Volga (2.5X coverage) using the pseudo-haploid academic sample with one that I diploid genotyped using ATLAS. Results are on my website at http://www.eurasiandna.com/2017/10/0...e-ancient-dna/ various figures including fig11 and fig12. Notice with the academic pseudo-haploid version Swedes show greater allele sharing with it than Tatars, and Sardinians greater than Estonians, and S Iranians less than Bedouins, whereas for the diploid version I genotyped the opposite is true which of course makes more sense.
    EurasianDNA.com - A study of the population history of West & South Asia.

  5. The Following 3 Users Say Thank You to Kurd For This Useful Post:

     lukaszM (11-05-2018),  Michał (11-05-2018),  Pribislav (11-05-2018)

  6. #24
    Gold Member Class
    Posts
    5,164
    Sex

    Duplicate
    Last edited by Kurd; 11-05-2018 at 12:48 PM.
    EurasianDNA.com - A study of the population history of West & South Asia.

  7. #25
    Gold Member Class
    Posts
    5,164
    Sex

    Continued from above:


    The confusion and hype about using the 20X as a threshold for accuracy probably started with the GATK team at Broad and that is where folks such as Nick have probably gotten their info. We have to remember that the concerns at GATK center around the work they do with rare diseases and cancers. In those cases it is critical to get the genotype correct to diagnose and treat, thus they require a very high amount of evidence for the genotype which in turn means high GQ and high coverage.

    In our world we deal with junk DNA, thus is not that important that our called genotypes be 99.99% accurate. Also, the published coverages are AVERAGE coverages and not INDIVIDUAL coverages at various positions. Thus you can have a published coverage of 2.5X and many positions covered by 5 or 10 reads. I for example can easily write software that identifies positions with say only 1 read depth and either filter those out, or changing them to pseudo-haploid individually without changing the whole genome to pseudo-haploid.




    This version of the program will take any genotype with a low quality heterozygote call (Q<30) and convert to the next best homozygote call.
    This to me wouldn't make sense to take a diploid genome with say 2X coverage and plenty of positions with say read depths of 5 or more and turn it into a pseudo-haploid.

    For everyone's benefit on this forum, here is some clarification on GQ. This is a snapshot of a VCF of academic pseudo-haploid ancients ; EHG Poltavka_I0440 Scythian_I0247 Srubnaya_I0232 Srubnaya_I0430 Yamna_I0443 Yamna_IO231 which I diploid genotyped using the GATK HaplotypeCaller pipeline with VQSR:

    #CHROM POS ID REF ALT QUAL FILTER FORMAT EHG Poltavka_I0440 Scythian_I0247 Srubnaya_I0232 Srubnaya_I0430 Yamna_I0443 Yamna_IO231
    1 752566 rs3094315 G A 2062.69 PASS GT:ADP:GQL 1/1:0,16:16:47:455,47,0 0/1:9,3:12:43:43,0,283 0/1:5,5:10:99:126,0,121 0/1:12,3:15:44:44,0,363 1/1:0,3:3:9:84,9,0 1/1:0,23:23:68:675,68,0 0/1:9,27:36:99:678,0,196
    1 797440 rs58013264 T C 42.08 PASS GT:ADP:GQL 0/0:2,0:2:6:0,6,61 0/0:1,0:1:3:0,3,30 0/0:1,0:1:3:0,3,29 1/1:0,2:2:6:73,6,0 ./.:0,0:0:.:0,0,0 ./.:0,0:0:.:0,0,0 ./.:0,0:0:.:0,0,0
    1 1648831 rs61777460 T
    C 1639.77 PASS GT:ADP:GQL 0/1:1,8:9:8:255,0,8 0/1:2,11:13:30:391,0,30 0/1:6,1:7:16:16,0,174 0/1:3,17:20:38:611,0,38 ./.:5,0:5:.:0,0,0 0/1:8,3:11:77:77,0,228 1/1:0,10:10:30:328,30,0



    Focus on the EHG 3rd line for now. The GT (genotype) is 0/1 meaning hetro site T/C. The DP (depth) is 9 with 8 C reads (variant) and 1 T reads (REF) shown by AD. GQ is 8 which in reality means that the next likely genotype for this position is CC shown by the number 8 in the PL field. The least likely genotype is HOM-REF (A/A) shown by 255 in the PL field. So what does the GQ of 8 mean. It means that the likelihood of CC being correct in lieu of our genotyped TC is 10^-0.8 which is 0.158 or 15.8% which is still a better position to be in than going with a pseudo-haploid call of CC or even worse TT, which has a likelihood of only 10^-25.5 which is infinitesimal

    EDIT: Blue smiley is DP and yellow one is PL
    Last edited by Kurd; 11-05-2018 at 12:50 PM.
    EurasianDNA.com - A study of the population history of West & South Asia.

  8. The Following 2 Users Say Thank You to Kurd For This Useful Post:

     lukaszM (11-05-2018),  Michał (11-05-2018)

  9. #26
    Registered Users
    Posts
    2,705
    Sex

    Quote Originally Posted by teepean47 View Post
    Sure. Do you have any preference in output format or options or do you want whatever plink outputs?
    PLINK binary (bed, bim, fam).

  10. #27
    Registered Users
    Posts
    20
    Sex

    As a first test sample I chose I1165 from Harney et al. 2018.

    https://reich.hms.harvard.edu/datasets

    For mtDNA I used a separate sample. Both BAMs can be downloaded from:

    https://www.ebi.ac.uk/ena/data/view/ERS2663520

    Haplogroup results using BAM Analysis Kit + haplogrep and yhaplo:

    mtDNA haplogroup:

    mtDNA haplogroup from the first sample:
    U5b3d

    And from the second sample:
    HV1a'b'c

    mtDNA haplogroup according to haplogrep:

    SampleID Range Haplogroup Overall_Rank
    sample 1-16569; HV1 0.5399

    Separate mtDNA sample:
    SampleID Range Haplogroup Overall_Rank
    sample 1-16569; H2a2a1 0.5000

    Y haplogroup according to yhaplo. BAM Analysis Kit didn't give any meaningful result.
    sample-sample T-CTS2214 T-CTS2214 T1a1a1
    PLINK files:

    https://drive.google.com/open?id=1YG...BSMeKBzTEwwCnX
    Last edited by teepean47; 11-08-2018 at 08:27 PM. Reason: Updated Plink files

  11. The Following 6 Users Say Thank You to teepean47 For This Useful Post:

     anglesqueville (11-08-2018),  Generalissimo (11-08-2018),  lukaszM (11-09-2018),  R.Rocca (11-08-2018),  Robert1 (11-09-2018),  vettor (11-10-2018)

  12. #28
    Registered Users
    Posts
    2,705
    Sex

    Quote Originally Posted by teepean47 View Post
    Where's the bed file?

  13. The Following User Says Thank You to Generalissimo For This Useful Post:

     anglesqueville (11-08-2018)

  14. #29
    Registered Users
    Posts
    20
    Sex

    Quote Originally Posted by Generalissimo View Post
    Where's the bed file?
    https://drive.google.com/open?id=1YG...BSMeKBzTEwwCnX

  15. The Following User Says Thank You to teepean47 For This Useful Post:

     anglesqueville (11-08-2018)

  16. #30
    Registered Users
    Posts
    2,705
    Sex

    OK, my version of PLINK can't run that bim file because it has identical alleles on the same lines.

    So, can you upload a 23andMe version of the same files?

Page 3 of 5 FirstFirst 12345 LastLast

Similar Threads

  1. 23andme YDNA Hg/SNP of Interest
    By asquecco in forum 23andMe
    Replies: 3
    Last Post: 11-13-2018, 09:15 PM
  2. Replies: 51
    Last Post: 09-03-2017, 10:55 PM
  3. South Asian? This may be of interest.
    By Amerijoe in forum Other
    Replies: 4
    Last Post: 07-17-2017, 04:23 AM
  4. What motivates your interest?
    By A Norfolk L-M20 in forum General
    Replies: 48
    Last Post: 10-13-2016, 02:17 AM
  5. Replies: 2
    Last Post: 04-16-2015, 01:29 AM

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •