SEQanswers

SEQanswers (http://seqanswers.com/forums/index.php)
-   Bioinformatics (http://seqanswers.com/forums/forumdisplay.php?f=18)
-   -   Accessing Genomic data in a large vcf.gz file (http://seqanswers.com/forums/showthread.php?t=89062)

XeroxHero69 04-29-2019 11:34 AM

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.

GenoMax 04-29-2019 03:37 PM

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

JTH-Genome 04-30-2019 03:50 AM

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:
https://www.nature.com/articles/s41598-018-29190-3
https://www.youtube.com/watch?v=Ojsv...ature=youtu.be
https://www.pacb.com/blog/german-shepherd-genome/
https://www.pacb.com/blog/dog-meet-d...anine-genomes/
#SeqLife

XeroxHero69 04-30-2019 09:09 AM

Quote:

Originally Posted by GenoMax (Post 225863)
I think "bcftools view" command should give you what you need. Check the help page here.

Thanks so much, I'll try this out :)

XeroxHero69 04-30-2019 11:04 AM

Quote:

Originally Posted by GenoMax (Post 225863)
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:)

questor2010 04-30-2019 07:53 PM

Quote:

Originally Posted by GenoMax (Post 225863)
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.

XeroxHero69 05-01-2019 09:16 AM

Quote:

Originally Posted by questor2010 (Post 225892)
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?

questor2010 05-01-2019 09:35 AM

Quote:

Originally Posted by XeroxHero69 (Post 225922)
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.

XeroxHero69 05-01-2019 10:22 AM

Quote:

Originally Posted by questor2010 (Post 225924)
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

questor2010 05-01-2019 10:51 AM

Quote:

Originally Posted by XeroxHero69 (Post 225927)
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.

XeroxHero69 05-06-2019 08:58 AM

Quote:

Originally Posted by questor2010 (Post 225928)
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?

questor2010 05-06-2019 10:39 AM

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).

XeroxHero69 05-06-2019 10:44 AM

Quote:

Originally Posted by questor2010 (Post 225992)
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!


All times are GMT -8. The time now is 03:44 AM.

Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.