This function reads tabixed VCF-files, as distributed from the 1000 Genomes project (human).
readVCF(filename, numcols, tid, frompos, topos,
samplenames=NA, gffpath = FALSE, include.unknown=FALSE, approx=FALSE,
out="", parallel=FALSE)the corresponding tabixed VCF-file
number of SNPs that should be read in as a chunk
which chromosome ? (character)
start of the region
end of the region
a vector of individuals
the corresponding GFF file
includ positions with unknown/missing nucleotides
see details !
a folder suffix where the temporary files should be saved
parallel computation using mclapply
The function creates an object of class "GENOME"
---------------------------------------------------------
The following slots will be filled in the "GENOME" object
---------------------------------------------------------
| Slot | Description | |
| 1. | n.sites |
total number of sites |
| 2. | n.biallelic.sites |
number of biallelic sites |
| 3. | region.data |
some detailed information about the data read |
| 4. | region.names |
names of regions |
The readVCF function expects a tabixed VCF file with a diploid GT field.
In case of haploid data, the GT field has to be transformed to a pseudo-diploid
field (such as 0 -> 0|0). An alternative is to use readData(..., format="VCF"),
which can read non-tabixed haploid and any kind of polyploid VCFs directly.
When approx=TRUE, the algorithm will apply a logical OR to the GT-field:
(0|0=0,1|0=1,0|1=1,1|1=1). Note, this is an approximation for diploid data, which will
speed up calculations. In case of haploid data, approx should be switched to TRUE.
If approx=FALSE, the full diploid information will be considered.
The ff-package PopGenome uses to store the SNP information limits total data size to
individuals * (number of SNPs) <= .Machine$integer.max
In case of very large data sets, the bigmemory package will be used;
this will slow down calculations (e.g. this package have to be installed first !!!).
Use the function vcf_handle <-.Call("VCF_open", filename)
to open a VCF-file and .Call("VCF_getSampleNames",vcf_handle)
to get and define the individuals which should be considered in the analysis.
See also readData(..., format="VCF") !
# NOT RUN {
# GENOME.class <- readVCF("...\chr1.vcf.gz", 1000, "1", 1, 100000)
# GENOME.class
# [email protected]
# GENOME.class <- neutrality.stats(GENOME.class,FAST=TRUE)
# show the result:
# get.sum.data(GENOME.class)
# [email protected]
# }
Run the code above in your browser using DataLab