knausb / vcfR

Tools to work with variant call format files
246 stars 55 forks source link

See sequence length #213

Open jenifferteles opened 6 months ago

jenifferteles commented 6 months ago

Hello, is there any function where I can see all the information from the final VCF, like the .log file in Stacks? I would like to find out the sequence size to perform demographic analyses. Thank you.

knausb commented 6 months ago

Hi @jenifferteles

I am not familiar with th.log file from Stacks, so I can't comment on that. The sequence lengths are optionally included in VCF files. I believe this is optional, so not all VCF files have this. You should be able to check in vcfR as follows.

library(vcfR)
#> 
#>    *****       ***   vcfR   ***       *****
#>    This is vcfR 1.15.0 
#>      browseVignettes('vcfR') # Documentation
#>      citation('vcfR') # Citation
#>    *****       *****      *****       *****
data("vcfR_test")
vcfR_test
#> ***** Object of Class vcfR *****
#> 3 samples
#> 1 CHROMs
#> 5 variants
#> Object size: 0 Mb
#> 0 percent missing data
#> *****        *****         *****

queryMETA(vcfR_test)
#>  [1] "INFO=ID=NS"                          
#>  [2] "INFO=ID=DP"                          
#>  [3] "INFO=ID=AF"                          
#>  [4] "INFO=ID=AA"                          
#>  [5] "INFO=ID=DB"                          
#>  [6] "INFO=ID=H2"                          
#>  [7] "FILTER=ID=q10"                       
#>  [8] "FILTER=ID=s50"                       
#>  [9] "FORMAT=ID=GT"                        
#> [10] "FORMAT=ID=GQ"                        
#> [11] "FORMAT=ID=DP"                        
#> [12] "FORMAT=ID=HQ"                        
#> [13] "1 contig=<IDs omitted from queryMETA"
queryMETA(vcfR_test, element = "contig")
#> [[1]]
#> [1] "contig=ID=20"                        
#> [2] "length=62435964"                     
#> [3] "assembly=B36"                        
#> [4] "md5=f126cdf8a6e0c7f379d618ff66beb2da"
#> [5] "species=Homo sapiens\""              
#> [6] "taxonomy=x"

Created on 2024-03-25 with reprex v2.1.0

You may have more than one contig, so using grep to find the contig name in vcfR-test@meta may help you. Please let me know if that resolves your issue.

Brian