module-8-variantcalling/Variant_calling_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7fb437a4",
   "metadata": {},
   "source": [
    "### Variant calling module pt 2\n",
    "\n",
    "**CMM262, Winter 2022**\n",
    "\n",
    "Kyle Gaulton, kgaulton@gmail.com\n",
    "<br>\n",
    "<br>\n",
    "Variant genotype and annotation filtering and visualization in R\n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "637ec0fb",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# load required packages\n",
    "\n",
    "library(tidyverse)\n",
    "library(reshape2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd2e9426",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load tab-delimited genotype file\n",
    "\n",
    "geno <- read.table(\"GMTS.geno.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8cc2e3e",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "head(geno)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5c1c367",
   "metadata": {},
   "outputs": [],
   "source": [
    "colnames(geno) <- c(\"chr\",\"raw_pos\",\"id\",\"raw_ref\",\"raw_alt\",\"qual\",\"filter\",\"AC\",\"AF\",\"AN\",\"BaseQRankSum\",\"ClippingRankSum\",\"DP\",\"DS\",\"END\",\"ExcessHet\",\"FS\",\"Hrun\",\"HaplotypeScore\",\"InbreedingCoeff\",\"MLEAC\",\"MLEAF\",\"MQ\",\"MQRankSum\",\"QD\",\"RAW_MQ\",\"ReadPosRankSum\",\"SOR\",\"SAMPLE\",\"AD\",\"DP.1\",\"GQ\",\"GT\",\"MIN_DP\",\"PGT\",\"PID\",\"PL\",\"RGQ\",\"SB\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f04ca72b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "head(geno)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77c9536e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create a unique variant ID\n",
    "geno <- geno %>% mutate(varID=paste(chr,raw_pos,raw_ref,raw_alt,sep=\":\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5f011b0",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "head(geno)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecc851a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#add SNP/indel variant type info\n",
    "\n",
    "geno<-geno %>% mutate(type=ifelse(nchar(as.character(raw_ref))>1 | nchar(as.character(raw_alt))>1,\"indel\",\"SNV\"))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eed13dfb",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "head(geno)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51998944",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set sample column to character\n",
    "\n",
    "geno$SAMPLE <- as.character(geno$SAMPLE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19e769ea",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "head(geno$SAMPLE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86a2a0b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of variants\n",
    "geno %>% select(varID) %>% unique() %>% nrow()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9ac9ee6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of samples\n",
    "geno %>% select(SAMPLE) %>% unique() %>% nrow()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fd3be49",
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of each genotype\n",
    "table(geno$GT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78ec9f9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# retain only het and hom alt genotypes\n",
    "\n",
    "geno <- geno %>% filter(GT==\"0/1\" | GT==\"1/1\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bf89577",
   "metadata": {},
   "outputs": [],
   "source": [
    "# retain just snp variants\n",
    "geno$raw_ref <- as.character(geno$raw_ref)\n",
    "geno$raw_alt <- as.character(geno$raw_alt)\n",
    "geno$AD <- as.character(geno$AD)\n",
    "geno2<-geno %>% filter(!str_detect(AD,\"\\\\.\"))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "226c5f5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#getting allelic fraction\n",
    "geno2<-geno2 %>% separate(AD,into=c(\"Rcnt\",\"Acnt\"),remove = F, convert=T, sep=\",\") %>% mutate(AF=Acnt/(Rcnt+Acnt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4cde03e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#number of variant of each GT and type per sample\n",
    "geno %>% select(SAMPLE,GT,type,varID) %>% \n",
    "  unique() %>% \n",
    "  group_by(SAMPLE,GT,type) %>% \n",
    "  tally() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3524234",
   "metadata": {},
   "outputs": [],
   "source": [
    "# average number of each variant type per samples\n",
    "\n",
    "geno %>% select(SAMPLE,type,varID) %>% \n",
    "  unique() %>% \n",
    "  group_by(SAMPLE,type) %>% \n",
    "  summarize(Nvariants=length(varID)) %>%\n",
    "  group_by(type) %>%\n",
    "  summarize(medVar=median(Nvariants),minVar=min(Nvariants),maxVar=max(Nvariants))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2060175",
   "metadata": {},
   "outputs": [],
   "source": [
    "geno3 <- geno %>% select(SAMPLE,type,varID) %>% \n",
    "  unique() %>% \n",
    "  group_by(SAMPLE,type) %>% \n",
    "  summarize(Nvariants=length(varID))\n",
    "\n",
    "ggplot(geno3,aes(x=Nvariants,y=SAMPLE,col=type)) + geom_point() + theme_classic()\n",
    "ggplot(filter(geno3,type==\"SNV\"),aes(x=Nvariants)) + geom_histogram() + theme_classic()\n",
    " ggplot(filter(geno3,type==\"indel\"),aes(x=Nvariants)) + geom_histogram() + theme_classic()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e00151d7",
   "metadata": {},
   "source": [
    "**Filtering and visualizing variant annotations**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "914961b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "var <- read.delim2(\"myanno.21.txt\",header=F,sep=\"\\t\")\n",
    "colnames(var)<-c(\"chr\",\"pos\",\"raw_pos\",\"ref\",\"raw_ref\",\"alt\",\"raw_alt\",\"mut_type\",\"function\",\"region_type\",\"region_name\",\"refGene.name\",\"refGene.name2\",\"dbSNP.name\",\"ExAC_r0_3_sites_vep.name\",\"ExAC_r0_3_sites_vep.AF\",\"dbNSFP.SIFT_pred\",\"dbNSFP.Polyphen2_HDIV_pred\",\"dbNSFP.MutationTaster_pred\",\"dbNSFP.MetaSVM_pred\",\"dbNSFP.clinvar_clnsig\",\"clinvar.CLNSIG\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83c02bd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "head(var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37720362",
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove leading and trailing white spaces\n",
    "\n",
    "var$chr <- paste0(\"chr\",trimws(var$chr))\n",
    "var$raw_pos<-as.numeric(trimws(var$raw_pos))\n",
    "var$raw_ref<-trimws(var$raw_ref)\n",
    "var$raw_alt<-trimws(var$raw_alt)\n",
    "var$refGene.name2<-trimws(var$refGene.name2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f96ffcdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create a unique variant ID\n",
    "\n",
    "var <- var %>% mutate(varID=paste(chr,raw_pos,raw_ref,raw_alt,sep=\":\"))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f5b82e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# summarize number of variants\n",
    "\n",
    "var %>% select(varID) %>% unique() %>% nrow()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea452acb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# summarize number of variants by functional type\n",
    "\n",
    "var %>% group_by(mut_type) %>% tally()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "304ab99e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# summarize number of variants by functional type for a specific gene\n",
    "\n",
    "var %>% filter(refGene.name2==\"DOPEY2\") %>% group_by(mut_type) %>% tally()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fcc2e3cb",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# summarize and report variants with potential clinical significance\n",
    "\n",
    "var %>% filter(dbNSFP.clinvar_clnsig==5) %>% select(dbSNP.name,mut_type,refGene.name2) %>% unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7f55d4c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# nonsyn variants per gene\n",
    "var2 <- var %>% filter(grepl(\"nonsynonymous\",mut_type))\n",
    "\n",
    "ggplot(var2,aes(refGene.name2))+geom_bar()+coord_flip()+theme(axis.text.y = element_text(size = 5)) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87dc4744",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "topGenes <- var2 %>% group_by(refGene.name2) %>% tally() %>% top_n(20,n) %>% select(refGene.name2)\n",
    "topGenes\n",
    "var_top<-left_join(topGenes,var2)\n",
    "ggplot(var_top,aes(refGene.name2))+geom_bar() +coord_flip()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "381e054a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# join genotype and annotation data\n",
    "\n",
    "geno_lite <- geno %>% select(varID,SAMPLE,type,qual,GT,DP)\n",
    "anno_lite <- var %>% select(varID,refGene.name2,mut_type,region_type,`function`)\n",
    "data <- left_join(geno_lite,anno_lite)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bad61a35",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# number of nonsynonymous variants per sample\n",
    "\n",
    "data %>% group_by(SAMPLE) %>% filter(grepl(\"frameshift\",mut_type)) %>% summarize(NS=length(varID))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a84a4a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# comparing quality by type and GT\n",
    "\n",
    "ggplot(data,aes(type,qual))+geom_boxplot()+facet_wrap(~SAMPLE)\n",
    "\n",
    "ggplot(data,aes(type,qual))+geom_boxplot(outlier.shape = NA)+geom_jitter(width=0.1)+facet_wrap(~SAMPLE)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b13c034e",
   "metadata": {},
   "outputs": [],
   "source": [
    "sdata <- data %>% filter(SAMPLE==\"201849403\") %>% filter(type==\"SNV\")\n",
    "\n",
    "ggplot(sdata,aes(type,qual))+geom_boxplot(outlier.shape = NA)+geom_jitter(width=0.1)+facet_grid(~mut_type)\n",
    "\n",
    "sdata2 <- data %>% filter(SAMPLE==\"201849403\") %>% filter(type==\"indel\")\n",
    "\n",
    "ggplot(sdata2,aes(type,qual))+geom_boxplot(outlier.shape = NA)+geom_jitter(width=0.1)+facet_grid(~mut_type)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be4dbdd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# coverage depth distribution\n",
    "\n",
    "ggplot(data,aes(type,log10(DP)))+geom_boxplot()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c46c743e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot variants per sample by functional category\n",
    "\n",
    "ggplot(data,aes(SAMPLE,fill=mut_type))+geom_bar()\n",
    "\n",
    "del_vars <- data %>% filter(grepl(\"^frameshift|^stop\",mut_type))\n",
    "\n",
    "ggplot(del_vars,aes(SAMPLE,fill=mut_type))+geom_bar()\n",
    "\n",
    "# plot "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c73d0c52",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5d4de09",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}