Selecting 65000 SNPs where AF is close to 0.5 in all or most populations












2












$begingroup$


I am evaluating the tool somalier (https://github.com/brentp/somalier) and I need to create a list of about 65,000 SNPs where the allele frequency (AF) is as close to 0.5 as possible across the most representative set of populations possible with today's data.



My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is available in the somalier website, which is in GRCh37.



I am starting as input with the vcfs from the 1000genomes project in GRCh38 positions, and I would like to come up with a filtering strategy that:
(1) Sorts all SNPs by the AF being close to 0.5 in as many populations as possible.
(2) Takes the first 65,000 SNPs from that list and produces a vcf file from them.



My simplistic strategy so far has been to take the global AF value if it's between 0.3-0.7 and print it if it's a SNP with this perl one-liner:



gunzip -c my.vcf.gz | grep -v '#' | perl -lne 'if ($_ =~ /;AF=0.[34567]/) { if ($_ =~ /w+td+tS+twtwt/) {print $_} }'


Any recommendations?






vcf snp genomics variation







share|improve this question











$endgroup$












  • $begingroup$
    If possible, I strongly urge you to use hg19 instead of GRCh38 because there are much better frequency data available. Basically, GnomAD uses hg19, so any gnomAD frequencies for hg38 will come from lifting over the positions and that doesn't work too well in my experience.
    $endgroup$
    – terdon♦
    3 hours ago










  • $begingroup$
    I haven't been able to reliably liftover the sites.vcf.gz file that can be downloaded from the github.com/brentp/somalier site. I will create another Q for the liftover issue.
    $endgroup$
    – 719016
    2 hours ago










  • $begingroup$
    What do you mean? Somalier's github page says it's for hg19, not hg38: "where sites is a VCF of variant sites (provided by somalier for hg19)."
    $endgroup$
    – terdon♦
    2 hours ago










  • $begingroup$
    My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is in GRCh37 (gh19).
    $endgroup$
    – 719016
    2 hours ago
















2












$begingroup$


I am evaluating the tool somalier (https://github.com/brentp/somalier) and I need to create a list of about 65,000 SNPs where the allele frequency (AF) is as close to 0.5 as possible across the most representative set of populations possible with today's data.



My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is available in the somalier website, which is in GRCh37.



I am starting as input with the vcfs from the 1000genomes project in GRCh38 positions, and I would like to come up with a filtering strategy that:
(1) Sorts all SNPs by the AF being close to 0.5 in as many populations as possible.
(2) Takes the first 65,000 SNPs from that list and produces a vcf file from them.



My simplistic strategy so far has been to take the global AF value if it's between 0.3-0.7 and print it if it's a SNP with this perl one-liner:



gunzip -c my.vcf.gz | grep -v '#' | perl -lne 'if ($_ =~ /;AF=0.[34567]/) { if ($_ =~ /w+td+tS+twtwt/) {print $_} }'


Any recommendations?






vcf snp genomics variation







share|improve this question











$endgroup$












  • $begingroup$
    If possible, I strongly urge you to use hg19 instead of GRCh38 because there are much better frequency data available. Basically, GnomAD uses hg19, so any gnomAD frequencies for hg38 will come from lifting over the positions and that doesn't work too well in my experience.
    $endgroup$
    – terdon♦
    3 hours ago










  • $begingroup$
    I haven't been able to reliably liftover the sites.vcf.gz file that can be downloaded from the github.com/brentp/somalier site. I will create another Q for the liftover issue.
    $endgroup$
    – 719016
    2 hours ago










  • $begingroup$
    What do you mean? Somalier's github page says it's for hg19, not hg38: "where sites is a VCF of variant sites (provided by somalier for hg19)."
    $endgroup$
    – terdon♦
    2 hours ago










  • $begingroup$
    My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is in GRCh37 (gh19).
    $endgroup$
    – 719016
    2 hours ago














2












2








2





$begingroup$


I am evaluating the tool somalier (https://github.com/brentp/somalier) and I need to create a list of about 65,000 SNPs where the allele frequency (AF) is as close to 0.5 as possible across the most representative set of populations possible with today's data.



My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is available in the somalier website, which is in GRCh37.



I am starting as input with the vcfs from the 1000genomes project in GRCh38 positions, and I would like to come up with a filtering strategy that:
(1) Sorts all SNPs by the AF being close to 0.5 in as many populations as possible.
(2) Takes the first 65,000 SNPs from that list and produces a vcf file from them.



My simplistic strategy so far has been to take the global AF value if it's between 0.3-0.7 and print it if it's a SNP with this perl one-liner:



gunzip -c my.vcf.gz | grep -v '#' | perl -lne 'if ($_ =~ /;AF=0.[34567]/) { if ($_ =~ /w+td+tS+twtwt/) {print $_} }'


Any recommendations?






vcf snp genomics variation







share|improve this question











$endgroup$




I am evaluating the tool somalier (https://github.com/brentp/somalier) and I need to create a list of about 65,000 SNPs where the allele frequency (AF) is as close to 0.5 as possible across the most representative set of populations possible with today's data.



My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is available in the somalier website, which is in GRCh37.



I am starting as input with the vcfs from the 1000genomes project in GRCh38 positions, and I would like to come up with a filtering strategy that:
(1) Sorts all SNPs by the AF being close to 0.5 in as many populations as possible.
(2) Takes the first 65,000 SNPs from that list and produces a vcf file from them.



My simplistic strategy so far has been to take the global AF value if it's between 0.3-0.7 and print it if it's a SNP with this perl one-liner:



gunzip -c my.vcf.gz | grep -v '#' | perl -lne 'if ($_ =~ /;AF=0.[34567]/) { if ($_ =~ /w+td+tS+twtwt/) {print $_} }'


Any recommendations?






vcf snp genomics variation




vcf snp genomics variation






share|improve this question















share|improve this question













share|improve this question




share|improve this question








edited 2 hours ago







719016

















asked 3 hours ago









719016719016

1,029315




1,029315












  • $begingroup$
    If possible, I strongly urge you to use hg19 instead of GRCh38 because there are much better frequency data available. Basically, GnomAD uses hg19, so any gnomAD frequencies for hg38 will come from lifting over the positions and that doesn't work too well in my experience.
    $endgroup$
    – terdon♦
    3 hours ago










  • $begingroup$
    I haven't been able to reliably liftover the sites.vcf.gz file that can be downloaded from the github.com/brentp/somalier site. I will create another Q for the liftover issue.
    $endgroup$
    – 719016
    2 hours ago










  • $begingroup$
    What do you mean? Somalier's github page says it's for hg19, not hg38: "where sites is a VCF of variant sites (provided by somalier for hg19)."
    $endgroup$
    – terdon♦
    2 hours ago










  • $begingroup$
    My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is in GRCh37 (gh19).
    $endgroup$
    – 719016
    2 hours ago


















  • $begingroup$
    If possible, I strongly urge you to use hg19 instead of GRCh38 because there are much better frequency data available. Basically, GnomAD uses hg19, so any gnomAD frequencies for hg38 will come from lifting over the positions and that doesn't work too well in my experience.
    $endgroup$
    – terdon♦
    3 hours ago










  • $begingroup$
    I haven't been able to reliably liftover the sites.vcf.gz file that can be downloaded from the github.com/brentp/somalier site. I will create another Q for the liftover issue.
    $endgroup$
    – 719016
    2 hours ago










  • $begingroup$
    What do you mean? Somalier's github page says it's for hg19, not hg38: "where sites is a VCF of variant sites (provided by somalier for hg19)."
    $endgroup$
    – terdon♦
    2 hours ago










  • $begingroup$
    My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is in GRCh37 (gh19).
    $endgroup$
    – 719016
    2 hours ago
















$begingroup$
If possible, I strongly urge you to use hg19 instead of GRCh38 because there are much better frequency data available. Basically, GnomAD uses hg19, so any gnomAD frequencies for hg38 will come from lifting over the positions and that doesn't work too well in my experience.
$endgroup$
– terdon♦
3 hours ago




$begingroup$
If possible, I strongly urge you to use hg19 instead of GRCh38 because there are much better frequency data available. Basically, GnomAD uses hg19, so any gnomAD frequencies for hg38 will come from lifting over the positions and that doesn't work too well in my experience.
$endgroup$
– terdon♦
3 hours ago












$begingroup$
I haven't been able to reliably liftover the sites.vcf.gz file that can be downloaded from the github.com/brentp/somalier site. I will create another Q for the liftover issue.
$endgroup$
– 719016
2 hours ago




$begingroup$
I haven't been able to reliably liftover the sites.vcf.gz file that can be downloaded from the github.com/brentp/somalier site. I will create another Q for the liftover issue.
$endgroup$
– 719016
2 hours ago












$begingroup$
What do you mean? Somalier's github page says it's for hg19, not hg38: "where sites is a VCF of variant sites (provided by somalier for hg19)."
$endgroup$
– terdon♦
2 hours ago




$begingroup$
What do you mean? Somalier's github page says it's for hg19, not hg38: "where sites is a VCF of variant sites (provided by somalier for hg19)."
$endgroup$
– terdon♦
2 hours ago












$begingroup$
My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is in GRCh37 (gh19).
$endgroup$
– 719016
2 hours ago




$begingroup$
My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is in GRCh37 (gh19).
$endgroup$
– 719016
2 hours ago










2 Answers
2






active

oldest

votes


















4












$begingroup$

Just use bcftools view for filtering:



$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz > output.vcf


To truncate this list to 65,000 SNPs count the header lines, sum them to the number of SNPs you like to have and use head -n.



$ bcftools view -h input.vcf.gz|wc -l

255



$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz | head -n 65255 > output.vcf





share|improve this answer









$endgroup$





















    3












    $begingroup$

    Using bcftools is the right way, but you could also greatly simplify the command you're using. You can do the same thing with a single zgrep call:



    zgrep -P '^[^#]+td+tS+twtwt.*;AF=0.[3-7]' my.vcf.gz


    And to limit to the first 65000 results:



    zgrep -m 65000 -P '^[^#]+td+tS+twtwt.*;AF=0.[3-7]' my.vcf.gz


    Finally, to keep the header of the original vcf unchanged, you could do:



    zgrep -m 65000 -P '^(#|[^#]+td+tS+twtwt.*;AF=0.[3-7])' my.vcf.gz


    I also suggest you use hg19 instead of hg38 data since the best available source for frequencies, the gnomAD database, is based on hg19 and the tool you want to test also provides a site file for hg19. But that may not be possible for you.



    If your VCF file has multiple AF entries for different sub-populations and you want to use those, you'll probably need to script something. We can't really help with that unless you show us some example input lines and start writing the basic script yourself though.






    share|improve this answer









    $endgroup$













      Your Answer





      StackExchange.ifUsing("editor", function () {
      return StackExchange.using("mathjaxEditing", function () {
      StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
      StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
      });
      });
      }, "mathjax-editing");

      StackExchange.ready(function() {
      var channelOptions = {
      tags: "".split(" "),
      id: "676"
      };
      initTagRenderer("".split(" "), "".split(" "), channelOptions);

      StackExchange.using("externalEditor", function() {
      // Have to fire editor after snippets, if snippets enabled
      if (StackExchange.settings.snippets.snippetsEnabled) {
      StackExchange.using("snippets", function() {
      createEditor();
      });
      }
      else {
      createEditor();
      }
      });

      function createEditor() {
      StackExchange.prepareEditor({
      heartbeatType: 'answer',
      autoActivateHeartbeat: false,
      convertImagesToLinks: false,
      noModals: true,
      showLowRepImageUploadWarning: true,
      reputationToPostImages: null,
      bindNavPrevention: true,
      postfix: "",
      imageUploader: {
      brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
      contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
      allowUrls: true
      },
      onDemand: true,
      discardSelector: ".discard-answer"
      ,immediatelyShowMarkdownHelp:true
      });


      }
      });














      draft saved

      draft discarded


















      StackExchange.ready(
      function () {
      StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fbioinformatics.stackexchange.com%2fquestions%2f6988%2fselecting-65000-snps-where-af-is-close-to-0-5-in-all-or-most-populations%23new-answer', 'question_page');
      }
      );

      Post as a guest















      Required, but never shown

























      2 Answers
      2






      active

      oldest

      votes








      2 Answers
      2






      active

      oldest

      votes









      active

      oldest

      votes






      active

      oldest

      votes









      4












      $begingroup$

      Just use bcftools view for filtering:



      $ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz > output.vcf


      To truncate this list to 65,000 SNPs count the header lines, sum them to the number of SNPs you like to have and use head -n.



      $ bcftools view -h input.vcf.gz|wc -l
      
255
      

      
$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz | head -n 65255 > output.vcf





      share|improve this answer









      $endgroup$


















        4












        $begingroup$

        Just use bcftools view for filtering:



        $ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz > output.vcf


        To truncate this list to 65,000 SNPs count the header lines, sum them to the number of SNPs you like to have and use head -n.



        $ bcftools view -h input.vcf.gz|wc -l
        
255
        

        
$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz | head -n 65255 > output.vcf





        share|improve this answer









        $endgroup$
















          4












          4








          4





          $begingroup$

          Just use bcftools view for filtering:



          $ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz > output.vcf


          To truncate this list to 65,000 SNPs count the header lines, sum them to the number of SNPs you like to have and use head -n.



          $ bcftools view -h input.vcf.gz|wc -l
          
255
          

          
$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz | head -n 65255 > output.vcf





          share|improve this answer









          $endgroup$



          Just use bcftools view for filtering:



          $ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz > output.vcf


          To truncate this list to 65,000 SNPs count the header lines, sum them to the number of SNPs you like to have and use head -n.



          $ bcftools view -h input.vcf.gz|wc -l
          
255
          

          
$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz | head -n 65255 > output.vcf






          share|improve this answer












          share|improve this answer



          share|improve this answer










          answered 2 hours ago









          finswimmerfinswimmer

          84229




          84229























              3












              $begingroup$

              Using bcftools is the right way, but you could also greatly simplify the command you're using. You can do the same thing with a single zgrep call:



              zgrep -P '^[^#]+td+tS+twtwt.*;AF=0.[3-7]' my.vcf.gz


              And to limit to the first 65000 results:



              zgrep -m 65000 -P '^[^#]+td+tS+twtwt.*;AF=0.[3-7]' my.vcf.gz


              Finally, to keep the header of the original vcf unchanged, you could do:



              zgrep -m 65000 -P '^(#|[^#]+td+tS+twtwt.*;AF=0.[3-7])' my.vcf.gz


              I also suggest you use hg19 instead of hg38 data since the best available source for frequencies, the gnomAD database, is based on hg19 and the tool you want to test also provides a site file for hg19. But that may not be possible for you.



              If your VCF file has multiple AF entries for different sub-populations and you want to use those, you'll probably need to script something. We can't really help with that unless you show us some example input lines and start writing the basic script yourself though.






              share|improve this answer









              $endgroup$


















                3












                $begingroup$

                Using bcftools is the right way, but you could also greatly simplify the command you're using. You can do the same thing with a single zgrep call:



                zgrep -P '^[^#]+td+tS+twtwt.*;AF=0.[3-7]' my.vcf.gz


                And to limit to the first 65000 results:



                zgrep -m 65000 -P '^[^#]+td+tS+twtwt.*;AF=0.[3-7]' my.vcf.gz


                Finally, to keep the header of the original vcf unchanged, you could do:



                zgrep -m 65000 -P '^(#|[^#]+td+tS+twtwt.*;AF=0.[3-7])' my.vcf.gz


                I also suggest you use hg19 instead of hg38 data since the best available source for frequencies, the gnomAD database, is based on hg19 and the tool you want to test also provides a site file for hg19. But that may not be possible for you.



                If your VCF file has multiple AF entries for different sub-populations and you want to use those, you'll probably need to script something. We can't really help with that unless you show us some example input lines and start writing the basic script yourself though.






                share|improve this answer









                $endgroup$
















                  3












                  3








                  3





                  $begingroup$

                  Using bcftools is the right way, but you could also greatly simplify the command you're using. You can do the same thing with a single zgrep call:



                  zgrep -P '^[^#]+td+tS+twtwt.*;AF=0.[3-7]' my.vcf.gz


                  And to limit to the first 65000 results:



                  zgrep -m 65000 -P '^[^#]+td+tS+twtwt.*;AF=0.[3-7]' my.vcf.gz


                  Finally, to keep the header of the original vcf unchanged, you could do:



                  zgrep -m 65000 -P '^(#|[^#]+td+tS+twtwt.*;AF=0.[3-7])' my.vcf.gz


                  I also suggest you use hg19 instead of hg38 data since the best available source for frequencies, the gnomAD database, is based on hg19 and the tool you want to test also provides a site file for hg19. But that may not be possible for you.



                  If your VCF file has multiple AF entries for different sub-populations and you want to use those, you'll probably need to script something. We can't really help with that unless you show us some example input lines and start writing the basic script yourself though.






                  share|improve this answer









                  $endgroup$



                  Using bcftools is the right way, but you could also greatly simplify the command you're using. You can do the same thing with a single zgrep call:



                  zgrep -P '^[^#]+td+tS+twtwt.*;AF=0.[3-7]' my.vcf.gz


                  And to limit to the first 65000 results:



                  zgrep -m 65000 -P '^[^#]+td+tS+twtwt.*;AF=0.[3-7]' my.vcf.gz


                  Finally, to keep the header of the original vcf unchanged, you could do:



                  zgrep -m 65000 -P '^(#|[^#]+td+tS+twtwt.*;AF=0.[3-7])' my.vcf.gz


                  I also suggest you use hg19 instead of hg38 data since the best available source for frequencies, the gnomAD database, is based on hg19 and the tool you want to test also provides a site file for hg19. But that may not be possible for you.



                  If your VCF file has multiple AF entries for different sub-populations and you want to use those, you'll probably need to script something. We can't really help with that unless you show us some example input lines and start writing the basic script yourself though.







                  share|improve this answer












                  share|improve this answer



                  share|improve this answer










                  answered 2 hours ago









                  terdon♦terdon

                  4,2211729




                  4,2211729






























                      draft saved

                      draft discarded




















































                      Thanks for contributing an answer to Bioinformatics Stack Exchange!


                      • Please be sure to answer the question. Provide details and share your research!

                      But avoid …



                      • Asking for help, clarification, or responding to other answers.

                      • Making statements based on opinion; back them up with references or personal experience.


                      Use MathJax to format equations. MathJax reference.


                      To learn more, see our tips on writing great answers.




                      draft saved


                      draft discarded














                      StackExchange.ready(
                      function () {
                      StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fbioinformatics.stackexchange.com%2fquestions%2f6988%2fselecting-65000-snps-where-af-is-close-to-0-5-in-all-or-most-populations%23new-answer', 'question_page');
                      }
                      );

                      Post as a guest















                      Required, but never shown





















































                      Required, but never shown














                      Required, but never shown












                      Required, but never shown







                      Required, but never shown

































                      Required, but never shown














                      Required, but never shown












                      Required, but never shown







                      Required, but never shown







                      -

                      Popular posts from this blog

                      Callistus I

                      Image
                      Callistus I Res apud Vicidata repertae: Nativitas : unknown value ; Roma Obitus : 223 ; Roma Pontificatus 217–222 Praecessor : Zephyrinus Successor : Urbanus I Sanctus Callistus qui ieiunium instituit; e libro manu scripto Francogallico saeculo XIV Callistus (natus die ignoto; obiit anno 222), fuit episcopus Romae et Papa Ecclesiae Catholicae Romanae ab anno 217. Vita | Secundum adversarium suum, presbyterem, fuit Hippolytus Callistus, servus argentarii Christiani. Papae Zephyrino, Callistus Romae arcessivit. Callistus fuit primus diaconus ac dux catacumbae apud Viam Appiam, quae hodie Sancti Callisti appellatur. Pontificatus | Nexus externi | Vicimedia Communia plura habent quae ad Callistus I spectant. Fontes Latinae de papis usque ad annum 530 (Papa Felix IV) Liber pontificalis Opera Omnia .mw-parser-output .stipula{padding:3px;background:#F7F8FF;border:1px solid grey;margin:auto}.mw-parser-output .stipula td.cell1{bac...
                      Read more

                      Tabula Rosettana

                      Image
                      Tabula multilinguis Rosettana in Museo Britannico ostenditur. Tabula Rosettana, [1] etiam titulo OGIS 90 agnita, est stela decreto de rebus sacris in Aegypto anno 196 a.C.n. lato inscripta. Tabula iuxta Rosettam Aegypti, urbem in delta Nili et ad oram maris Mediterranei iacentem, anno 1799 a milite Francico reperta est. Inventio stelae, linguis duabus scripturisque tribus inscriptae, eruditis Instituti Aegypti statim nuntiata est; ibi enim iussu imperatoris Napoleonis eruditi omnium scientiarum (sub aegide Commissionis Scientiarum et Artium) properaverant cum expeditione Francica. Qua a Britannis mox debellata, tabula Rosettana Londinium missa hodie apud Museum Britannicum iacet. Textus Graecus cito lectus interpretationi textuum Aegyptiorum (in formis hieroglyphica et demotica expressorum) gradatim adiuvit. Denique textum plene interpretatus est Ioannes Franciscus Champollion. Ab opere eruditorum cumulativo coepit hodiernus scripturae hieroglyphicae linguaeque Aegyptiae a...
                      Read more

                      How to label and detect the document text images

                      Image
                      1 $begingroup$ This is what I mean as document text image: I want to label the texts in image as separate blocks and my model should detect these labels as classes. NOTE: This is how the end result should be like: The labels like Block 1, Block 2, Block 3,.. should be Logo, Title, Date,.. Others, etc. Work done: First approach : I tried to implement this method via Object Detection, it didn't work. It didn't even detect any text. Second approach : Then I tried it using PixelLink. As this model is build for scene text detection, it detected each and every text in the image. But this method can detect multiple lines of text if the threshold values are increased. But I have no idea how do I add labels to the text blocks. PIXEL_CLS_WEIGHT_all_ones = 'PIXEL_CLS_WEIGHT_all_ones' PIXEL_C...
                      Read more