The VariantAnnotation package allows .vcf files to be imported. The readVcf function can be used and requires the file of a vcf file. You also need to specify a genome name. As we have seen already with other Bioconductor objects, typing the name of the object will print a summary to the screen. In the case of a “CollapsedVCF” object this is very detailed
We will use the file combined.chr20.subset.freebayes.vcf that you should have generated in the previous section.
The “header” information can be extracted using header. This will contain the definitions of all the per-genotype (INFO) and per-sample information stored in the file
INFOAs described above, the INFO column in a .vcf file gives per-variant metadata as a series of key-value pairs. We can interrogate these data using the info function, which returns a data-frame-like object. Consequently, we can use the $ operator to select a particular column of interest.
The name of each row is derived from the genomic location and base-change for each variant. There are too many columns to go through in detail, so we will just describe a few that are likely to be common to different callers.
infoMatrix <- info(hapmap.calls)
infoMatrixThe called genotypes can be accessed using the geno function. Recall that in the .vcf file, we have one column of genotype information for each sample, with each column consisting of key:value pairs. Using geno we can access all the values for a particular keys and be able to compare across samples.
The column names of the data frame returned by geno are the same as the FORMAT in the .vcf description. Thus we can use a $ operator to access a particular set of values; GT in this case
geno(hapmap.calls)
head(geno(hapmap.calls)$GT)The output allows us to compare the genotype for a particular variant across all samples; which we could not easily do from the .vcf.
Usually each entry is 0/0 for a homozygous reference, 0/1 for a heterozygous call and 1/1 for a homozyous alternate allele. An entry of . indicates a position where no call could be made due to insufficient data. Moreover, we can also find 0/2 and 1/2 in rare cases where a second alternative allele was found.
With table we can tabulate the calls between one sample and another
table(geno(hapmap.calls)$GT[,1])
table(geno(hapmap.calls)$GT[,2])
table(geno(hapmap.calls)$GT[,1], geno(hapmap.calls)$GT[,2])
In this section we will see how we can overlap the calls we have made with other genomic features. For example, we are often interested in how many calls were made within a particular gene of interest.
For this section, we are going to use all calls made on chromosome 20 for a single sample; NA12878.
The rowRanges function will retrieve the positions of our variants as a familiar GRanges object. Along with the usual positional information, we also get extra “metadata” (mcols) about the base-change, a quality-score from freebayes and a placeholder for a filter (which has not been applied in this case).
NA12878.calls <- readVcf("NA12878.chr20.freebayes.vcf","hg19")
NA12878.calls <- readVcf("NA12878.chr20.freebayes.vcf","hg19")
NA12878.calls.ranges <- rowRanges(NA12878.calls)
NA12878.calls.ranges
GRanges object with 74812 ranges and 5 metadata columns:
seqnames ranges strand |
<Rle> <IRanges> <Rle> |
20:61795_G/T 20 [61795, 61795] * |
20:63244_A/C 20 [63244, 63244] * |
20:65900_G/A 20 [65900, 65900] * |
20:66370_G/A 20 [66370, 66370] * |
20:67500_TT/TTGGTATCTAGT 20 [67500, 67501] * |
... ... ... ... .
20:62961014_CTTTTTTTTTTTTTA/CTTTTTTTTTTTTTTA 20 [62961014, 62961028] * |
20:62961318_G/A 20 [62961318, 62961318] * |
20:62961724_C/T 20 [62961724, 62961724] * |
20:62962130_C/T 20 [62962130, 62962130] * |
20:62962891_C/T 20 [62962891, 62962891] * |
paramRangeID REF
<factor> <DNAStringSet>
20:61795_G/T <NA> G
20:63244_A/C <NA> A
20:65900_G/A <NA> G
20:66370_G/A <NA> G
20:67500_TT/TTGGTATCTAGT <NA> TT
... ... ...
20:62961014_CTTTTTTTTTTTTTA/CTTTTTTTTTTTTTTA <NA> CTTTTTTTTTTTTTA
20:62961318_G/A <NA> G
20:62961724_C/T <NA> C
20:62962130_C/T <NA> C
20:62962891_C/T <NA> C
ALT QUAL
<DNAStringSetList> <numeric>
20:61795_G/T T 95.0616
20:63244_A/C C 88.3208
20:65900_G/A A 44.2954
20:66370_G/A A 57.5310
20:67500_TT/TTGGTATCTAGT TTGGTATCTAGT 12.1582
... ... ...
20:62961014_CTTTTTTTTTTTTTA/CTTTTTTTTTTTTTTA CTTTTTTTTTTTTTTA 7.2667200
20:62961318_G/A A 46.1486000
20:62961724_C/T T 0.0518025
20:62962130_C/T T 57.5250000
20:62962891_C/T T 18.2903000
FILTER
<character>
20:61795_G/T .
20:63244_A/C .
20:65900_G/A .
20:66370_G/A .
20:67500_TT/TTGGTATCTAGT .
... ...
20:62961014_CTTTTTTTTTTTTTA/CTTTTTTTTTTTTTTA .
20:62961318_G/A .
20:62961724_C/T .
20:62962130_C/T .
20:62962891_C/T .
-------
seqinfo: 86 sequences from hg19 genomeIf we happen to know the genomic region corresponding to a particular gene, we can restrict our list of variants using standard R syntax.
chr20:4,700,556-4,711,106 on the human genome version hg19NA12878.calls.ranges[start(NA12878.calls.ranges) > 4700556 & end(NA12878.calls.ranges) < 4711106]However, we might want to something more sophisticated and only consider variants in coding regions. For this we can take advantage of some pre-built packages in Bioconductor.
Aside from the many useful software packages, Bioconductor also provides numerous annotation resources that we can utilise in our analysis. Firstly, we have a set of organism-level packages that can translate between different types of identifer. The package for humans is called org.Hs.eg.db. The advantage of such a package, rather than services such as biomaRt, is that we can do queries offline. The packages are updated every 6 months, so we can always be sure of what version of the relevant databases are being used.
library(org.Hs.eg.db)
org.Hs.eg.dbThere are several types of “key” we can use to make a query, and we have to specify one of these names.
keytypes(org.Hs.eg.db)
For the given keytype we have chosen, we can also choose what data we want to retrieve. We can think of these as columns in a table, and the pre-defined values are given by:-
columns(org.Hs.eg.db)
eg <- select(org.Hs.eg.db, keys=c("BRCA1","PTEN"), keytype = "SYMBOL",columns = c("REFSEQ","ENSEMBL"))You should see that the above command prints a message to the screen:- select() returned 1:many mapping between keys and columns. This is not an error message and R has still been able to generate the output requested.
egIn this case, we have “many”" (well, two) values of ENSEMBL for the gene PTEN. In practice this means we probably want to think carefully about merging this data with other tables.
org.Hs.eg.db package to retrieve the Entrez Gene ID for the Gene PRNDYou might expect to be able to retrieve information about the coordinates for a particular gene using the same interface. This was supported until recently, but the recommended approach now is to use another class of packages which describe the structure of genes in more detail.
The packages with the prefix TxDb.... represent the structure of all genes for a given organism in an efficient manner. For humans, we can use the package TxDb.Hsapiens.UCSC.hg19.knownGene to tell us about transcripts for the hg19 build of the genome. The package was generated using tables from the UCSC genome browser
As with the org.Hs.eg.db package we can load the package and inspect the kind of mappings available to us.
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
columns(txdb)keytypes(txdb)You’ll see that all the mappings are regarding the coordinates and IDs of various genomic features. There is only one type of identifier used, in this case Entrez ID. If we know a genes Entrez ID, we can get the exon coordinates with the following query.
org.Hs.eg.db package to retrieve the exon coordinates for the Gene PRND
mygene <- It is useful to be able to retrive the coordinates in this manner. However, we should now be familiar with the way intervals can be represented using GRanges. We have the ability to create a GRanges object from the result:-
my.gr <- GRanges(mygene$EXONCHROM, IRanges(mygene$EXONSTART,mygene$EXONEND))
my.grThe txdb packages also allow us to construct GenomicRanges representations for the exon structure of all genes. The results is a list, which the names of the list being a gene ID in Entrez format
all.exons <- exonsBy(txdb, "gene")
all.exons
my.gene <- all.exons[["23627"]]
my.geneWe are almost ready to do the overlap. However, there is an inconsistency in the naming conventions of the two sets of regions we are trying to overlap;
seqlevelsStyle(NA12878.calls.ranges)
seqlevelsStyle(my.gene)Fortunately, we can rename the chromosome names of the variants with a single call. We can also modify the object so only information about chromosome 20 is retained (later-on we would receive an error due to the MT chromosomes being different length)
seqlevelsStyle(NA12878.calls.ranges) <- "UCSC"
NA12878.calls.ranges <- keepSeqlevels(NA12878.calls.ranges, "chr20")The %over% function can be comapre two sets of ranges and produce a logical vector. Each entry in this vector being whether a particular location in the first set of ranges is present in the other. The vector can then be used to subset the variants in the usual manner.
in.gene <- NA12878.calls.ranges %over% my.gene
NA12878.calls.ranges[in.gene]We can even write-out a vcf file containing just the positions we have identified
writeVcf(NA12878.calls[in.gene],filename = "selected.variants.vcf")If all we cared about was the number of variants, we could use countOverlaps
countOverlaps(my.gene, NA12878.calls.ranges)regions.of.interest.bed.bed file using the import function in the rtracklayer packageregions.of.interest.bed contains all exon coordinates of genes on chromosome 20. How could you generate such a file?
unlist(all.exons) will give a GRanges object with one entry per-exon (not in the list structure)export in rtracklayer will write out various file types from GRanges objectsexonsBy to see what other options are available for extracting genomic featuresfreebayes to call SNVs from a set of healthy individuals
.vcf format contains a rich description of the called variants.vcf filesGenomicRanges to overlap our calls with other genomic features of interest.vcf would probably involve other non-R tools
Commands used to generate bam files for genotype calling
samtools view -h ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/other_exome_alignments/NA19239/exome_alignment/NA19239.mapped.solid.mosaik.YRI.exome.20111114.bam 22 | samtools view -bS - > /data/hapmap/NA19239.chr22.bam
samtools index /data/hapmap/NA19239.chr22.bam
wget ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/phase3/data/NA12878/alignment/NA12878.chrom20.ILLUMINA.bwa.CEU.low_coverage.20121211.bam -O NA12878.chr20.bam
samtools index NA12878.chr20.bam
wget ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/phase3/data/NA12874/alignment/NA12874.chrom20.ILLUMINA.bwa.CEU.low_coverage.20130415.bam -O NA12874.chr20.bam
samtools index NA12874.chr20.bam