Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Accessing Genomic data in a large vcf.gz file

    Hello all, I am very new to bioinformatics; I'm a high school intern. I have a .vcf.gz file that is 332 GB. It contains the whole genomes of around 800 dogs and all I need to do is check for a mutation for each of them. But before I can do that, I need to figure out how to crack the file open. I have heard of ways to possibly access them without decompressing the .vcf.gz because someone told me its possible that the uncompressed file could be up to three terabytes.

    If anyone has any suggestions on how to proceed, I would be eternally grateful as I am nearing the end of my internship and need to figure this out.

  • #2
    I think "bcftools view" command should give you what you need. Check the help page here.

    Comment


    • #3
      Canine genomes

      I liked GenoMax's suggestion, so go with that if you haven't already. Also be sure to use long read (PacBio) based assembly for your analysis against a reference genome.
      If you are interested in canine genomes, here are some articles that may be of interest:

      Enjoy the videos and music you love, upload original content, and share it all with friends, family, and the world on YouTube.

      A proposal to sequence the genome of the German Shepherd dog by the University of Wisconsin - Madison was selected as winner of the 2019 Plant and Animal SMRT Grant after garnering more than 3,000 votes in a close online video competition. We caught up with members of the Comparative Genetics Research Laboratory to find out more about the project.

      Scientists are doing deep dives into the genomes of a range of canine cousins along the evolutionary chain.

      #SeqLife

      Comment


      • #4
        Originally posted by GenoMax View Post
        I think "bcftools view" command should give you what you need. Check the help page here.
        Thanks so much, I'll try this out

        Comment


        • #5
          Originally posted by GenoMax View Post
          I think "bcftools view" command should give you what you need. Check the help page here.
          So I have bcftools installed. I want to check for a single-replacement mutation on a specific gene in each of the genomes in the file. Do you know of a way to do this using bcftools? I am new to bioinformatics and this is going way over my head. Thanks

          Comment


          • #6
            Originally posted by GenoMax View Post
            I think "bcftools view" command should give you what you need. Check the help page here.
            Note that for very large vcf files, using -t/-T instead of -r/-R is going to be much faster. And use the --threads option, that will help as well.

            Comment


            • #7
              Originally posted by questor2010 View Post
              Note that for very large vcf files, using -t/-T instead of -r/-R is going to be much faster. And use the --threads option, that will help as well.
              Thanks for your reply. Could you give me an example of the format of what a command that would use bcftools view to find a base at a position would look like?

              Comment


              • #8
                Originally posted by XeroxHero69 View Post
                Thanks for your reply. Could you give me an example of the format of what a command that would use bcftools view to find a base at a position would look like?
                Sure - here's an example:

                bcftools view -f PASS --threads 8 -T target.bed -o gnomad.genomes.r2.1.sites.target.vcf.gz -O z gnomad.genomes.r2.1.sites.vcf.bgz

                In this case, I'm using a bed file, instead of a single region. It is pulling out the FILTER=PASS variants that intersect the bed file into a new compressed vcf file. The source vcf file in this case is 465GB. If you have a single variant, you could use -t 1:11022. It might be best to specify a short range (1:11015-11030) if you're looking at indels - variant callers represent indels in different ways and you want to be sure you properly intersect them.

                Comment


                • #9
                  Originally posted by questor2010 View Post
                  Sure - here's an example:

                  bcftools view -f PASS --threads 8 -T target.bed -o gnomad.genomes.r2.1.sites.target.vcf.gz -O z gnomad.genomes.r2.1.sites.vcf.bgz

                  In this case, I'm using a bed file, instead of a single region. It is pulling out the FILTER=PASS variants that intersect the bed file into a new compressed vcf file. The source vcf file in this case is 465GB. If you have a single variant, you could use -t 1:11022. It might be best to specify a short range (1:11015-11030) if you're looking at indels - variant callers represent indels in different ways and you want to be sure you properly intersect them.
                  Thanks so much. I tried this command:

                  bcftools view -f PASS --threads 8 -r chr9:55252802-55252802 -o 722g.990.SNP.INDEL.chrAll.vcf.gz -O z 722g.990.SNP.INDEL.chrAll.vcf.gz

                  and it returned:
                  [W::hts_idx_load2] The index file is older than the data file: 722g.990.SNP.INDEL.chrAll.vcf.gz.tbi
                  [W::hts_idx_load2] The index file is older than the data file: 722g.990.SNP.INDEL.chrAll.vcf.gz.tbi
                  [W::bgzf_read_block] EOF marker is absent. The input is probably truncated

                  Do you think you could point out what I did wrong?
                  Edit: I think i see my error, i made it output the data to the same file and now I think i ruined that data file... its only 16.1 kb now
                  Last edited by XeroxHero69; 05-01-2019, 10:38 AM.

                  Comment


                  • #10
                    Originally posted by XeroxHero69 View Post
                    Thanks so much. I tried this command:

                    bcftools view -f PASS --threads 8 -r chr9:55252802-55252802 -o 722g.990.SNP.INDEL.chrAll.vcf.gz -O z 722g.990.SNP.INDEL.chrAll.vcf.gz

                    and it returned:
                    [W::hts_idx_load2] The index file is older than the data file: 722g.990.SNP.INDEL.chrAll.vcf.gz.tbi
                    [W::hts_idx_load2] The index file is older than the data file: 722g.990.SNP.INDEL.chrAll.vcf.gz.tbi
                    [W::bgzf_read_block] EOF marker is absent. The input is probably truncated

                    Do you think you could point out what I did wrong?
                    Edit: I think i see my error, i made it output the data to the same file and now I think i ruined that data file... its only 16.1 kb now
                    Ouch. Yes - that's the problem. I hope you have another copy readily available. The index file warnings are often caused by file transfers. They transfer first for large files and then have an older timestamp. You can use "touch" to fix that, or reindex, but that takes a long time with large vcf files.

                    Comment


                    • #11
                      Originally posted by questor2010 View Post
                      Ouch. Yes - that's the problem. I hope you have another copy readily available. The index file warnings are often caused by file transfers. They transfer first for large files and then have an older timestamp. You can use "touch" to fix that, or reindex, but that takes a long time with large vcf files.
                      So to make sure I do it right this time, what do I do to make the .vcf.bgz file?

                      Comment


                      • #12
                        The -o switch specifies the output file name, the -O switch specifies the format:

                        -O, --output-type <b|u|z|v> b: compressed BCF, u: uncompressed BCF, z: compressed VCF, v: uncompressed VCF

                        So -O z implies that the name should be basename.vcf.gz (or basename.vcf.bgz).

                        If you have a single variant you just want to examine, you could use -v and output to a standard vcf file (text).

                        Comment


                        • #13
                          Originally posted by questor2010 View Post
                          The -o switch specifies the output file name, the -O switch specifies the format:

                          -O, --output-type <b|u|z|v> b: compressed BCF, u: uncompressed BCF, z: compressed VCF, v: uncompressed VCF

                          So -O z implies that the name should be basename.vcf.gz (or basename.vcf.bgz).

                          If you have a single variant you just want to examine, you could use -v and output to a standard vcf file (text).
                          So do I need to create a new file for the output prior to running this command? how would the command look if I use -v without overwriting my file again. Sorry if these are simple questions. Thanks!

                          Comment

                          Latest Articles

                          Collapse

                          • seqadmin
                            Strategies for Sequencing Challenging Samples
                            by seqadmin


                            Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
                            03-22-2024, 06:39 AM
                          • seqadmin
                            Techniques and Challenges in Conservation Genomics
                            by seqadmin



                            The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

                            Avian Conservation
                            Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
                            03-08-2024, 10:41 AM

                          ad_right_rmr

                          Collapse

                          News

                          Collapse

                          Topics Statistics Last Post
                          Started by seqadmin, Yesterday, 06:37 PM
                          0 responses
                          10 views
                          0 likes
                          Last Post seqadmin  
                          Started by seqadmin, Yesterday, 06:07 PM
                          0 responses
                          9 views
                          0 likes
                          Last Post seqadmin  
                          Started by seqadmin, 03-22-2024, 10:03 AM
                          0 responses
                          49 views
                          0 likes
                          Last Post seqadmin  
                          Started by seqadmin, 03-21-2024, 07:32 AM
                          0 responses
                          67 views
                          0 likes
                          Last Post seqadmin  
                          Working...
                          X