gbs_prelimPCA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### initial analysis for species confirmation and identification of bad samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#input file\n",
    "vcffile=\"prelim.recode.vcf\"\n",
    "\n",
    "#\"databases\"\n",
    "samplefile=\"/home/megan/megan/research/eucalyptus/eucalyptus_data/Emelliodora_PlantsSamples.csv\"\n",
    "accessionfile=\"/home/megan/megan/research/eucalyptus/eucalyptus_data/Emelliodora_Accessions.csv\"\n",
    "\n",
    "#count header lines in command line (grep -c \"##\" batch_1.vcf)\n",
    "skipr=8\n",
    "#non sample columns\n",
    "skipc=9\n",
    "\n",
    "thresh.missing.loci=90  #higher number allows more missing data/keeps more loci\n",
    "thresh.missing.samples=90  #higher number allows more missing data/keeps more samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "library(oz)\n",
    "options(jupyter.plot_mimetypes = \"image/png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#read in vcf, skupping ##header rows\n",
    "vcf=read.delim(vcffile, sep=\"\\t\", header=T, skip=skipr)\n",
    "#remove non sample columns\n",
    "vcf=vcf[c(-skipc:-1)]\n",
    "#get just genotypes\n",
    "genos=apply(vcf,2,substr,1,3)\n",
    "#transpose \n",
    "genos=t(genos)\n",
    "#clean sample names\n",
    "row.names(genos)=substr(row.names(genos),9,14)\n",
    "#genos[1:5,1:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#add sample metadata\n",
    "sampleinfo=read.csv(samplefile, header=T)  #read in sample information\n",
    "samplematches=match(sampleinfo$SampleID, row.names(genos))  #subset for just samples in this analysis\n",
    "samples=sampleinfo[!is.na(samplematches),][order(na.omit(samplematches)),]\n",
    "\n",
    "#incorporate latitude and longitude from accession\n",
    "accessioninfo=read.csv(accessionfile, header=T)  #read in accession information\n",
    "accessionmatches=match(samples$AccessionID, accessioninfo$AccessionID)\n",
    "samples=cbind(samples,accessioninfo$Latitude[accessionmatches],accessioninfo$Longitude[accessionmatches],\n",
    "              accessioninfo$PopulationName[accessionmatches])\n",
    "names(samples)[names(samples)==\"accessioninfo$Latitude[accessionmatches]\"]=\"Latitude\"\n",
    "names(samples)[names(samples)==\"accessioninfo$Longitude[accessionmatches]\"]=\"Longitude\"\n",
    "names(samples)[names(samples)==\"accessioninfo$PopulationName[accessionmatches]\"]=\"PopulationName\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#map samples\n",
    "oz(ylim=c(min(samples$Latitude)-1,max(samples$Latitude)+1),\n",
    "   xlim=c(min(samples$Longitude)-1, max(samples$Longitude)+1))\n",
    "text(samples$Longitude,samples$Latitude, \n",
    "     label=samples$PopulationName, cex = .6) #add population names\n",
    "abline(v=130) #use to define boundary for geographic outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#identify geographic outliers\n",
    "minlong=130  #set minimimum longitude for filtering out geographic outliers\n",
    "geog.out=as.character(samples[which(samples$Longitude < minlong),]$SampleID)\n",
    "length(geog.out)\n",
    "geog.out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dim(genos)\n",
    "\n",
    "#identify missing data\n",
    "missing.filt=genos==\"./.\"                 \n",
    "\n",
    "#identify bad loci\n",
    "missPerLocus.filt=colSums(missing.filt)   #count missing data per locus\n",
    "locus.filt=which(missPerLocus.filt < thresh.missing.loci / 100 * nrow(genos))\n",
    "\n",
    "#identify bad samples\n",
    "missPerSample.filt=rowSums(missing.filt)\n",
    "ind.filt=which(missPerSample.filt < thresh.missing.samples / 100 * ncol(genos))\n",
    "\n",
    "#filter\n",
    "genos.filt=genos[ind.filt,locus.filt]\n",
    "dim(genos.filt)\n",
    "failed=names(which(missPerSample.filt > thresh.missing.samples / 100 * ncol(genos)))\n",
    "length(failed)\n",
    "failed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#put in matrix format\n",
    "genos.filt[genos.filt == \"./.\"] <- \"NA\"\n",
    "genos.filt[genos.filt == \"1/1\"] <- \"0\"\n",
    "genos.filt[genos.filt == \"0/1\"] <- \"1\"\n",
    "genos.filt[genos.filt == \"1/0\"] <- \"1\"\n",
    "genos.filt[genos.filt == \"0/0\"] <- \"2\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#calculate PCA\n",
    "geno.dist=dist(genos.filt)  #generate distance matrix\n",
    "sum(is.na(geno.dist))  #count NAs in it\n",
    "pca=cmdscale(geno.dist,20, eig=T) #PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#plot PCA\n",
    "#pdf(file=\"pca.pdf\")\n",
    "pca=cmdscale(geno.dist,20,eig=T)\n",
    "per_expl=round(pca$eig/sum(pca$eig)*100,digits=1)\n",
    "plot(pca$points[,1:2], pch=\".\", \n",
    "    xlab=paste(\"PCA axis 1 (\",per_expl[1],\"%)\"),\n",
    "    ylab=paste(\"PCA axis 2 (\",per_expl[2],\"%)\"))\n",
    "text(pca$points[,1:2], label=rownames(pca$points), cex = .6)\n",
    "abline(h=-30,v=-30)\n",
    "#dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#identify outlier species\n",
    "minpc1=-30  #remove samples below this value on pc1\n",
    "minpc2=-30  #remove samples below this value on pc2\n",
    "sp.out=names(which(pca$points[,1]<minpc1 | pca$points[,2]<minpc2))\n",
    "length(sp.out)\n",
    "sp.out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#filter to remove failed, geographic outliers, and species outliers\n",
    "remove_ids=sort(unique(c(failed,geog.out,sp.out)))\n",
    "length(remove_ids)\n",
    "remove_ids\n",
    "\n",
    "#remove samples\n",
    "samples=samples[!samples$SampleID %in% remove_ids,]\n",
    "dim(samples)[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#map remaining samples\n",
    "#pdf(file=\"map.pdf\")\n",
    "oz(ylim=c(min(samples$Latitude)-1,max(samples$Latitude)+1),\n",
    "   xlim=c(min(samples$Longitude)-1, max(samples$Longitude)+1))\n",
    "text(samples$Longitude,samples$Latitude, \n",
    "     label=samples$PopulationName, cex = .6) #add population names\n",
    "#dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "##plot other PCA dimentions\n",
    "#par(mfrow=c(3,3), mar=c(4,4,2,1))\n",
    "#plot(pca$points[,c(1,2)], pch=\"\")\n",
    "#text(pca$points[,c(1,2)],label=rownames(pca$points), cex = .6)\n",
    "#plot(pca$points[,c(1,3)], pch=\"\")\n",
    "#text(pca$points[,c(1,3)], label=rownames(pca$points), cex = .6)\n",
    "#plot(pca$points[,c(1,4)], pch=\"\")\n",
    "#text(pca$points[,c(1,4)], label=rownames(pca$points), cex = .6)\n",
    "#plot(pca$points[,c(2,1)], pch=\"\")\n",
    "#text(pca$points[,c(2,1)], label=rownames(pca$points), cex = .6)\n",
    "#plot(pca$points[,c(2,3)], pch=\"\")\n",
    "#text(pca$points[,c(2,3)], label=rownames(pca$points), cex = .6)\n",
    "#plot(pca$points[,c(2,4)], pch=\"\")\n",
    "#text(pca$points[,c(2,4)], label=rownames(pca$points), cex = .6)\n",
    "#plot(pca$points[,c(3,1)], pch=\"\")\n",
    "#text(pca$points[,c(3,1)], label=rownames(pca$points), cex = .6)\n",
    "#plot(pca$points[,c(3,2)], pch=\"\")\n",
    "#text(pca$points[,c(3,2)], label=rownames(pca$points), cex = .6)\n",
    "#plot(pca$points[,c(3,4)], pch=\"\")\n",
    "#text(pca$points[,c(3,4)], label=rownames(pca$points), cex = .6)\n",
    "\n",
    "##plot dendrogram\n",
    "#hc <- hclust(geno.dist)\n",
    "#plot(hc, cex=.5)\n",
    "\n",
    "##determine outliers\n",
    "#assigns=cutree(hc, k=3)\n",
    "#outliers2=subset(assigns, assigns==2)\n",
    "#outliers3=subset(assigns, assigns==3)\n",
    "#sp=subset(assigns, assigns==1)\n",
    "#length(outliers2)\n",
    "#outliers2\n",
    "#length(outliers3)\n",
    "#outliers3\n",
    "#length(sp)\n",
    "#assigns"
   ]
  }
 ],
 "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
}