Accessing Genomic data in a large vcf.gz file

[XeroxHero69]

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]

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

[JTH-Genome]

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=OjsvbDqfp7w&feature=youtu.be

https://www.pacb.com/blog/german-shepherd-genome/

https://www.pacb.com/blog/dog-meet-dog-world-exploring-canine-genomes/

#SeqLife

[XeroxHero69]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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!