Working with VCF files and Trees

About this notebook

This is R Markdown Notebook continue my studies of graphing trees using ggtree. I have been working with Variant Call Formats (VCFs, See SAM/BAM and related specifications). The data used is from ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/pilot_data/release/2010_07/exon/snps/ file CEU.exon.2010_03.genotypes.vcf.gz which is part of the 1000genomes project.

Samtools Bcftools or GATK HaplotypeCaller can be used for variant calling. These programs will generate vcf files that will need to be filtered for SNPs and other paramaters. For GATK see From FastQ Data to High-Confidence Variant Calls: The Genome Analysis Toolkit Best Practices Pipeline (https://doi.org/10.1002/0471250953.bi1110s43) by Van der Auwera et. al.

This document was written in RStudio using R Markdown and thus it may contains typos. I take no responsibility for any damage you may do based upon knowledge gained from this document. Please read all instructions before and when installing any computer program(s), application(s) and module(s) / package(s). Please read all instructions before using any program(s), application(s) and module(s) / package(s).

## Loading the following R libraries
library(gdsfmt)
library(SNPRelate) 
library(ggplot2)
library(ggtree)
library(ape)

SNPRelate Code

### Save the path to the GATK vcfs to  variables

exon2010_03genotypes <- "CEU.exon.2010_03.genotypes.vcf"

### turn the VCF file into a less data intensive form (GDS) for easier computing

snpgdsVCF2GDS(exon2010_03genotypes,"CEUexon2010_03genotype.gds",method ="biallelic.only")

### Preparing the data so it is formatted correctly to create a dissimilarity matrix.

genofileExon2010_03 <- snpgdsOpen("CEUexon2010_03genotype.gds")

set.seed(100) ## making the code reproducible

ibs_Exon2010_03 <- snpgdsHCluster(snpgdsIBS(genofileExon2010_03,num.thread=2, autosome.only=FALSE))

Identity-By-State (IBS) analysis on genotypes:
Excluding 49 SNPs (monomorphic: TRUE, MAF: NaN, missing rate: NaN)
Working space: 90 samples, 3,440 SNPs
    using 2 (CPU) cores
IBS:    the sum of all selected genotypes (0,1,2) = 488121
Sun Dec 15 16:24:52 2019    (internal increment: 65536)

[..................................................]  0%, ETC: ---        
[==================================================] 100%, completed, 0s
Sun Dec 15 16:24:52 2019    Done.

rvExon2010_03 <- snpgdsCutTree(ibs_Exon2010_03)

Determine groups by permutation (Z threshold: 15, outlier threshold: 5):
Create 1 groups.

### Saving the dendrograms to new Variables

treeExon2010_03 = rvExon2010_03$dendrogram

plot(rvExon2010_03$dendrogram,horiz=T, main ="CEU.exon.2010_03.genotypes.vcf SNP Tree" )

ggtree and ape Code

tree2 = ggtree(as.phylo(as.hclust(treeExon2010_03)), layout="circular",color='darkgreen', branch.length="branch.length")+ geom_tiplab(size=2.5, aes(angle=angle))+ ggtitle("CEU.exon.2010_03.genotypes.vcf SNP Tree")
tree2

### converting dendrograms to class hclust and to new variables

hcExon2010_03 = as.hclust(rvExon2010_03$dendrogram)

### Making the hclust object into a phylo object in ape

thisExon2010_03 <- as.phylo(hcExon2010_03) 

#### Writing to a newick tree file
write.tree(phy=thisExon2010_03 , file="CEUexon2010_03genotypes.newick") 

Newick Format of trees

The newick format of the saved file

(NA11930:0.06076188698,(NA12043:0.06033848847,(NA11994:0.05995671748,(NA07051:0.05928583721,(NA12776:0.05909499358,(NA12287:0.05893776241,(NA12546:0.05884010489,(NA12413:0.05845146464,(NA12154:0.05827547638,(NA11920:0.0581097818,(NA12045:0.05805355315,(NA06986:0.05780013849,(NA07347:0.05759102976,(NA11918:0.05737944516,(NA11992:0.05719160744,(NA12003:0.05692680589,(NA11831:0.05676047966,(NA12383:0.0566187369,(NA12829:0.05651646832,(NA07346:0.05637855456,(NA12234:0.05632008556,(NA11832:0.05614709088,((NA11830:0.05534591195,NA12144:0.05534591195):0.0005054755299,(NA12005:0.05573779242,((NA07000:0.05450844391,NA12828:0.05450844391):0.001120513186,(NA12762:0.05545329394,(NA12775:0.05532717655,((NA07037:0.05365178262,NA12750:0.05365178262):0.001607203476,(NA12814:0.05491784382,((NA07357:0.05336911151,NA12155:0.05336911151):0.0009418059885,(NA12340:0.05427117222,((NA12273:0.05218100847,(NA12891:0.04903433595,(NA12878:0.03687796209,NA12892:0.03687796209):0.01215637387):0.003146672515):0.001705880089,((NA12347:0.05210896309,NA12400:0.05210896309):0.001425151892,(NA12843:0.05336717085,(NA12890:0.05311198727,((NA10847:0.05055305911,NA12760:0.05055305911):0.00249025893,(NA12749:0.05289726265,(NA11829:0.05270030097,(NA11881:0.05259554217,(NA12283:0.05250511042,((NA12044:0.04792248603,NA12286:0.04792248603):0.004377337486,((NA11843:0.0506318306,NA12748:0.0506318306):0.001373398837,((NA12249:0.04870010985,NA12414:0.04870010985):0.003204675922,(NA12751:0.05170800322,(NA12716:0.05151973142,((NA11893:0.04950325454,NA12718:0.04950325454):0.001522371039,(NA12275:0.05080338904,(NA12348:0.0505532199,((NA12342:0.04755163225,NA12830:0.04755163225):0.00264395034,(NA12058:0.05005136838,((NA06985:0.04616372981,NA12006:0.04616372981):0.00364515562,(NA07048:0.04960065804,((NA12717:0.04779411765,NA12874:0.04779411765):0.001326106937,(NA06984:0.04858367492,(NA11995:0.04829443111,(NA11840:0.04805128445,((NA12272:0.04611194591,NA12815:0.04611194591):0.001379913366,((NA12812:0.04613970588,NA12872:0.04613970588):0.0007757612955,(NA12763:0.04662911816,(NA12341:0.04584688474,((NA10851:0.04375951294,NA12761:0.04375951294):0.001126877505,(NA12004:0.04473378893,((NA12489:0.0414746894,(NA12873:0.04057040186,(NA06989:0.03636363636,NA12156:0.03636363636):0.004206765499):0.0009042875394):0.002076120524,(NA06994:0.04250422171,(NA12282:0.03914448173,(NA11919:0.0382151604,(NA12842:0.03152866242,NA12889:0.03152866242):0.006686497981):0.0009293213316):0.003359739972):0.00104658822):0.001182979005):0.0001526015118):0.000960494293):0.0007822334215):0.0002863490207):0.0005763920931):0.0005594251839):0.0002431466516):0.0002892438108):0.0005365496666):0.0004804334596):0.0002082273857):0.0002424829515):0.0001442142043):0.0003576373128):0.0002501691423):0.000222236538):0.0004941058371):0.000188271806):0.0001967825513):0.0001004436655):0.0002945940816):0.000205286903):9.043174815e-05):0.0001047587986):0.0001969616809):0.0001460553881):6.866923629e-05):0.0002551835736):0.0001669441369):0.0003527735697):0.0003842836661):3.974527651e-05):0.0006069263235):0.0003411422794):6.819045188e-05):0.0001261173869):0.0001756631553):0.0001088353249):0.0001135950605):0.0002957034002):0.0001729946798):5.846899577e-05):0.0001379137604):0.0001022685835):0.0001417427572):0.0001663262341):0.0002648015536):0.0001878377111):0.0002115846046):0.0002091087326):0.0002534146622):5.622864423e-05):0.0001656945766):0.0001759882635):0.0003886402467):9.765752529e-05):0.000157231165):0.0001908436331):0.000670880274):0.0003817709915):0.000423398503);

The newick file format can be used with many programs / software that are capable of generating trees like the following:

• ete3, ETE Toolkit (python module for phylogenetic analysis) : ETE 3: Reconstruction, analysis and visualization of phylogenomic data. Jaime Huerta-Cepas, Francois Serra and Peer Bork. Mol Biol Evol 2016; https://doi.org/10.1093/molbev/msw046
• FigTree

References

Tutorials used:

• Tutorials for the R/Bioconductor Package SNPRelate by Xiuwen Zheng

• Data Integration, Manipulation and Visualization of Phylogenetic Trees by Guangchuang Yu

R/Bioconductor Libraries for data analysis and visualization:

• R
• RStudio
• ggtree
• ape
• SNPRelate
• RStudio
• R Markdown