\documentclass{beamer} \usetheme{Warsaw} \author{Sean Davis} \title{Variant Annotation \\ and Other Exercises} \section{Background} <>= library(knitr) opts_chunk$set(cache=TRUE,out.width='0.5\\textwidth',fig.align='center',warnings=FALSE,messages=FALSE) @ \begin{document} \begin{frame} \titlepage \end{frame} \begin{frame}[fragile]{A Example Variant} \begin{figure} \begin{center} \includegraphics[width=1\textwidth]{figure/BRAFV600E-SK-MEL-5} \caption{A true variant in a cancer sample as viewed in the Integrated Genomic Viewer (IGV) software} \end{center} \end{figure} \end{frame} \begin{frame}{Tasks} \begin{itemize} \item Read in a VCF file and locate all variants with respect to genes \item Predict the coding effect of the variants in a VCF file \item Investigate mutation spectrum in two cancer samples that have undergone exome sequencing \end{itemize} \end{frame} \section{Locate Variants with Respect to Genes} \begin{frame}[fragile]{Possible locations for variants} \begin{table}[h!] \begin{center} \begin{tabular}{l|l} \hline Location & Details \\ \hline coding & falls \emph{within} a coding region \\ fiveUTR & falls \emph{within} a 5' untranslated region \\ threeUTR & falls \emph{within} a 3' untranslated region \\ intron & falls \emph{within} an intron region \\ intergenic & does not fall \emph{within} a transcript associated with a gene \\ spliceSite & overlaps any portion of the first 2 or last 2 nucleotides of an intron \\ promoter & falls \emph{within} a promoter region of a transcript \\ \hline \end{tabular} \end{center} \caption{Variant locations} \label{table:location} \end{table} \end{frame} \begin{frame}[fragile]{Necessary Pieces of Information to Locate Variants} \begin{itemize} \item Gene models and transcripts \item Variants for annotation \end{itemize} \end{frame} \begin{frame}[fragile]{Actually Locating Variants} \scriptsize <>= library(VariantAnnotation) fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation") vcf <- readVcf(fl, genome="hg19") library(TxDb.Hsapiens.UCSC.hg19.knownGene) txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene codingvar <- locateVariants(vcf, txdb, CodingVariants()) seqlevels(vcf)=paste0('chr',seqlevels(vcf)) codingvar <- locateVariants(vcf, txdb, CodingVariants()) head(codingvar, 3) allvar <- locateVariants(vcf, txdb, AllVariants()) head(allvar) @ \end{frame} \section{Amino acid coding changes} \begin{frame}[fragile]{Predicting Coding Effects of Variants} \begin{itemize} \item Gene models and transcripts \item Variants for annotation \item \textit{Sequence information} \end{itemize} \end{frame} \begin{frame}[fragile]{Actually Predicting Variants} <>= library(BSgenome.Hsapiens.UCSC.hg19) coding <- predictCoding(vcf, txdb, seqSource=Hsapiens) coding[5:7] @ When the resulting \textit{varCodon} is not a multiple of 3 it cannot be translated. The consequence is considered a \textit{frameshift} and \textit{varAA} will be missing. <>= ## CONSEQUENCE is 'frameshift' where translation is not possible coding[mcols(coding)$CONSEQUENCE == "frameshift"] @ \end{frame} \begin{frame}[fragile]{Visualizing Variants in Genomic Context Using Gviz} \tiny <>= library(Gviz) txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene geneTrack = GeneRegionTrack(txdb,chromosome='chr22',start=50e6,end=51e6) variantTrack = AnnotationTrack(rowData(vcf)) plotTracks(list(geneTrack,variantTrack),from=50.5e6,to=50.6e6) @ \end{frame} \section{Investigating Some Cancer Samples} \begin{frame}{Transitions and Transversions} \begin{figure} \centering \includegraphics[width=0.6\textwidth]{TransitionTransversion} \label{fig:transversion1} \end{figure} \end{frame} \begin{frame}{Transitions and Transversions} \begin{figure} \centering \includegraphics[width=0.3\textwidth]{TransitionTransversion} \caption{DNA substitution mutations are of two types. Transitions are interchanges of two-ring purines (A--G) or of one-ring pyrimidines (C--T): they therefore involve bases of similar shape. Transversions are interchanges of purine for pyrimidine bases, which therefore involve exchange of one-ring and two-ring structures. Note that transversions (tv) are twice as likely as transitions(ti), leading to an expected ratio of ti/tv of 0.5.} \label{fig:transversion2} \end{figure} \end{frame} \begin{frame}{Spontaneous Mutation Alters ti/tv Ratio} \begin{figure} \centering \includegraphics[width=0.6\textwidth]{Transition_mutation} \label{fig:spontaneousmutation1} \end{figure} \end{frame} \begin{frame}{Spontaneous Mutation Alters ti/tv Ratio} \begin{figure} \centering \includegraphics[width=0.3\textwidth]{Transition_mutation} \caption{In the original double-stranded DNA molecule, A in the standard (amino) form pairs with T. During replication, the two strands separate. In the upper diagram, T pairs with A as usual, which replicates the wild-type sequence. In the lower diagram, A has undergone a tautomeric shift to the non-standard (imino) form A', \textit{which pairs with C}. In the next round of replication, the imino A' shifts back to the amino A form, which pairs with T, which again reproduces the wild-type sequence. Replication of the other strand pairs C with G. By comparison with the original molecule, the result is a T--C mutation. A tautomeric shift in one strand has produced a transition mutation in the complementary strand. This process leads to a ti/tv ratio that is typically much larger than the expected 0.5.} \label{fig:spontaneousmutation2} \end{figure} \end{frame} \begin{frame}{ti/tv ratio and false positive variants} \begin{itemize} \item If all observed variants are true positives, then we should observe approximately the ti/tv for the exome, which is reliably estimated at 3.3. \item If all observed variants are false positives, then we should observe the naive expected ti/tv, or 0.5. \item In reality, the observed ti/tv ratio is a mixture of true and false positives. \end{itemize} \end{frame} \begin{frame}{ti/tv ratio and false positive variants} Given the observed ti/tv ($titv_{obs}$), the false positive ti/tv ($titv_{fp}$), the true positive ti/tv ($titv_{tp}$), and $\alpha$, the proportion of SNVs that are true positives, we can write: \begin{equation} \alpha ( titv_{tp} ) + ( 1 - \alpha ) ( titv_{fp} ) = titv_{obs} \end{equation} Solving for 1 - $\alpha$ gives an estimate of the false positive rate. \begin{equation} \label{eq:fpr} FPR = 1 - \alpha = 1 - (titv_{obs} - titv_{fp})/(titv_{tp} - titv_{fp}) \end{equation} \end{frame} \begin{frame}[fragile]{What do we need to find ti/tv} \begin{centering} \begin{verbatim} alt ref A C G T A 0 520 3195 328 C 592 0 770 3598 G 3782 756 0 568 T 359 3349 542 0 \end{verbatim} \end{centering} \end{frame} \begin{frame}{Investigating Cancer Sample Tasks} \begin{itemize} \item Load two VCF files representing a cancer sample \item Define a function to produce a data frame of SNPs with two columns: \begin{itemize} \item Reference allele \item Alternate allele (the variant) \end{itemize} \item Define a function to calculate the ti/tv ratio for a vcf file \item Subset our two VCF files to include only the ``Type 2'' variants \item Calculate the ti/tv ratio for each sample \item Make a plot of the mutation spectrum for each sample type \end{itemize} \end{frame} \begin{frame}[fragile]{Examining our VCF file} \scriptsize <>= # you may need to change the file name to load in the data vcfMel = readVcf('../data/SK-MEL-5.vcf.gz',genome='hg19') vcfLung = readVcf('../data/A549_ATCC.vcf.gz',genome='hg19') nrow(vcfMel) nrow(vcfLung) vcfMel vcfLung rowData(vcfMel) info(vcfMel) rowData(vcfMel) header(vcfMel) ref(vcfMel) alt(vcfMel) @ \end{frame} \begin{frame}[fragile]{Our functions} \pause \scriptsize <>= vcf2snpDF = function(vcf) { refall = as.character(ref(vcf)) altall = as.character(unlist(alt(vcf))[start(PartitioningByEnd(alt(vcf)))]) tmpDF = data.frame(ref=refall,alt=altall,stringsAsFactors=FALSE) # snps only tmpDF = tmpDF[nchar(tmpDF$ref)==1 & nchar(tmpDF$alt)==1,] return(tmpDF) } @ \pause <>= titv = function(vcf) { variantDF = vcf2snpDF(vcf) tbl = table(variantDF) ti = tbl['C','T']+tbl['T','C']+tbl['A','G']+tbl['G','A'] tv = tbl['A','C']+tbl['C','A']+tbl['A','T']+tbl['T','A']+tbl['C','G']+tbl['G','C']+tbl['G','T']+tbl['T','G'] return(list(ti=ti,tv=tv,tbl=tbl,titv=ti/tv)) } @ \end{frame} \begin{frame}[fragile]{Subset to include only ``Type 2'' variants} \scriptsize <>= vcfMelT2 = vcfMel[grep('Type2',rowData(vcfMel)$FILTER)] vcfLungT2 = vcfLung[grep('Type2',rowData(vcfLung)$FILTER)] @ <>= titv(vcfMelT2) titv(vcfLungT2) @ \end{frame} \begin{frame}[fragile]{And a plot of the mutation spectrum} \scriptsize <>= library(ggplot2) ggplot(vcf2snpDF(vcfMelT2),aes(x=ref,fill=alt))+geom_bar() ggplot(vcf2snpDF(vcfLungT2),aes(x=ref,fill=alt))+geom_bar() @ \end{frame} \begin{frame}{Remnants of DNA Damage in Cancer Cells} In lung cancer, polycyclic aromatic hydrocarbons (PAH) cause an increase in transversions (G-T and C-A), leading to a decrease in ti/tv. In melanoma, thymine dimers associated with UV damage lead to increased transitions, leading to an increase in ti/tv. \end{frame} \end{document}