Used in analyzing snp calls for Axiom_aegypti1 chip

SNPolish.R

library("SNPolisher")

cwd <- getwd()
summ <- file.path(cwd,"AxiomGT1.summary.txt")
conf <- file.path(cwd,"AxiomGT1.confidences.txt")
post <- file.path(cwd,"AxiomGT1.snp-posteriors.txt")
call <- file.path(cwd,"AxiomGT1.calls.txt")

# Ps_Metrics
ps.metrics <- Ps_Metrics(posteriorFile=post,
  callFile=call,
  output.metricsFile=file.path(cwd,"Output/metrics.txt"))

ps.metrics[1:5,]


# Ps_Classification
Ps_Classification(metrics.file=file.path(cwd,"Output/metrics.txt"), 
  output.dir=file.path(cwd,"Output"), SpeciesType="Diploid", CR.cut=95)

ps.performance <- read.table(file.path(cwd,"Output/Ps.performance.txt"), header=T)


# Ps_Visualization
Ps_Visualization(pidFile=file.path(cwd,"Output/PolyHighResolution.ps"), 
  output.pdfFile=file.path(cwd,"Output/Cluster_PolyHighResolution.pdf"), 
  summaryFile=summ, 
  callFile=call, 
  confidenceFile=conf, 
  posteriorFile=post, 
  temp.dir=file.path(cwd,"Output/Temp"), keep.temp.dir=FALSE, plot.prior=T, 
  match.cel.file.name=TRUE, max.num.SNP.draw=6)

# changing colors
Ps_Visualization(pidFile=file.path(cwd,"Output/PolyHighResolution.ps"), 
  output.pdfFile=file.path(cwd,"Output/Cluster_PolyHighResolution_blue.pdf"), 
  summaryFile=summ, 
  callFile=call,
  confidenceFile=conf, 
  posteriorFile=post, 
  temp.dir=file.path(cwd,"Output/Temp"), keep.temp.dir=FALSE, plot.prior=T, 
  match.cel.file.name=TRUE, max.num.SNP.draw=6,
  geno.col=c("red","dodgerblue2","blue","gray","cyan","green", "darkgreen","purple"))

# highlighting samples, no reference genotypes
# Ps_Visualization(pidFile=file.path(cwd,"Output/PolyHighResolution.ps"),
#   output.pdfFile=file.path(cwd,"Output/Cluster_PolyHighResolution_highlight.pdf"),
#   summaryFile=summ,
#   callFile=call,
#   confidenceFile=conf,
#   posteriorFile=post,
#   temp.dir=file.path(cwd,"Output/Temp"), keep.temp.dir=FALSE, plot.prior=T, 
#   match.cel.file.name=TRUE, max.num.SNP.draw=6, 
#   geno.col=c("red","dodgerblue2","blue","gray","cyan","green", "darkgreen","purple"),
#   sampFile=file.path(cwd,"sample_list.txt")



# OTV_Caller
OTV_Caller(summaryFile=summ,
  posteriorFile=post,
  pidFile=file.path(cwd,"Output/OffTargetVariant.ps"),
  output.dir=file.path(cwd,"OTV/"))

# Ps_Visualization on OTV results
Ps_Visualization(pidFile=file.path(cwd,"OTV/OTV.list.ps"), 
  output.pdfFile=file.path(cwd,"OTV/Cluster_OTV.pdf"), 
  summaryFile=file.path(cwd,"OTV/OTV.summary.txt"), 
  callFile=file.path(cwd,"OTV/OTV.calls.txt"),
  confidenceFile=file.path(cwd,"OTV/OTV.confidences.txt"),
  posteriorFile=file.path(cwd,"OTV/OTV.snp-posteriors.txt"),
  temp.dir=file.path(cwd,"OTV/Temp"), keep.temp.dir=FALSE, plot.prior=T, 
  match.cel.file.name=TRUE, max.num.SNP.draw=6)

# Ps_CallAdjust preparation
Ps_Visualization(pidFile=file.path(cwd,"Output/CallRateBelowThreshold.ps"),
  output.pdfFile=file.path(cwd,"Output/Cluster_Below.pdf"),
    summaryFile=summ,
    callFile=call,
    confidenceFile=conf,
    posteriorFile=post,
    temp.dir=file.path(cwd,"Output/Temp"), keep.temp.dir=FALSE, 
    plot.prior=T, match.cel.file.name=TRUE, max.num.SNP.draw=6)

# Ps_CallAdjust
Ps_CallAdjust(callFile=call, confidenceFile=conf, threshold=0.01,
  outputFile=file.path(cwd,"Call_Adjust/call_adjust_calls.txt"))

# redo analysis with new calls
ps.metrics <- Ps_Metrics(posteriorFile=file.path(cwd,"AxiomGT1.snp-posteriors.txt"),
  callFile=file.path(cwd,"Call_Adjust/call_adjust_calls.txt"),
  output.metricsFile=file.path(cwd,"Call_Adjust/call_adjust_metrics.txt"))

Ps_Classification(metrics.file=file.path(cwd,"Call_Adjust/call_adjust_metrics.txt"),
  output.dir=file.path(cwd,"Call_Adjust"), SpeciesType="Diploid", CR.cut=95)

Ps_Visualization(pidFile=file.path(cwd,"Call_Adjust/CallRateBelowThreshold.ps"),
  output.pdfFile=file.path(cwd,"Call_Adjust/Cluster_Below_new.pdf"),
    summaryFile=summ,
    callFile=file.path(cwd,"Call_Adjust/call_adjust_calls.txt"),
    confidenceFile=conf,
    posteriorFile=post,
    temp.dir=file.path(cwd,"Call_Adjust/Temp"),
    keep.temp.dir=FALSE, plot.prior=T, match.cel.file.name=TRUE, max.num.SNP.draw=6)

# redo plots with new calls
Ps_Visualization(output.pdfFile=file.path(cwd,"Call_Adjust/Cluster_Below_new_2.pdf"),
  summaryFile=summ,
  callFile=file.path(cwd,"Call_Adjust/call_adjust_calls.txt"),
  confidenceFile=conf,
  posteriorFile=post,
  temp.dir=file.path(cwd,"Call_Adjust/Temp"),
  keep.temp.dir=FALSE, plot.prior=T, match.cel.file.name=TRUE, max.num.SNP.draw=6)