gbs_filtering_thresholds.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### VCF filtering thresholds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#read in gatk table from VariantsToTable\n",
    "stats=read.delim(\"stats.table\",sep=\"\\t\", header=T)\n",
    "nsamples=376\n",
    "dim(stats)\n",
    "colnames(stats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for (i in 3:14)\n",
    "    {\n",
    "        plot(density(stats[,i], na.rm=T), main=colnames(stats)[i])\n",
    "        print(colnames(stats)[i])\n",
    "        print(sum(is.na(stats[,i])))\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#ExcessHet\n",
    "#ExcessHet>13 => pvalue<.05\n",
    "excesshet_thresh=13\n",
    "10^(-excesshet_thresh/10)\n",
    "stats_filt=subset(stats, ExcessHet<excesshet_thresh | is.na(ExcessHet))\n",
    "dim(stats_filt)\n",
    "dim(stats)[1]-dim(stats_filt)[1]\n",
    "plot(density(stats_filt$ExcessHet, na.rm=T), main=\"ExcessHet\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#inbreeding\n",
    "ib_thresh=.3\n",
    "stats_filt=subset(stats, InbreedingCoeff<ib_thresh & InbreedingCoeff>-ib_thresh | is.na(InbreedingCoeff))\n",
    "dim(stats_filt)\n",
    "dim(stats)[1]-dim(stats_filt)[1]\n",
    "plot(density(stats$InbreedingCoeff, na.rm=T), main=\"InbreedingCoeff\")\n",
    "abline(v=ib_thresh)\n",
    "abline(v=-ib_thresh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#MQ\n",
    "mq_thresh=40\n",
    "stats_filt=subset(stats, MQ>mq_thresh | is.na(MQ))\n",
    "dim(stats_filt)\n",
    "dim(stats)[1]-dim(stats_filt)[1]\n",
    "plot(density(stats$MQ, na.rm=T), main=\"MQ\")\n",
    "abline(v=mq_thresh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#MQRankSum\n",
    "mqrs_thresh=12.5\n",
    "stats_filt=subset(stats, MQRankSum<mqrs_thresh & MQRankSum>-mqrs_thresh | is.na(MQRankSum))\n",
    "dim(stats_filt)\n",
    "dim(stats)[1]-dim(stats_filt)[1]\n",
    "plot(density(stats$MQRankSum, na.rm=T), main=\"MQRankSum\")\n",
    "abline(v=mqrs_thresh)\n",
    "abline(v=-mqrs_thresh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#QD\n",
    "qd_thresh=2\n",
    "stats_filt=subset(stats, QD>qd_thresh | is.na(QD))\n",
    "dim(stats_filt)\n",
    "dim(stats)[1]-dim(stats_filt)[1]\n",
    "plot(density(stats$QD, na.rm=T), main=\"QD\")\n",
    "abline(v=qd_thresh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#NO.CALL proportion (NO.CALL=number of samples with missing genotypes)\n",
    "#0 don't allow missing data\n",
    "#1 keep all, regardless of how much missing\n",
    "nocall_thresh=.8\n",
    "stats_filt=subset(stats, NO.CALL<(nocall_thresh*nsamples))\n",
    "dim(stats_filt)\n",
    "dim(stats)[1]-dim(stats_filt)[1]\n",
    "plot(density(stats$NO.CALL/nsamples, na.rm=T), main=\"NO.CALL\")\n",
    "abline(v=nocall_thresh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#excesshet_thresh=\n",
    "#mq_thresh=\n",
    "#mqrs_thresh=\n",
    "#qd_thresh=\n",
    "#nocall_thresh=\n",
    "\n",
    "stats_filt=subset(stats, NO.CALL<(nocall_thresh*nsamples))\n",
    "stats_filt=subset(stats_filt, ExcessHet<excesshet_thresh | is.na(ExcessHet))\n",
    "stats_filt=subset(stats_filt, InbreedingCoeff<ib_thresh & InbreedingCoeff>-ib_thresh | is.na(InbreedingCoeff))\n",
    "stats_filt=subset(stats_filt, MQ>mq_thresh | is.na(MQ))\n",
    "stats_filt=subset(stats_filt, MQRankSum<mqrs_thresh & MQRankSum>-mqrs_thresh | is.na(MQRankSum))\n",
    "stats_filt=subset(stats_filt, QD>qd_thresh | is.na(QD))\n",
    "stats_filt=subset(stats_filt, QD>qd_thresh | is.na(QD))\n",
    "\n",
    "dim(stats_filt)[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for (i in 3:14)\n",
    "    {\n",
    "        plot(density(stats_filt[,i], na.rm=T), main=colnames(stats_filt)[i])\n",
    "    }"
   ]
  }
 ],
 "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.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}